Predicting properties of the stationary probability currents for two-species reaction systems without solving the Fokker-Planck equation

We derive methods for estimating the topology of the stationary probability current ${\vec{j}}_s$ of the two-species Fokker-Planck equation (FPE) without the need to solve the FPE. These methods are chosen such that they become exact in certain limits, such as infinite system size or vanishing coupling between species in the diffusion matrix. The methods make predictions about the fixed points of ${\vec{j}}_s$ and their relation to extrema of the stationary probability distribution and to fixed points of the convective field, which is related to the deterministic drift of the system. Furthermore, they predict the rotation sense of ${\vec{j}}_s$ around extrema of the stationary probability distribution. Even though these methods cannot be proven to be valid away from extrema, the boundary lines between regions with different rotation senses are obtained with surprising accuracy. We illustrate and test these methods, using simple reaction systems with only one coupling term between the two species as well as a few generic reaction networks taken from the literature. We use it also to investigate the shape of nonphysical probability currents occurring in reaction systems with detailed balance due to the approximations involved in deriving the Fokker-Planck equation.

Figure 1

The different contributions to ${\vec{j}}_s$ for the rotationally symmetric system, (51), for the parameter values $\omega =\lambda =d=N=1$: (a) stationary probability density $p_s(x,y)$; (b) stationary probability current ${\vec{j}}_s(x,y)$; (c) curl of the stationary probability current $|\vec{\nabla }\times {\vec{j}}_s(x,y)|$ (as all vectors point along the positive $z$ axis, only the absolute value is shown); (d) cut along the $x$ axes in (c).

Figure 2

Stationary probability current ${\vec{j}}_s$ for the coupled system (58) with (63) and (61) for the parameter set $r_1=r_2=50, s_1=s_2=1,d=100,N=1$, and different values for the coupling terms: (a) $a_1=0.1,a_2=0.2,b_i=c_i=0$; (b) $a_i=b_i=c_2=0,c_1=0.005$; (c) $a_i=c_i=0,b_1=0.001,b_2=0.002$; (d) $a_i=c_i=0,b_1=b_2=0.001$. The green (red) background color indicates the predicted positive (negative) rotation sense as of Eq. (59), i.e., based on the sign of the curl of $\vec{\alpha }$. Red dots indicate the positions of the fixed points of $\vec{\alpha }$; white rectangles, the fixed points of ${\vec{j}}_s$.

Figure 3

Stationary probability current for the example system, (63), with linear coupling in the diagonal terms of the diffusion, (72), with $D_{11}^u\left(x\right)=D_{22}^u\left(y\right)=d=100$ and different coupling strengths: (a) $a_1=1,a_2=0$; (b) $a_1=1,a_2=5$. The green (red) background color indicates the positive (negative) rotation sense predicted based on the sign of Eq. (73). Red dots indicate the positions of the fixed points of $\vec{\alpha }$; white rectangles, the fixed points of ${\vec{j}}_s$.

Figure 4

Stationary probability current (blue arrows) for system (77) with (a) symmetric and (b) asymmetric parameter values. The green (red) background color indicates the predicted positive (negative) sense of rotation according to Eq. (39). Red dots indicate the positions of the fixed points of $\vec{\alpha }$; white rectangles, the fixed points of ${\vec{j}}_s$. Parameter values: (a) $m_x=m_y=1,n_x=n_y=4,{\mathrm{\Theta }}_x={\mathrm{\Theta }}_y=4,d_x=d_y=0.2,b_x=b_y=0.34,N=20$; (b) $m_x=m_y=1,n_x=4,n_y=3,{\mathrm{\Theta }}_x={\mathrm{\Theta }}_y=4,d_x=0.12,d_y=0.22,b_x=0.22,b_y=0.15,N=20$.

Figure 5

Color gradient: height profile of $p_s$. Arrows: direction of ${\vec{j}}_s$. (a) A saddle point in $p_s$ between two maxima with equal rotation sense in ${\vec{j}}_s$ leads to the formation of a saddle in ${\vec{j}}_s$ at the location of the $p_s$ saddle. (b) If a $p_s$ saddle lies between two maxima with opposing rotation senses in ${\vec{j}}_s$, the formation of a ${\vec{j}}_s$ saddle is prevented.

Figure 6

Stationary probability density (density plot in blue and yellow) and current (vector/streamline plot in red) for system (81) with parameter sets that (a) create a stable fixed point or (b) lead to an unstable fixed point enclosed by a limit cycle. The darker-red vectors in (b) indicate the direction and strength of the vector field $\vec{j}$, whereas the lighter-red streamlines indicate only its direction in places where its magnitude is too weak to use the former visualization. Furthermore, the red dots in both panels indicate the fixed point in $\vec{\alpha }$, the white rectangle in (b) marks the minimum of the probability density, and the black line in (b) shows the location of the limit cycle in $\vec{\alpha }$. Parameters: (a) $A=2/3,\beta =0.8,d=0.65,c=0.65,b=1,\delta =0.2,q=0.01,\mathrm{\Omega }=1000$; (b) $A=2/3,\beta =0.9,d=0.65,c=0.65,b=1,\delta =0.1,q=0.01,\mathrm{\Omega }=1000$.

Figure 7

Estimated rotation sense of ${\vec{j}}_s$ based on the sign of $\vec{\nabla }\times {\mathbf{D}}^{-1}\vec{\alpha }$ (background color: green, positive; red, negative) and numerical results for ${\vec{j}}_s$ (blue vectors). Red dots indicate the fixed points of the convective field which coincide with the corresponding maxima of the probability density. (a) Model (85), with parameters chosen as in [20], i.e., in our notation $k_1=1,k_{-1}=3,k_2=2,k_{-2}=500,k_3=75,k_{-3}=0.1,n=2\times {10}^4$. (b) Model (90), with parameters $r=50,s=\lambda =1,N=1$.