(a) Schematic representation of the diffusion-reaction operator-splitting. The final value after diffusion process at time t + Δt 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 d = 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).
(a) Histogram of ln (1/r), where r 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.
A + B → C case study: (a) Initially, species A and B exist only in left-hand side. All A and B molecules and their product P 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.
A + B → C case study: Effect of diffusion constant, D (m 2/s), on (a) τ D (or Δt) and (b) computational time for our method and the GMP method. As D increases, τ D (or Δt) decreases and computational time increases for both algorithms. For D = 10−14 m2/s, the system transitions from diffusion-controlled (Δt = k 2τ D ; k 2 = 2) to reaction-controlled regime during the time-course. For D ⩾ 10−13 m2/s, the system becomes reaction-controlled (Δt = 10τ D ), explaining the increase in the absolute value of the slope of Δt or computational time vs. D plots at D = 10−13 m2/s for our method.
A + B → C case study: (a) and (b) As cell size becomes finer, the results from our approach converge to the numerical solution. L is the number of cells along each direction (larger L represents finer cell size). The vertical axis shows the number of C 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 Δt 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.
Gene expression case study: (a) Dash-dotted line shows the result of Gillespie algorithm which deals with only reaction process. For D = 10−12 m2/s with L = 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 D = 10−15 m2/s with L = 5 exhibits long time lag to reach the steady-state value due to diffusion effect and has much larger fluctuations. (b) In case of L = 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 L is increased from 5 to 20.
Gene expression case study: (a) The result of GMP algorithm for various L and the corresponding Δt (=τ D ) values. Diffusion constant has a fixed value, D = 10−15 m2/s. For L = 20 ( s and τ D = 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 Δt. (b) Our algorithm performs well for both cases of L 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.
CheY diffusion case study: (a) E. coli 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 ).
The tunable parameter k 1 is used as a criterion to decide if the system is diffusion- or reaction-controlled. As the probability of τ R being less than τ D increases, the system becomes more diffusion-controlled. The other parameter, k 2, is related to the probability of a reaction taking place during Δt. As k 2 increases, the probability of a reaction occurrence during Δt increases. In our algorithm, k 1 = 0.5, , k 2 = 2, and are used. Please refer to Fig. 2 .
A synthetic reaction-diffusion system A + B → C with diffusion constant D = 10−12 m2/s. The system is reaction-controlled for L = 4 and 8. In these cases, we set the average value of all Δt to . Computational time increases with smaller Δt or larger L (smaller cell size). For L = 16, the system undergoes transitions between the mixed and diffusion-controlled regimes. In this case, . As L increases, τ D (or Δt) decreases and computational time increases for both algorithms. However, for any L, our method is faster than the GMP method. The relative error-rate is also shown (see Eq. (12) and related text). As L 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 L, as D increases, τ D (or Δt) decreases, and computational time increases for both algorithms.
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 Δt in the GMP method equals τ D , whereas Δt in our method can vary according to the system classification as reaction- or diffusion-controlled. Since the cases of L = 5 and 10 are diffusion-controlled, we set . For L = 20, the system changes from diffusion- to reaction-controlled as time progresses.
Article metrics loading...
Full text loading...