Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                
  • Open Access

Reaction-diffusion kinetics on lattice at the microscopic scale

Wei-Xiang Chew, Kazunari Kaizu, Masaki Watabe, Sithi V. Muniandy, Koichi Takahashi, and Satya N. V. Arjunan
Phys. Rev. E 98, 032418 – Published 28 September 2018

Abstract

Lattice-based stochastic simulators are commonly used to study biological reaction-diffusion processes. Some of these schemes that are based on the reaction-diffusion master equation (RDME) can simulate for extended spatial and temporal scales but cannot directly account for the microscopic effects in the cell such as volume exclusion and diffusion-influenced reactions. Nonetheless, schemes based on the high-resolution microscopic lattice method (MLM) can directly simulate these effects by representing each finite-sized molecule explicitly as a random walker on fine lattice voxels. The theory and consistency of MLM in simulating diffusion-influenced reactions have not been clarified in detail. Here, we examine MLM in solving diffusion-influenced reactions in three-dimensional space by employing the spatiocyte simulation scheme. Applying the random walk theory, we construct the general theoretical framework underlying the method and obtain analytical expressions for the total rebinding probability and the effective reaction rate. By matching Collins-Kimball and lattice-based rate constants, we obtained the exact expressions to determine the reaction acceptance probability and voxel size. We found that the size of voxel should be about 2% larger than the molecule. The theoretical framework of MLM is validated by numerical simulations, showing good agreement with the off-lattice particle-based method, enhanced Green's function reaction dynamics (egfrd). MLM run time is more than an order of magnitude faster than egfrd when diffusing macromolecules with typical concentrations observed in the cell. MLM also showed good agreements with egfrd and mean-field models in case studies of two basic motifs of intracellular signaling, the protein production-degradation process and the dual phosphorylation-dephosphorylation cycle. In addition, when a reaction compartment is populated with volume-excluding obstacles, MLM captures the nonclassical reaction kinetics caused by anomalous diffusion of reacting molecules.

  • Figure
  • Figure
  • Figure
  • Figure
  • Figure
  • Figure
  • Figure
1 More
  • Received 4 June 2018
  • Corrected 21 February 2019

DOI:https://doi.org/10.1103/PhysRevE.98.032418

Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.

Published by the American Physical Society

Physics Subject Headings (PhySH)

Physics of Living SystemsInterdisciplinary Physics

Corrections

21 February 2019

Correction: The JSPS contract number contained an error and has been fixed.

Authors & Affiliations

Wei-Xiang Chew1,2, Kazunari Kaizu1, Masaki Watabe1, Sithi V. Muniandy2, Koichi Takahashi1, and Satya N. V. Arjunan1,*

  • 1Laboratory for Biologically Inspired Computing, RIKEN Center for Biosystems Dynamics Research, Suita, Osaka, Japan
  • 2Department of Physics, Faculty of Science, University of Malaya, 50603, Kuala Lumpur, Malaysia

  • *satya@riken.jp

Article Text

Click to Expand

References

Click to Expand
Issue

Vol. 98, Iss. 3 — September 2018

Reuse & Permissions
Author publication services for translation and copyediting assistance advertisement

Authorization Required


×

Images

  • Figure 1
    Figure 1

    A voxel on the hcp lattice has 12 nearest neighbor voxels. The distance between the centers of two adjacent voxels is the voxel size l.

    Reuse & Permissions
  • Figure 2
    Figure 2

    The rebinding time of a reactive pair that is initially in contact. The rebinding time is sampled from simulations with ka/kD = 0.1, 1, and 100. Markers show the simulation results of spatiocyte while solid lines depict the analytical results from the continuum-based scheme (7). The vertical dashed line marks the characteristic diffusion time step td. Simulation parameters: l=0.01μm, volume = (10000l)3 with periodic boundary, runs = 104, DA=1μm2s1, DB=0, α=1/Pa for the diffusion-influenced case.

    Reuse & Permissions
  • Figure 3
    Figure 3

    Survival probability and time-dependent rate coefficient. (a) Survival probability of A in A+BB with ka/kD = 0.1, 1, and 100. (b) Simulated time-dependent rate coefficients of the reaction and the corresponding long-time approximation of Collins-Kimball (CK) theory in Eq. (4). Simulation parameters: volume = (3.5μm)3 with periodic boundary, R=0.01μm, l=0.01×1.0209μm, DA=0, DB=1μm2s1, Na=Nb=4000, duration = 0.05 s, runs = 3×104, α=1/Pa for the diffusion-influenced case.

    Reuse & Permissions
  • Figure 4
    Figure 4

    3D diffusion performance of particle-based methods. Vertical axis T shows the run times to diffuse molecules with diffusion coefficient Dx = 1 μm2s1 in volume V for 10 s. Bottom horizontal axis N represents the number of diffusing molecules, while the top axis shows the corresponding concentration at V=30μm3. (a) Molecules are represented as hard-sphere particles with radius r=l/2 = 2.5 nm. (b) Molecules are dimensionless point particles that can overlap one another. egfrd does not support point particle diffusion and, conversely, fast BD here can only diffuse point particles. smoldyn simulation interval is set to the step interval td (4.17 μs) of spatiocyte and fast BD for comparison. The egfrd algorithm uses variable time steps. Each model was simulated for a predefined run time tr and the resulting simulated time ts was recorded. We calculated T, the run time in seconds for 10 s of simulated time with T=10tr/ts. tr was set such that at least hundreds of simulation steps have been completed. The resulting range of tr was between 1 hour to several days. Solid lines depict the ideal scaling for spatiocyte. Vertical dashed lines indicate the typical concentration range of proteins in the cytoplasm (0.1 to 10 μM). All simulations were executed on the same server with Intel Xeon Platinum 8180 2.5 GHz (max 3.80 GHz) CPU, 768 GB memory, and Ubuntu 18.04 LTS operating system.

    Reuse & Permissions
  • Figure 5
    Figure 5

    Particle simulation performance of the Michaelis-Menten reaction. Original benchmark model from [66, 68] was modified with volume (V) = 90.9 μm3, diffusion coefficient (Dx) = 1 μm2s1, k1 = 0.01 μm3s1, k2 = k3 = 0.1 s1. Molecule or voxel radius (r), simulation or diffusion step interval (Δt), and run time (T) are as indicated. All simulations were executed on the same workstation with Intel Xeon X5680 3.33 GHz CPU, 48 GB memory, and Ubuntu 16.04 LTS operating system.

    Reuse & Permissions
  • Figure 6
    Figure 6

    Production-degradation response of A. (a) Time-series profile of A in Eq. (36) simulated with spatiocyte (using intrinsic rate ka), egfrd, and RDME. DA=DB=0.1μm2s1, molecule radius rA=rB=r{0.005,0.02,0.05}μm. Note that r represents half of the subvolume size in RDME, and the actual molecule radius in spatiocyte and egfrd. For comparison, solid line shows the well-mixed model. (b) Mean equilibrium concentration of A from spatiocyte, egfrd, RDME, and RDME with modified propensity (RDMEm) with DA=DB=Dx{0.01,0.02,0.1}μm2s1. Solid and dashed lines represent expected results according to the well-mixed model and the microscopic theory, respectively. (c) Steady-state distribution of A from spatiocyte and egfrd with r=0.05μm and DA=DB=Dx{0.1,0.02}μm2s1. RDMEm simulated with r=1 and D=0.02 is also shown for comparison. The frequency is normalized such that the sum over the bin is unity. Dotted line represents the well-mixed model simulated using the Gillespie method. Simulation parameters: k1=0.1μm3s1, k2=0.02μm3s1, [B]=1μm3, runs = 700, duration >104 s to achieve steady state, volume = 100μm3 with periodic boundary.

    Reuse & Permissions
  • Figure 7
    Figure 7

    Effects of rebinding in dual phosphorylation cycle. (a) Reaction model showing MAPK (K) is first activated into Kp and then Kpp by MAPKK (KK) in two phosphorylation steps. Kpp is also deactivated by phosphatase (P) in two dephosphorylation steps to become K again. Enzymes KK and P become inactive immediately after reacting with their respective substrates and then relax back to the active state after some delay τrel. (b) Fraction of Kpp in response to MAPKK/phosphatase ratio at steady state. Circle and square markers denote simulation result using spatiocyte with Dx=4μm2s1 and Dx=0.06μm2s1, respectively. Dashed and solid lines represent distributive and processive mechanism models, respectively. Cross and plus markers show the results from the original spatiocyte scheme, wherein the voxel and molecule sizes are exactly the same. We used a short reactivation time τrel=1μs, relative to td (for comparison td1μs when Dx=4μm2s1, td70μs when Dx=0.06μm2s1) with the total number of substrates Ktotal=120. Hysteresis responses from mean-field distributive model with fivefold substrate concentration (Ktotal=600) are indicated by dotted and dashed-dotted lines with initial conditions [Kpp]/[K]total=1 and [Kpp]/[K]total=0, respectively. Diamond and triangle markers represent spatiocyte responses with fivefold substrate concentration at the indicated diffusion coefficient Dx. Simulation parameters: molecule size l=0.0025×1.0209μm, diffusion coefficient Dx, [KK]+[P]=60, duration = 200 s, volume = 1μm3 with periodic boundary.

    Reuse & Permissions
  • Figure 8
    Figure 8

    Diffusion and bimolecular reaction kinetics in crowded compartment. (a) Time-dependent diffusion coefficient of tracer molecules in the presence of immobile obstacles at volume occupancy ϕ. The diffusion coefficient at a time point is determined from the mean-squared displacement of simulated particle trajectories. Dashed lines denote the diffusion coefficient at long time as predicted by D=D0(1ϕ/ϕp). (b) Survival probability of E in E+S at ϕ. (c) The corresponding time-dependent reaction rates (dashed lines) at ϕ. Solid lines represent Collins-Kimball theory with the long-time diffusion coefficient calculated in (a). Simulation parameters: compartment volume = (1μm)3 with periodic boundary, R=0.01μm, l=0.01×1.0209μm, DE=DS=1μm2s1, ka=10kD, [S]=5[E]=0.001Nv, duration = 0.02 s.

    Reuse & Permissions
×

Sign up to receive regular email alerts from Physical Review E

Reuse & Permissions

It is not necessary to obtain permission to reuse this article or its components as it is available under the terms of the Creative Commons Attribution 4.0 International license. This license permits unrestricted use, distribution, and reproduction in any medium, provided attribution to the author(s) and the published article's title, journal citation, and DOI are maintained. Please note that some figures may have been included with permission from other third parties. It is your responsibility to obtain the proper permission from the rights holder directly for these figures.

×

Log In

Cancel
×

Search


Article Lookup

Paste a citation or DOI

Enter a citation
×