banner image
No data available.
Please log in to see this content.
You have no subscription access to this content.
No metrics data to plot.
The attempt to load metrics for this article has failed.
The attempt to plot a graph for these metrics has failed.
Stochastic operator-splitting method for reaction-diffusion systems
Rent this article for


Image of FIG. 1.
FIG. 1.

(a) Schematic representation of the diffusion-reaction operator-splitting. The final value after diffusion process at time + Δ is used as the initial value for the reaction process. Final value of reaction process is the final value at the end of diffusion-reaction process. (b) and (c) Cellular automata neighborhoods in = 2 dimension: in the von Neumann automata the probability of staying in a cell or diffusing to its neighbors is 1/5 (b), in the Moore automata this probability is 1/9 (c).

Image of FIG. 2.
FIG. 2.

(a) Histogram of ln (1/), where is a uniformly distributed random variable in [0,1]. (b) Cumulative fraction of counts out of total counts (one million). About 63% of the numbers have values less than 1 and 86% of the numbers are less than 2.

Image of FIG. 3.
FIG. 3.

+ case study: (a) Initially, species and exist only in left-hand side. All and molecules and their product diffuse with the same diffusion constant. (b) Comparison of results from analytical solution, cellular automata (CA) and Brownian dynamics (BD). The Brownian dynamics results agree with the analytical solution, while the cellular automata results do not. The increasing curves represent the number of molecules in the right half of the domain and the decreasing ones in the left half.

Image of FIG. 4.
FIG. 4.

+ case study: Effect of diffusion constant, ( 2/), on (a) τ (or Δ) and (b) computational time for our method and the GMP method. As increases, τ (or Δ) decreases and computational time increases for both algorithms. For = 10−14 m2/s, the system transitions from diffusion-controlled (Δ = τ; = 2) to reaction-controlled regime during the time-course. For ⩾ 10−13 m2/s, the system becomes reaction-controlled (Δ = 10τ), explaining the increase in the absolute value of the slope of Δ or computational time vs. plots at = 10−13 m2/s for our method.

Image of FIG. 5.
FIG. 5.

+ case study: (a) and (b) As cell size becomes finer, the results from our approach converge to the numerical solution. is the number of cells along each direction (larger represents finer cell size). The vertical axis shows the number of molecules produced in the left half of the domain. The results are based on average of 8 realizations. (a) shows the results of reaction-controlled system, whereas (b) is for a diffusion-controlled system. (c) and (d) Dashed line is the result of reaction-first and diffusion-later order and dotted line is the reverse order. Diffusion-reaction ordering has better agreement with PDE solution than the reaction-diffusion order. It is because molecules cannot diffuse during the time-step if we treat reaction first. This ordering is effective for both diffusion-controlled and reaction-controlled systems. (e) Black line denotes the result of PDE solution. Three gray lines are for the GMP method. In the GMP method, time-step Δ is related to the cell size. So, it takes longer time to simulate the system with finer cell size. (f) As the initial number of molecules gets lower, the difference rate between deterministic solution and our stochastic solution increases. The fluctuations also increase. This result proves that the stochastic solution approaches the deterministic solution as the number of molecules (or concentration if volume is fixed) increases. In other words, randomness or stochasticity is less important at higher concentrations.

Image of FIG. 6.
FIG. 6.

Gene expression case study: (a) Dash-dotted line shows the result of Gillespie algorithm which deals with only reaction process. For = 10−12 m2/s with = 5, our results are similar to those obtained with the Gillespie algorithm, because it is a reaction-limited process so that diffusion does not have serious impact on the system. On the contrary, the case of = 10−15 m2/s with = 5 exhibits long time lag to reach the steady-state value due to diffusion effect and has much larger fluctuations. (b) In case of = 20, the mesh is much finer than the above cases. The results are similar to those in (a) as our algorithm is able to adjust time-steps according to the system characteristics even though is increased from 5 to 20.

Image of FIG. 7.
FIG. 7.

Gene expression case study: (a) The result of GMP algorithm for various and the corresponding Δ (=τ) values. Diffusion constant has a fixed value, = 10−15 m2/s. For = 20 ( s and τ = 0.417 s), the number of P molecules does not reach its steady-state value of around 1000 because (Table III ). This means reactions cannot fully take place during Δ. (b) Our algorithm performs well for both cases of because it classifies the system as diffusion- or reaction-controlled and decides the appropriate time-steps accordingly. (c) Fluctuations become smaller as the number of realizations become larger. In comparison to the Gillespie algorithm, our method has much higher fluctuations because it considers spatial randomness as well as randomness due to the small number of molecules.

Image of FIG. 8.
FIG. 8.

CheY diffusion case study: (a) has length [2.48 0.88 0.88] μm along x, y, and z direction. R denotes receptor cluster located on the anterior wall and is represented by [x_min x_max; y_min y_max; z_min z_max] = [0 0.08; 0.16 0.64; 0.16 0.64] μm. CheY molecules are phosphorylated in the receptor cluster and diffuse into the cytoplasm. M1∼M4 show the location of the four flagellar motors on side walls and are located in M1 = [0.48 0.56; 0.40 0.48; 0 0.08] μm, M2 = [0.96 1.04; 0 0.08; 0.40 0.48] μm, M3 = [1.44 1.52; 0.40 0.48; 0.80 0.88] μm, and M4 = [1.92 2.00; 0.80 0.88; 0.40 0.48] μm. The remaining domain is considered as cytoplasm. (b) and (c) The simulation results from both the GMP method and our method are in good agreement although there are some differences in rise-time because the GMP method and our method use cellular automata and Brownian dynamics, respectively, to model the diffusion process (as explained in Fig. 3 ).


Generic image for table
Generic image for table
Table I.

The tunable parameter is used as a criterion to decide if the system is diffusion- or reaction-controlled. As the probability of τ being less than τ increases, the system becomes more diffusion-controlled. The other parameter, , is related to the probability of a reaction taking place during Δ. As increases, the probability of a reaction occurrence during Δ increases. In our algorithm, = 0.5, , = 2, and are used. Please refer to Fig. 2 .

Generic image for table
Table II.

A synthetic reaction-diffusion system + with diffusion constant = 10−12 m2/s. The system is reaction-controlled for = 4 and 8. In these cases, we set the average value of all Δ to . Computational time increases with smaller Δ or larger (smaller cell size). For = 16, the system undergoes transitions between the mixed and diffusion-controlled regimes. In this case, . As increases, τ (or Δ) decreases and computational time increases for both algorithms. However, for any , our method is faster than the GMP method. The relative error-rate is also shown (see Eq. (12) and related text). As increases, the relative error-rate decreases for both methods. However, our algorithm is more accurate than the GMP method. In a similar way, for a given , as increases, τ (or Δ) decreases, and computational time increases for both algorithms.

Generic image for table
Table III.

Gene expression case study: Reaction time is averaged over 256 realizations of a simplified gene expression process. As cell sizes become smaller reaction times increase, since fewer molecules in each cell imply lower probability for reactions to take place within a cell. The time-step Δ in the GMP method equals τ, whereas Δ in our method can vary according to the system classification as reaction- or diffusion-controlled. Since the cases of = 5 and 10 are diffusion-controlled, we set . For = 20, the system changes from diffusion- to reaction-controlled as time progresses.


Article metrics loading...


Full text loading...

This is a required field
Please enter a valid email address
752b84549af89a08dbdd7fdb8b9568b5 journal.articlezxybnytfddd
Scitation: Stochastic operator-splitting method for reaction-diffusion systems