† Corresponding author. E-mail:

Project supported by the Thailand Research Fund and Mahidol University (Grant No. TRG5880157), the Thailand Center of Excellence in Physics (ThEP), CHE, Thailand, and the Development Promotion of Science and Technology.

Most biochemical processes in cells are usually modeled by reaction–diffusion (RD) equations. In these RD models, the diffusive process is assumed to be Gaussian. However, a growing number of studies have noted that intracellular diffusion is anomalous at some or all times, which may result from a crowded environment and chemical kinetics. This work aims to computationally study the effects of chemical reactions on the diffusive dynamics of RD systems by using both stochastic and deterministic algorithms. Numerical method to estimate the mean-square displacement (MSD) from a deterministic algorithm is also investigated. Our computational results show that anomalous diffusion can be solely due to chemical reactions. The chemical reactions alone can cause anomalous sub-diffusion in the RD system at some or all times. The time-dependent anomalous diffusion exponent is found to depend on many parameters, including chemical reaction rates, reaction orders, and chemical concentrations.

Diffusion is likely to be the most important transportation process inside a living cell. However, the cell cytoplasm is not a simple solvent; it is crowded and consists of many chemically reactive molecules. Therefore, various factors can hinder diffusion in living cells.^{[1,2]} In particular, the chemical reactions or biochemical conversions in metabolic and signal transduction pathways can affect the cellular transportation processes.^{[1,3]} In general, the reaction diffusion processes in a living cell can be described by a set of reaction–diffusion equations,

*c*

_{i}(

*,*r

*t*),

*D*

_{i}, and

*R*

_{i}are the concentration, diffusion constant, and reaction term of the

*i*-th substance, respectively. For systems in which particles or molecules diffuse freely via Brownian dynamics, the mean-square displacement (MSD) of the particles increases linearly with time. Some researchers also hold that this reaction diffusion process can be investigated under the assumption of Brownian or normal diffusion.

^{[4,5]}However, a number of recent experimental studies based on single-particle tracking (SPT) and fluorescence correlation spectroscopy (FCS) mapped out the intracellular molecules’ trajectories to determine the relationship of MSD with time and revealed that the diffusions in cells and on cell membranes are often anomalous.

^{[6–9]}This was initially suggested by the experimental results of Cajal bodies in the nucleus of HeLa cells.

^{[10]}Two other notable studies of the discovery of anomalous sub-diffusion in living cells were the diffusion of telomeres in the nucleus of mammalian cells and the diffusion of protein in HeLa cells and polytene cells of Drosophila larval salivary glands.

^{[11,12]}It was reported that the diffusions in different cell types cause different anomalous diffusions and the anomalous behavior of intracellular diffusion can be presumably the result of the obstruction of molecular crowding in the cytoplasm as well as chemical kinetics.

^{[1,3,6,13,14]}However, the experimental data on diffusion in living cells are difficult to identify what factors exactly cause the anomalous diffusion within cells. Therefore, theoretical and simulation studies are necessarily required to add more understanding of diffusion and transport mechanisms of biological processes in living cells.

Anomalous diffusion is characterized by the MSD of diffusing particles, which grows with the power law of time:

*α*represents an anomalous diffusion exponent that indicates the diffusion behavior.

^{[15]}If

*α*> 1, the diffusion process is characterized by super-diffusion, whereas

*α*=1 corresponds to normal diffusion.

^{[1,13,16]}Anomalous sub-diffusion occurs when 0 <

*α*< 1 and has been found in many intracellular systems.

^{[3,6]}However, the underlying physical or chemical causes of anomalous sub-diffusion in cells are not fully understood.

^{[1,3]}

The ability of chemical reactions to cause anomalous sub-diffusion has not yet been experimentally investigated. However, it has been observed and clearly identified in many other cases. Many researches demonstrated that the continuous time random walk (CTRW) model is one of the most successful theoretical models of anomalous diffusion.^{[13,17,18]} Also theoretical models based on the CTRW and fractional reaction–diffusion equations confirmed that chemical components can cause anomalous sub-diffusion.^{[13,17,19–23]} Likewise, Haugh used the framework of Green’s function analysis to demonstrate that the reaction diffusion mechanism could not be assumed to be normal diffusion.^{[3]} Using Monte Carlo simulations on a simple lattice model, Saxton also showed that the binding mechanism alone can cause anomalous diffusion at some or all times.^{[24]} In addition, by using lattice gas automata, researchers have recently found that the diffusion in the Michaelis-Menten mechanism exhibits anomalous diffusion in a transient period.^{[25,26]}

In this work, we further investigate the diffusive anomaly of a reaction diffusion system by using MCell, a 3D spatially realistic Monte Carlo simulator that tracks individual molecules.^{[27]} Although the MCell program has been optimized for Monte Carlo simulation, a deterministic approach is computationally more efficient and easier to implement for most problems. Thus, we propose alternative methods using a deterministic algorithm to investigate the anomalous behaviors in reaction diffusion systems which have been numerously elucidated by stochastic simulation in published versions. The effects of reaction rate and reaction type on the diffusion behavior are also investigated.

Reaction–diffusion equations are a mathematical model that explains the changes of one or more substances in space and time in chemical reaction and diffusion processes. In addition to the study of chemistry, these equations can also be used to describe many intracellular processes.^{[28,29]} In the present work, we study a *Ca*^{2+} buffering system in cells. This *Ca*^{2+} buffering system can control various processes, including neurotransmitter release during synaptic transmission. The calcium buffering dynamics can be written in the form of a chemical equation as follows:

*B*

_{i}and

*CaB*

_{i}represent the free buffers and calcium bound buffers, respectively, and

*i*=

*s*or

*m*represents stationary or mobile buffer, respectively.

^{[5]}In our model, an initial number of

*Ca*

^{2+}ions are released in the center of an unbounded space, where either mobile or stationary buffers are homogeneously distributed.

We investigate anomalous diffusion in the *Ca*^{2+} buffering system by using MCell version 3.1, a Monte Carlo (MC) algorithm that tracks individual molecules.^{[27]} The MCell simulations are performed in three-dimensional (3D) space, and all of the *Ca*^{2+} trajectories are recorded during the simulations. The recorded trajectories are then used to calculate the mean squared displacement (MSD) of *Ca*^{2+}. In all simulations, the dimensions of the computational geometry are 4 μm × 4 μm × 4 μm. Initially, 54902 ions of *Ca*^{2+} are held in a spherical shell with a radius of 7.5 nm located at the origin and are then released at time *t* = 0, whereas 549020 ions of buffers are uniformly distributed throughout the computational box. These numbers of ions correspond to a calcium concentration of 1.425 μM and a buffer concentration of 14.25 μM in our computational domain. All surfaces of the computational box are reflective to all molecules. To study only the reaction diffusion dynamics and avoid the possible effects of the boundary on calcium ion dynamics, we ensure that calcium ions and calcium bound buffers do not reach the boundaries of our computational box. All parameter values used in MCell simulations are adjusted to reduce the simulation time, as given in Table

Deterministic simulation is an alternative approach to the study of the diffusive behaviors in the system in order to obtain longer time courses of the MSD. Assuming mass action kinetics and Fickian diffusion, the concentration changes in Eq. (

*D*

_{Ca},

*D*

_{Bi}, and

*D*

_{CaBi}are the diffusion constants of free calcium, free buffers, and

*Ca*

^{2+}bound buffers, respectively, and

*i*, respectively. In deterministic simulations, we treat the computational space as a visual one-dimensional infinite system, within which the initial pulse of

*Ca*

^{2+}is released as the Dirac delta function ([

*Ca*

^{2+}]

_{0}=

*δ*(

*x*)) and buffers initially distributed homogeneously. The reaction diffusion equations in Eq. (

The deterministic calculations that solve the differential equations in Eq. (^{2+} is estimated from either the half width at half maximum (HWHM) or numerical integration of the PDFs. Both MSD estimation methods are generic and can be used with any PDFs. In the first method, the shape of the PDF is characterized by the HWHM value. This approach has previously been used to investigate anomalous diffusion as described in Sagi, *et al.*’s work.^{[31]} The HWHM is estimated in MATLAB using the linear interpolation function, interp1, which requires at least two known data points. The syntax of interp1 is interp1 (*x*, *y*, *x*_{i}), where *x* is a vector that contains independent variable data, *y* is a vector that contains dependent variable data and *x*_{i} is the array of data points to be interpolated.

For the integration method, the MSD (〈*x*^{2}〉) can be directly calculated from definition,^{[32]}

*P*(

*x,t*) is a PDF. The equation is numerically integrated using the trapezoidal rule.

Different types of reactions are investigated to systematically study the effects of reactions in reaction diffusion systems.

First, we used MCell, a particle-based Monte Carlo algorithm, to investigate anomalous diffusion in reversible bimolecular reaction, which can be found in most living cells. The particle trajectories in three-dimensional space are recorded and used to calculate their mean squared displacement (MSD or 〈*r*^{2}〉) as shown in Figs. *r*^{2}〉 as a function of time, obtained from the free diffusion system and a reaction diffusion system that exhibits anomalous diffusion. The distinction between normal diffusion and anomalous diffusion becomes more apparent when the data are plotted as log (〈*r*^{2}〉/*t*) versus log*t* as shown in Fig. *α* − 1, normal diffusion produces a horizontal line, which corresponds to *α* = 1. The result shows that the anomalous diffusion exponent is not constant (see more results in Section 3). To avoid boundary effects, the box size must be sufficiently large to prevent calcium ions and calcium-bound buffers from reaching the surfaces of the box. This constraint causes the simulations to require a very long time to reach the steady state. Obtaining 10 s of the simulation result in Fig.

We cannot track the movement of individual molecules in the deterministic simulation; thus, the MSD cannot be directly computed from particle trajectories, and an alternative method to compute the MSD becomes necessary. We propose two alternative approaches to extract the MSD from the PDFs. In the deterministic calculations, we obtain PDFs at different times (*P*(*x*,*t*)) by solving the differential equations in Eq. (

For the free diffusion in an infinite system with an initial [*Ca*^{2+}] given by the Dirac delta function, the PDF can be written in the form of a Gaussian function^{[34]} as verified in Fig. *f* (*x*) = *A* exp(−*x*^{2}/2*σ*^{2}), where *A* is the amplitude and *σ* is the standard deviation,^{[4]} can be utilized to determine the relationship between the MSD and the HWHM because *f* (*HWHM*) = *A* exp(−*HWHM*^{2}/2*σ*^{2}) = *A*/2. This calculation yields *x*^{2}〉, is the variance of the Gaussian distribution, *σ*^{2}, which grows linearly with time, i.e., *σ*^{2} = 2*Dt*.^{[34]} Thus, *HWHM*^{2} = 2 ln 2〈*x*^{2}〉. In *n*-dimensional space, this relationship can also be generalized to the following expression:

^{[13,18–20]}Hence, its general form has not yet been found, although many approaches have been introduced to elucidate the PDFs of anomalous dynamics.

^{[1,18,35]}An explicit relationship between the MSD and the HWHM is lacking for this case. Therefore, in this study Mote Carlo simulations are used to numerically determine the relationship between the HWHM and the MSD for the reaction diffusion process. The MSD is directly calculated from individual particle trajectories, and the HWHM is estimated from the corresponding [

*Ca*

^{2+}] curves. We find that the MSD of a reaction diffusion system remains linearly proportional to HWHM

^{2}as shown in Fig.

*HWHM*

^{2}∝

*t*

^{α}.

The MSD values obtained from HWHM interpolation and numerical integration for normal and anomalous diffusion are compared in Figs.

The reaction diffusion dynamics is henceforth investigated in one-dimensional space by using deterministic simulations in MATLAB. To study the anomalous behaviors of the system, the anomalous diffusion exponent (*α*) in the equation 〈*x*^{2}〉 = *Γ**t*^{α} is computed. As described above, the numerical integration method is selected to investigate the dynamics of reaction diffusion mechanisms and analyze the anomalous diffusion exponents (*α*). In this section, we focus on the reversible reaction of the *Ca*^{2+} buffering system. We investigate the effects of buffer concentration and diffusion coefficient on the diffusive behavior of *Ca*^{2+} in a buffer solution; the parameter values for this system are shown in Table

Figure *x*^{2}〉/*t*) versus log*t* for different buffer concentrations and diffusion coefficients. Because the slope of the straight line in this figure is *α* − 1, normal diffusion yields a horizontal line. This figure indicates that the diffusion of reversible bimolecular reaction is characterized by three time regimes. The diffusion is normal at early times and then rapidly changes to anomalous diffusion at an intermediate time. The diffusion subsequently becomes normalized when the dynamics reaches a steady state as illustrated in Figs.

Because the MSD can be expressed as

*Γ*is a constant, the value of

*α*can be obtained from the slope of the plot of log(〈

*x*

^{2}〉/

*t*) versus log

*t*. The results in Fig.

*α*is not constant over time. Therefore, we calculate the instantaneous anomalous exponent (

*α*

_{ins}), an anomalous diffusion exponent during an infinitesimal time interval as shown in Fig.

*α*= 1 at early times to

*α*< 1 at intermediate times and a return to

*α*= 1 again at late times, are clearly evident, which suggests that the diffusion in reversible bimolecular reactions follows transient anomalous sub-diffusion. Furthermore, figure

Figure *Ca*^{2+} diffusion in mobile buffers and different dissociation rates (*k*^{−}) and constant association rates (*k*^{+}). The values of other parameters used in this simulation are shown in Table *k*^{−}) decreases. The overall dynamics continues to follow the transient anomalous sub-diffusion, similar to the results in Section 3. At high values of *t* when the reaction reaches a steady state, the diffusion reverts to the normal diffusion. However, if *k*^{−} approaches to zero, the reaction becomes an irreversible second-order reaction, *r* = *k*^{+}[*Ca*^{2+}][*B*_{m}]. In this case, diffusion remains normal until the *Ca*^{2+} concentration is extremely close to zero at very high values of *t* as shown in Figs. *B*_{m}]_{0} ≫ [*Ca*^{2+}]_{0}; thus, the concentration of mobile buffers can be considered to be constant and *r* = *k*^{+}[*Ca*^{2+}][*B*_{m}] = *k*′[*Ca*^{2+}], where *k*′ = *k*^{+}[*B*_{m}] is constant. Therefore, the second-order reaction reduces to a pseudo-first order reaction in this case, which results in normal diffusion.^{[31]}

To demonstrate the effect of the initial *Ca*^{2+} concentration, we also change [*Ca*^{2+}]_{0} to 0.623 μM. The results in Fig. *k*− and the slope in the anomalous region inversely correlate. Generally, reducing the dissociation rate exacerbates the anomalous nature of the diffusion. The diffusion remains anomalous even when *k*^{−} = 0, unlike the case when [*Ca*^{2+}]_{0} = 6.23 nM (see Fig. *k*^{−} = 0 result in an anomalous sub-diffusion throughout the simulation time, even when the *Ca*^{2+} concentration reaches very low values at very late times as shown in Fig.

For the first-order reaction

*k*

^{+}and

*D*

_{Ca}are constants. In this case, we can analytically show that the diffusion dynamics of

*Ca*

^{2+}is normal at all times. Now we write the [

*Ca*

^{2+}] in the following form:

*C*

_{0}(

*x*,

*t*) is a function. Substituting Eq. (

*C*

_{0}(

*x*,

*t*),

*C*

_{0}(

*x*,

*t*) is a solution of the diffusion equation. Therefore, the corresponding solution of the first-order reaction diffusion equation with the initial condition

*C*(

*x*,0) =

*δ*(

*x*) is in the form of the Gaussian function multiplied by the exponential function of time. At a particular time, this exponential function exp(−

*k*

^{+}

*t*) does not change the shape of the Gaussian function, except its height. Thus the MSD of

*Ca*

^{2+}in the first-order reaction diffusion process is not affected by the reaction. Therefore

*Ca*

^{2+}shows normal diffusion at all times. This analytical result agrees well with our computational result shown in Fig.

*α*

_{ins}) for different reaction types are summarized in Figs.

To calculate the average anomalous diffusion exponent (*α*_{avr}) in the period of anomalous sub-diffusion, we first find crossovers between normal and anomalous diffusion by plotting log(〈*x*^{2}〉/*t*) versus log(*t*). As shown in previous sections, in the case of a reversible second-order reaction, there are two major transitions at short time and at long time respectively. At the first crossover (*x*_{1},*y*_{1}) diffusion changes from normal to sub-diffusion and at the second crossover (*x*_{2},*y*_{2}) diffusion changes from sub-diffusion to normal diffusion again. To find these crossovers, we locate the intersections between two straight lines in anomalous and normal regions at early time and at long time as shown in the upper inset of Fig. *x*_{1},*y*_{1}) = (log*t*_{cr1},log*D*_{Ca}) and (*x*_{2},*y*_{2}) = (log*t*_{cr2},log*D*(∞)) where *t*_{cr1} and *t*_{cr2} are the first and second crossover time, respectively, and *D*(∞) is an effective diffusion coefficient of calcium at very long time.^{[6]} Figure *t*_{cr1} and *t*_{cr2} each as a function of buffer concentration for various calcium concentrations. It indicates that *t*_{cr1} and *t*_{cr2} are insensitive to the change of the initial *Ca*^{2+} concentration. This also corresponds with the result of *α*_{ins} with [*Ca*^{2+}] varying as shown in the lower inset of Fig.

To get an insight into *t*_{cr1} and *t*_{cr2}, we calculate the half-life of calcium concentration in a well-mixed condition where both calcium and buffer are uniformly distributed. For the reversible bimolecular reactions,

*Ca*

^{2+}(

*t*

_{hl}) is found to be as follows (see Appendix A, for detailed description of the derivation):

*k*=

*k*

^{−}/

*k*

^{+}; [

*Ca*

^{2+}]

_{0}; and [

*B*]

_{0}are the initial concentrations of calcium and buffer, respectively. The calcium half-life calculated from Eq. (

*Ca*

^{2+}concentration (Fig.

*t*

_{cr1}and

*t*

_{cr2}) in Fig.

*t*

_{cr1}is relatively close to the

*t*

_{hl}as illustrated in Fig.

*t*

_{cr1},

*t*

_{cr2}and

*t*

_{hl}are all involved with the buffer concentrations in power-law relations. We also find that

*t*

_{cr1},

*t*

_{cr2}, and

*t*

_{hl}are all insensitive to the change of dissociation rate

*k*

^{−}(Fig.

Likewise, we obtain the calcium effective diffusion coefficient *D*_{eff} (or *D*(∞)), from the simulation results as plotted in Fig. *α*_{avr} in the anomalous sub-diffusive region can be then calculated from the slope of the graph of log(〈*x*^{2}〉/*t*) versus log*t*. Consequently, we obtain the empirical form of *α*_{avr} as a function of buffer concentration as follows:

*α*

_{avr}in the period of anomalous sub-diffusion are illustrated in Fig

In this work, the roles of biochemical reactions in diffusion dynamics are computationally investigated by using both deterministic and stochastic algorithms. Because we cannot track individual particle trajectories in the deterministic simulations, methods to estimate the MSD from deterministic calculations are also investigated. We find that the numerical integration of the PDF best estimates the MSD. Computational and theoretical analyses show that the first-order reaction diffusion dynamics are normal at all times, whereas the second-order reaction diffusion dynamics is more complex and depends on reaction parameters. However, the second-order dynamics is found to follow the anomalous sub-diffusion at least at some times. This work clearly confirms that anomalous diffusion can be a result of only chemical reactions and only second-order reactions can give rise to anomalous diffusion dynamics in a reaction diffusion system.

**Reference**

1 | Biophys. 95 2049 |

2 | Biophys. 107 2579 |

3 | Biophys. 97 435 |

4 | Neuron 18 473 |

5 | Biophys. 67 447 |

6 | Biophys. 92 1178 |

7 | Biophys. 104 1652 |

8 | Chem. Phys. 284 389 |

9 | Magn. Reson. Imaging 31 359 |

10 | Nat. Cell Biol. 4 502 |

11 | Phys. Rev. Lett. 103 018102 |

12 | J. Phys. Chem. 117 1241 |

13 | |

14 | |

15 | Chin. Phys. 22 038701 |

16 | Chin. Phys. 20 040502 |

17 | Chin. Phys. 23 110206 |

18 | Phys. Rep. 339 1 |

19 | Phys. Rev. 73 031102 |

20 | Phys. Rev. 90 032136 |

21 | Chem. Phys. 344 90 |

22 | Chem. Phys. 284 169 |

23 | J. Comput. Appl. Math. 275 216 |

24 | Biophys. 70 1250 |

25 | Phys. Chem. Chem. Phys. 16 4492 |

26 | Biophys. 105 2064 |

27 | SIAM. J. Sci. Comput. 30 3126 |

28 | Angew. Chem. Int. Ed. 49 4170 |

29 | Chin. Phys. Lett. 27 048201 |

30 | Phys. Biol. 7 026008 |

31 | |

32 | Magn. Reson. Med. 63 1323 |

33 | J. Chem. Phys. 135 174104 |

34 | J. R. Soc. Interf. 5 813 |

35 | J. Phys.: Condens. Matter 19 065118 |