Volume 7 Supplement 4

## Symposium of Computations in Bioinformatics and Bioscience (SCBB06)

# Particle simulation approach for subcellular dynamics and interactions of biological molecules

- Ryuzo Azuma
^{1}Email author, - Tetsuji Kitagawa
^{2}, - Hiroshi Kobayashi
^{3}and - Akihiko Konagaya
^{1, 2}

**7(Suppl 4)**:S20

**DOI: **10.1186/1471-2105-7-S4-S20

© Azuma et al; licensee BioMed Central Ltd 2006

**Published: **12 December 2006

## Abstract

### Background

Spatio-temporal dynamics within cells can now be visualized at appropriate resolution, due to the advances in molecular imaging technologies. Even single-particle tracking (SPT) and single fluorophore video imaging (SFVI) are now being applied to observation of molecular-level dynamics. However, little is known concerning how molecular-level dynamics affect properties at the cellular level.

### Results

We propose an algorithm designed for three-dimensional simulation of the reaction-diffusion dynamics of molecules, based on a particle model. Chemical reactions proceed through the interactions of particles in space, with activation energies determining the rates of these chemical reactions at each interaction. This energy-based model can include the cellular membrane, membranes of other organelles, and cytoskeleton. The simulation algorithm was tested for a reversible enzyme reaction model and its validity was confirmed. Snapshot images taken from simulated molecular interactions on the cell-surface revealed clustering domains (size ~0.2 μm) associated with rafts. Sample trajectories of raft constructs exhibited "hop diffusion". These domains corralled the diffusive motion of membrane proteins.

### Conclusion

These findings demonstrate that our approach is promising for modelling the localization properties of biological phenomena.

## Background

We propose here a general method of simulating the dynamics and interactions of molecules based on a particle framework. Analyses of subcellular localization are now quite important to know how the cellular properties of interest are regulated. Experimental techniques helped unraveling these properties. For example, single particle tracking (SPT) and single fluorophore video imaging (SFVI) techniques have enabled observation of how individual molecules actually move and interact in space and time. SPT and SFVI have been used to investigate the dynamics of receptors in the plasma membrane [1, 2] and mRNAs in the nucleus [3–5]. These techniques have also enabled the sizes of microdomain structures [6, 7] to be measured. While some of these observations provided quantitative data, as in the case of SPT/SFVI studies, most merely provided a multitude of qualitative information focused primarily on biological significance. Moreover, the length and time scales in these experiments were approximately intermediate in scale between those that may be addressed by typical microscopic and macroscopic simulations (i.e., molecular dynamics and rate equations). Hence, our objective was to provide a simulation tool useful for integrating and examining experimental data quantitatively at the "mesoscopic" scale.

Our simulation method incorporates Brownian motion of molecules in 3D space. Interactions of molecules in the space are due to the production of complexes. The association and dissociation of these complexes, coupled with the modification of substrates to products, are accepted at certain rates based on a Monte Carlo (MC) algorithm which takes account of changes in their energy states. Although our setup is based on this type of molecular-level description, its ensemble average has been confirmed to correctly reproduce predictions derived from thermodynamic and rate equation theories. By using this method we succeeded to demonstrate the production of clustering patterns on the cellular membrane associated with cholesterol-rich detergent-resistant membranes (DRM) termed "rafts" [8]. Furthermore, there we obtained important observations implying (1) molecular trajectories of raft constructs give rise to characteristic types of diffusion associated with "hop diffusion" [9] (2) the production of membrane protein complexes are facilitated by entering a clustering domain (3) the escaping rate of protein complexes from the clustering domain seems to be less than that of unbound molecules (4) hence the membrane proteins are corralled in these domains most of the time, and in turn appear to stabilize the clustering domains. These results suggest the usefulness of our simulation approach.

This paper is organized as follows. In **Methods** we explain the random walk, binding, dissociation, and catalysis processes in our MC simulation algorithm. Then we discuss the correspondence of our probability constants to kinetic parameters of the corresponding rate equations. In **Results** we compare our simulation result of a reversible enzyme reaction model with the prediction from the rate equation theory. Then we show the result of cell-surface clustering simulation. The final section is devoted to our conclusions.

## Methods

### The random walk process

In this process, each particle moves along a cubic lattice, taking random steps to reach one of the six nearest neighbor sites with equal probability. The steps are of length *l*, so that each subsequent particle position can only take the values (*n*_{
x
}*l*, *n*_{
y
}*l*, *n*_{
z
}*l*), where *n*_{
x
}, *n*_{
y
}, and *n*_{
z
}are integers. This process allows the particle to take steps with probability *d* per unit time (τ), meaning that the particle waits at each point for a variable amount of time. A master equation theory can show that at the limit *l* → 0, this type of random walk can be considered a Wiener process with a time-dependent Gaussian distribution with the following diffusion coefficient [10]:

The fastest rate of diffusion here is equivalent to *d* = 1/6[τ^{-1}], meaning that the particle must take the random walk step.

### The binding process

- 1.
Particle S moves upward, as indicated by the pointer, causing another particle T to enter the interaction range. Here, the symbols {φ} attached to these particles denote an empty variable and indicate the absence of bound particles. The variable is hereinafter referred to as the bond variable (Figure 1A).

- 2.
The ability of S to bind with T is checked by referring to the predefined table for S. More precisely, when the table has a combination of single S and single T (represented by S-T in (b)), particle S is able to bind with particle T (Figure 1B).

- 3.
A uniform random number ξ (0 ≤ ξ < 1) is generated and compared with probability

*P*(ST|S+T) = exp(-Δ*E*_{1}/*k*_{B}*T*) ≡*p*_{1}, where Δ*E*_{1}/*k*_{B}*T*is a dimensionless activation energy. When ξ <*p*_{1}, the movement of particle S is accepted; when ξ ≥*p*_{1}, it is rejected (Figure 1C). Here we assume that the rate of the production of ST complex is related to the conditional probability*p*_{1}. - 4.
When the movement is accepted, we assign T to the bond variable of particle S, and vice versa (Figure 1D).

### Conservation of stoichiometry

- 1.
The movement trial drives particle S upward, thereby causing particle T to enter the interaction range (Figure 2A).

- 2.
By referring to the predefined table, it is found that molecule S can bind with molecule T (Figure 2B).

- 3.
The bond variable in molecule T is checked, yielding species U. This prevents the assignment of T to the bond variable of particle S, since S cannot be assigned to the variable of T (Figure 2C).

### The dissociation process

- 1.
A uniform random number ξ (0 ≤ ξ < 1) is generated and compared with probability

*P*(S+T|ST) = exp(-Δ*E*_{2}/*k*_{B}*T*) ≡*p*_{2}, where Δ*E*_{2}/*k*_{B}*T*is a dimensionless activation energy. When ξ <*p*_{2}, the movement of particle S is accepted; when ξ ≥*p*_{2}, it is rejected (Figure 3A). The basis of this mechanism is Kramers' theory, which predicts that the reaction rate of S+T→ST can be written as an exponential function of Δ*E*_{2}/*k*_{B}*T*[11].

- 2.
When the movement is accepted, T in the bond variable of particle S is cleared, and vice versa (Figure 3B).

### The modification process

- 1.
A uniform random number ξ (0 ≤ ξ < 1) is generated and compared with probability

*P*(VT|ST) = exp(-Δ*E*_{3}/*k*_{B}*T*) ≡*p*_{3}, where Δ*E*_{3}/*k*_{B}*T*is a dimensionless activation energy. When ξ <*p*_{3}, the modification of particle S to V is accepted; when ξ ≥*p*_{3}, it is rejected (Figure 4A).

- 2.
When the modification is accepted, S is replaced by V (Figure 4B).

### Time and length scales

^{5}τ and 1 μm ≡ 181.9

*l*, or (b) when we study the long-time behavior of reactions within a relatively large volume, 1 sec ≡ 5 × 10

^{3}τ and 1 μm ≡ 18.19

*l*(Table 1).

Relationships between probability constants and kinetic coefficients.

(a) | (b) | |
---|---|---|

1 sec | 5 × 10 | 5 × 10 |

1 μm | 181.9 | 18.19 |

| 45 | 45 |

| 38 | 38 |

| 5 × 10 | 5 × 10 |

| 5 × 10 | 5 × 10 |

These two combinations of scale conversions are selected so that both give *D* ≈ 7.6 [μm^{2}/sec] when *d* = 1/6[τ^{-1}], i.e, the fastest diffusion. This is readily checked by the following calculation. Since the diffusion coefficient is estimated using *D* ≡ 3*l*^{2}*d*[τ^{-1}] as an approximation of Equation (1), and by noting that *l* = 5.498 × 10^{-3} [μm] and *d* = 5 × 10^{5}/6 [sec^{-1}] for case (a), and *l* = 5.498 × 10^{-2} [μm] and *d* = 5 × 10^{3}/6 [sec^{-1}] for case (b), we obtain *D* = 3*l*^{2}*d*[τ^{-1}] ≈ 7.6 [μm^{2}/sec] equally for (a) and (b). This value approximates experimental values for globular proteins in the cytoplasm [12].

### Kinetic coefficients

Table 1 summarizes the correspondence between the probability constants and kinetic coefficients theoretically derived for cases (a) and (b). Here, *k*_{1}, *k*_{-1}, and *k*_{2} denote the kinetic coefficients for reactions S+T→ST, ST→S+T, and ST→VT, respectively. The second-order kinetic coefficient (*k*_{1}) is approximated by 4π*R* × 2*D* × *p*_{1}, as predicted from diffusion-limited reaction rate theory in three-dimensional solutions in the limit case of *t* → ∞ [13, 14]. Here, *R* is the reaction radius and approximated by *l*. The factor 1/3 in *k*_{-1} indicates that the total number of movements with unbinding for all possible binding conformations of S for fixed T is one-third of the total number of movements, with or without unbinding for all possible binding conformations of S for fixed T.

## Results

### Dynamics of bound molecules

*p*

_{2}= 0 (i.e., permanent binding). Within this time range, molecules S and E remain continually bound to each other, while undergoing random movements in space (a movie presentation of this can be obtained in [15]).

### Reversible enzyme reaction

To determine in more quantitative fashion whether our technique satisfies the theoretical requirements for the binding, unbinding, and modification processes, we conducted a series of simulations for the following reversible enzyme reactions:

Then we compared these reactions with the results obtained from rate equations. Here, the variables in each reaction denote the combination of a kinetic coefficient and dimensionless activation energy (e.g., *a*_{1} and ΔEf1 in S+E→SE reaction).

_{0}and two sets of activation energies listed in Table 2(a) and 2(b) (open symbols in Figure 6A and 6B, respectively). We confirmed that the ratio [P]/[S] of steady-state values correctly reproduces exp(-ΔG/

*k*

_{B}

*T*) (e.g., coincidence between steady state values for [S]

_{0}= 10 (or 4) with plotting symbol ▽ (or ◊) in Figure 6A and 6B), where ΔG = ΔEf1-ΔEr1-ΔEf2+ΔEr2+ΔEf3-ΔEr3 = 0.7

*k*

_{B}

*T*.

Parameter values of the reversible enzyme reaction model in Eq. 2. The translation of activation energies into kinetic coefficients are performed based on Table 1(a). *K*_{eq} = *a*_{1}*a*_{2}*b*_{-1}/*a*_{-1}*b*_{2}*b*_{1}

Activation energy [ | (a) (Figure 6A) | (b) (Figure 6B) | Kinetic coefficients | (c) (Figure 6A) | (d) (Figure 6B) |
---|---|---|---|---|---|

ΔEf1 | 0.1 | 1.7 |
| 5.73 | 1.16 |

ΔEr1 | 0.1 | 1.7 |
| 5.73 | 1.16 |

ΔEf2 | 2.3 | 3.9 |
| 2.79 × 10 | 5.62 × 10 |

ΔEr2 | 4.6 | 4.6 |
| 2.79 × 10 | 2.79 × 10 |

ΔEf3 | 10.8 | 10.8 |
| 1.02 × 10 | 1.02 × 10 |

ΔEr3 | 12.4 | 10.8 |
| 2.06 | 1.02 × 10 |

ΔG | 0.7 | 0.7 |
| 0.50 | 0.50 |

To make comparisons with rate equations, we applied the scale conversion rule in Table 1(a) to the activation energies in Table 2(a) and 2(b), thereby obtained the corresponding kinetic parameters in Table 2(c) and 2(d), respectively. The curved lines in Figure 6 are then drawn by numerically solving the following rate equations:

with the following constraints:

where, *x*, *y*, and *z* denote [SE], [PE], and [P], respectively.

### Extension to higher-order chemical reactions

The simulations thus far described have included only chemical reactions with complexes consisting of at most two different species. For general application of our simulation method, however, we must incorporate algorithms that address higher-order chemical reaction models that involve many-body interactions with more than two species. To do this, we developed an extended version of the algorithm by incorporating "bond tables" in which the indexes and species of bound partners for each particle are continually updated. We will not describe further details of the algorithm here because of their complexity, and will note only that they conform in essence to what we have described in **Methods**.

### Clustering on the cellular membrane

The advantages of this type of strict stoichiometric processing of interactions are enormous. One of the best demonstrations of the usefulness of the extended algorithm is simulation of clustering on the cellular membrane. This simulation examines the effects of associating LAT, a transmembrane adaptor protein, with T-cell receptor TCR on the lifetime of raft clustering. The T-cell receptor and LAT are shown to preferentially associate with cholesterol-rich microdomains on the cell surface.

Cholesterol has very important effects on the clustering, since it can change the properties of plasma membranes that mainly consist of sphingomyelins. The addition of cholesterol is shown to drive a solid-ordered (SO) phase transition into a liquid-ordered (LO) one. In the intermediate level of this addition, the LO phase and liquid-disordered (LD) phase may coexist. This may give rise to lateral heterogeneity in the plasma membrane, thereby segregating cholesterols into cholesterol-rich domains [8]. The effects of other factors relevant to clustering in the outer leaflet of plasma membranes are also emphasized. It has been proposed that the length and saturation of alkyl chains are responsible for clustering: glycosphingolipids, sphingomyelins, and phospholipids with long saturated alkyl chains may constitute rafts by entering the cholesterol-rich domains [16–19].

## Conclusions

Our simulation method represents molecular movements and interactions with the use of a particle model. The movements are expressed as a random walk in 3D space, and interactions are expressed as inter-particle binding, unbinding, and modification processes. A characteristic feature of our approach is that it enables complexes to be modeled as bound particles. The binding, unbinding, and modification processes of these complexes are based on energetic considerations, in which these processes proceed according to transition probabilities determined by the activation energy. This does not require setup of sub-volumes assuming well-stirred surrounding media, as incorporated in previous reaction-diffusion methods [25, 26] based on Gillespie's algorithm [27, 28]. Instead, our method requires precise treatment of the stoichiometry of these processes.

Moreover, it is important to examine the magnitude and effects of stochastic fluctuations, since our particle simulation model is essentially probabilistic. Our quantification of simulation results was intended to demonstrate that it correctly reproduced theoretical curves from a rate equation theory for the reversible enzyme reaction model. However, by observing raw time courses obtained for individual samples, continual fluctuations can be observed in raw time courses around the average value. This type of fluctuation corresponds to intrinsic noise, whereas fluctuation due to external input is referred to as extrinsic noise [29, 30]. We are now performing systematic analyses of quantification of intrinsic noise by examining its dependence on kinetic constants such as *K*_{m} and *V*_{max}.

Using this simulation method, we successfully demonstrated the existence of a clustering pattern that may be associated with receptor-cluster rafts and the "fluid mosaic model" [31, 32] in the plasma membrane. In immune cell signalling, rafts are hypothesized to be "platforms" for raft-philic adaptor proteins such as LAT in T-cells [33, 34], with these proteins insulated from others with different signals. Using simulation to verify this will require (1) automatic identification of clustering domains or rafts, and (2) analyses of the net flow transferred between these domains. These topics will be addressed in the near future.

## Declarations

### Acknowledgements

The authors wish to acknowledge useful discussions with T. Yamamoto. The computations in this work were performed using the Computer/Informatics Facilities, RIKEN Genomic Sciences Center, RIKEN super combined cluster (RSCC), and TSUBAME at Tokyo Institute of Technology Global Scientific Information and Computing Center.

This article has been published as part of *BMC Bioinformatics* Volume 7, Supplement 4, 2006: Symposium of Computations in Bioinformatics and Bioscience (SCBB06). The full contents of the supplement are available online at http://www.biomedcentral.com/1471-2105/7?issue=S4.

## Authors’ Affiliations

## References

- Daumas F, Destainville N, Millot C, Lopez A, Dean D, Salome L:
**Confined diffusion without fences of a g-protein-coupled receptor as revealed by single particle tracking.***Biophys J*2003,**84**(1):356–66.PubMed CentralView ArticlePubMedGoogle Scholar - Ritchie K, Kusumi A:
**Single-particle tracking image microscopy.***Methods Enzymol*2003,**360:**618–34.View ArticlePubMedGoogle Scholar - Shav-Tal Y, Darzacq X, Shenoy SM, Fusco D, Janicki SM, Spector DL, Singer RH:
**Dynamics of single mRNPs in nuclei of living cells.***Science*2004,**304**(5678):1797–800. 10.1126/science.1099754View ArticlePubMedGoogle Scholar - Shav-Tal Y, Singer RH, Darzacq X:
**Imaging gene expression in single living cells.***Nat Rev Mol Cell Biol*2004,**5**(10):855–61. 10.1038/nrm1494View ArticlePubMedGoogle Scholar - Fusco D, Accornero N, Lavoie B, Shenoy SM, Blanchard JM, Singer RH, Bertrand E:
**Single mRNA molecules demonstrate probabilistic movement in living mammalian cells.***Curr Biol*2003,**13**(2):161–7. 10.1016/S0960-9822(02)01436-7View ArticlePubMedGoogle Scholar - Murase K, Fujuware T, Umemura Y, Suzuki K, Iino R, Yamashita H, Saito M, Murakoshi H, Ritchie K, Kusumi A:
**Ultrafine membrane compartments for molecular diffusion as revealed by single molecule techniques.***Biophys J*2004,**86:**4075–4093. 10.1529/biophysj.103.035717PubMed CentralView ArticlePubMedGoogle Scholar - Kusumi A, Ike H, Nakada C, Murase K, Fujiwara T:
**Single-molecule tracking of membrane molecules: plasma membrane compartmentalization and dynamic assembly of raft-philic signaling molecules.***Semin Immunol*2005,**17:**3–21. 10.1016/j.smim.2004.09.004View ArticlePubMedGoogle Scholar - Barenholz Y:
**Sphingomyelin and cholesterol: from membrane biophysics and rafts to potential medical applications.***Subcell Biochem*2004,**37:**167–215.View ArticlePubMedGoogle Scholar - Fujiwara T, Ritchie K, Murakoshi H, Jacobson K, Kusumi A:
**Phospholipids undergo hop diffusion in compartmentalized cell membrane.***J Cell Biol*2002,**157:**1071–1081. 10.1083/jcb.200202050PubMed CentralView ArticlePubMedGoogle Scholar - Gardiner CW:
*Handbook of stochastic methods*. Springer, Berlin; 2004.View ArticleGoogle Scholar - Kramers HA:
**Brownian motion in a field of force and the diffusion model of chemical reactions.***Physica*1940,**7:**284–304. 10.1016/S0031-8914(40)90098-2View ArticleGoogle Scholar - Arrio-Dupont M, Foucault G, Vacher M, Devaux PF, Cribier S:
**Translational diffusion of globular proteins in the cytoplasm of cultured muscle cells.***Biophys J*2000,**78:**901–907.PubMed CentralView ArticlePubMedGoogle Scholar - Smoluchouski MV:
**Versuch einer mathematischen theorie der koagulationskinetik kolloider losungen.***Z Physic Chem*1917,**92:**129–168.Google Scholar - Collins FC, Kimbal GE:
**Diffusion-controlled reactions in liquid solutions.***Industrial and Engineering Chemistry*1949,**41:**2551–2553. 10.1021/ie50479a040View ArticleGoogle Scholar - [http://big.gsc.riken.jp/index_html/Members/azuma/folder.2006–06–13.0548915430/002.mp4]
- Subczynski WK, Antholine WE, Hyde JS, Kusumi A:
**Microimmiscibility and three-dimensional dynamic structures of phosphatidylcholine-cholesterol membranes: translational diffusion of a copper complex in the membrane.***Biochemistry*1990,**29:**7936–7945. 10.1021/bi00486a023View ArticlePubMedGoogle Scholar - Pasenkiewicz-Gierula M, Subczynski WK, Kusumi A:
**Influence of phospholipid unsaturation on the cholesterol distribution in membranes.***Biochemistry*1991,**73:**1311–1316. 10.1016/0300-9084(91)90094-HView ArticleGoogle Scholar - Gil T, Ipsen JH, Mouritsen OG, Sabra MC, Sperotto MM, Zuckermann MJ:
**Theoretical analysis of protein organization in lipid membranes.***Biochem Biophys Acta*1998,**1376:**245–266.PubMedGoogle Scholar - Bretscher MS, Munro S:
**Cholesterol and the Golgi apparatus.***Science*1993,**261:**1280–1281. 10.1126/science.8362242View ArticlePubMedGoogle Scholar - [http://big.gsc.riken.jp/index_html/Members/azuma/folder.2006–06–13.0548915430/AdditionalFile2.pdf]
- [http://big.gsc.riken.jp/index_html/Members/azuma/folder.2006–06–13.0548915430/112.mp4]
- Kusumi A, Koyama-Honda I, Suzuki K:
**Molecular dynamics and interactions for creation of stimulation-induced stabilized rafts from small unstable rafts.***Traffic*2004,**5:**213–230. 10.1111/j.1600-0854.2004.0178.xView ArticlePubMedGoogle Scholar - [http://big.gsc.riken.jp/index_html/Members/azuma/folder.2006–06–13.0548915430/112trj.mp4]
- Ike H, Kosugi A, Kato A, Iino R, Hirano H, Fujiwara T, Ritchie K, Kusumi A:
**Mechanism of Lck recruitment to the T-cell receptor cluster as studied by single-molecule-fluorescence video imaging.***Chemphyschem*2003,**4**(6):620–6. 10.1002/cphc.200300670View ArticlePubMedGoogle Scholar - Stundzia AB, Lumsden CJ:
**Stochastic simulation of coupled reaction-diffusion processes.***J Compt Phys*1996,**127:**196–207. 10.1006/jcph.1996.0168View ArticleGoogle Scholar - Elf J, Doncic A, Ehrenberg M:
**Mesoscopic reaction-diffusion in intracellular signaling.**In*Fluctuations and noise in biological, biophysical and biomedical systems*. Edited by: Bezrukov S et al. SPIE, Bellingham, WA; 2003:114–124.View ArticleGoogle Scholar - Gillespie DT:
**A general method for numerically simulating the stochastic time evolution of coupled chemical reactions.***J Compt Phys*1976,**22:**403–434. 10.1016/0021-9991(76)90041-3View ArticleGoogle Scholar - Gillespie DT:
**Exact stochastic simulation of coupled chemical reactions.***J Phys Chem*1977,**81**(25):2340–2361. 10.1021/j100540a008View ArticleGoogle Scholar - Elowitz MB, Levine AJ, Siggia ED, Swain PS:
**Stochastic gene expression in a single cell.***Science*2002,**297:**1183–1186. 10.1126/science.1070919View ArticlePubMedGoogle Scholar - Swain PS, Elowitz MB, Siggia ED:
**Intrinsic and extrinsic contributions to stochasticity in gene expression.***Proc Natl Acad Sci USA*2002,**99**(20):12795–12800. 10.1073/pnas.162041399PubMed CentralView ArticlePubMedGoogle Scholar - Singer SJ, Nicolson GL:
**The fluid mosaic model of the structure of cell membranes.***Science*1972,**175:**720–730. 10.1126/science.175.4023.720View ArticlePubMedGoogle Scholar - Barenholz Y, Cevc G:
**Structure and properties of membranes.**In*Physical Chemistry of Biological Surfaces*. Edited by: Baszkin A, Norde W. Marcel Dekker, New York; 2000:171–241.Google Scholar - Veillette A:
**Specialized adaptors in immune cells.***Curr Opin Cell Biol*2004,**16:**146–155. 10.1016/j.ceb.2004.01.002View ArticlePubMedGoogle Scholar - Leo A, Schraven B:
**Networks in signal transduction: the role of adaptor proteins in platelet activation.***Platelets*2000,**11:**429–445. 10.1080/09537100020027815View ArticlePubMedGoogle Scholar

## Copyright

This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.