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  and single fluorophore video imaging  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  <cit>  and mrnas in the nucleus  <cit> .
these techniques have also enabled the sizes of microdomain structures  <cit>  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 .
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  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  termed "rafts"  <cit> .
furthermore, there we obtained important observations implying  molecular trajectories of raft constructs give rise to characteristic types of diffusion associated with "hop diffusion"  <cit>   the production of membrane protein complexes are facilitated by entering a clustering domain  the escaping rate of protein complexes from the clustering domain seems to be less than that of unbound molecules  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 , where nx, ny, and nz 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 →  <dig>  this type of random walk can be considered a wiener process with a time-dependent gaussian distribution with the following diffusion coefficient  <cit> :
the fastest rate of diffusion here is equivalent to d = 1/ <dig>  meaning that the particle must take the random walk step.
the binding process
suppose that a particle of chemical species s has just entered the interaction range of another particle of t species through the movement trial described above, and that these particles may bind to each other.
in this binding process, whether the particles can form a st complex must be determined.
specifically, we need to obtain both the candidates for complexes in a predefined table and the complexes that have already formed.
figure  <dig> shows a typical example of the case of a binary complex.
for simplicity, we consider only the dynamics projected onto a particular plane.
here, s, t, and u denote the chemical species of particles.
the area enclosed by the dash-dot circle indicates the interaction range for particle s, which is defined as a sphere of radius √ <dig>  in this case, the binding process consists of the series of steps listed below.
<dig>  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 .
<dig>  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 ), particle s is able to bind with particle t .
<dig>  a uniform random number ξ  is generated and compared with probability p = exp ≡ p <dig>  where Δe1/kbt is a dimensionless activation energy.
when ξ <p <dig>  the movement of particle s is accepted; when ξ ≥ p <dig>  it is rejected .
here we assume that the rate of the production of st complex is related to the conditional probability p <dig>
<dig>  when the movement is accepted, we assign t to the bond variable of particle s, and vice versa .
conservation of stoichiometry
it is not necessarily true that all pairs of s and t particles within the interaction range can bind to each other.
figure  <dig> shows this exceptional case, in which molecule t has already bound to a molecule of species u and therefore cannot bind to molecule s. in this case, these particles undergo the series of steps listed below.
<dig>  the movement trial drives particle s upward, thereby causing particle t to enter the interaction range .
<dig>  by referring to the predefined table, it is found that molecule s can bind with molecule t .
<dig>  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 .
the dissociation process
in this process, the species s and t in the bond variables are assigned to particles t and s, respectively, and are cleared upon the acceptance of unbinding.
<dig>  a uniform random number ξ  is generated and compared with probability p = exp ≡ p <dig>  where Δe2/kbt is a dimensionless activation energy.
when ξ <p <dig>  the movement of particle s is accepted; when ξ ≥ p <dig>  it is rejected .
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 Δe2/kbt  <cit> .
<dig>  when the movement is accepted, t in the bond variable of particle s is cleared, and vice versa .
the modification process
this process allows each particle to undergo a trial in which its chemical species is replaced by a different one.
the steps of the procedure are as follows:
<dig>  a uniform random number ξ  is generated and compared with probability p = exp ≡ p <dig>  where Δe3/kbt is a dimensionless activation energy.
when ξ <p <dig>  the modification of particle s to v is accepted; when ξ ≥ p <dig>  it is rejected .
<dig>  when the modification is accepted, s is replaced by v .
time and length scales
we define unit time τ of a simulation as a cycle in which every particle has undergone the movement trial only once.
also the trial of the modification process is performed only once for every particle per single unit time.
we then relate the unit time and unit length to real ones.
specifically, our choices of scale conversions are:  when we are interested in relatively fast dynamics within a small volume,  <dig> sec ≡  <dig> × 105τ and  <dig> μm ≡  <dig> l, or  when we study the long-time behavior of reactions within a relatively large volume,  <dig> sec ≡  <dig> × 103τ and  <dig> μm ≡  <dig> l .
these two combinations of scale conversions are selected so that both give d ≈  <dig>   when d = 1/ <dig>  i.e, the fastest diffusion.
this is readily checked by the following calculation.
since the diffusion coefficient is estimated using d ≡ 3l2d as an approximation of equation , and by noting that l =  <dig>  × 10- <dig>  and d =  <dig> × 105/ <dig>  for case , and l =  <dig>  × 10- <dig>  and d =  <dig> × 103/ <dig>  for case , we obtain d = 3l2d ≈  <dig>   equally for  and .
this value approximates experimental values for globular proteins in the cytoplasm  <cit> .
kinetic coefficients
RESULTS
dynamics of bound molecules
to show that our algorithm and its implementation – especially for the binding process – correctly function as designed in a simulation, figure  <dig> displays a trajectory of an se complex with p <dig> =  <dig> .
within this time range, molecules s and e remain continually bound to each other, while undergoing random movements in space .
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 .
in our monte carlo  simulations, we estimated average time evolutions over  <dig> samples, starting with a random initial distribution of particles.
figure  <dig> plots the average of time evolution of  for various initial concentrations  <dig> and two sets of activation energies listed in table  <dig> and  <dig> .
we confirmed that the ratio / of steady-state values correctly reproduces exp  with plotting symbol ▽  in figure 6a and 6b), where Δg = Δef1-Δer1-Δef2+Δer2+Δef3-Δer <dig> =  <dig> kbt.
to make comparisons with rate equations, we applied the scale conversion rule in table  <dig> to the activation energies in table  <dig> and  <dig>  thereby obtained the corresponding kinetic parameters in table  <dig> and  <dig>  respectively.
the curved lines in figure  <dig> are then drawn by numerically solving the following rate equations:
with the following constraints:
where, x, y, and z denote , , and , 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  phase transition into a liquid-ordered  one.
in the intermediate level of this addition, the lo phase and liquid-disordered  phase may coexist.
this may give rise to lateral heterogeneity in the plasma membrane, thereby segregating cholesterols into cholesterol-rich domains  <cit> .
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  <cit> .
taking these factors into account, we incorporated key components and their mutual interaction parameters other than those related to constituents in the bulk domain.
we considered the interactions among cholesterols , glycosphingolipids , tcr , and lat  molecules: a lower magnitude of coupling is assumed for cg, gt, and gl complexes, and a higher affinity is assumed for the tl complex .
with this setup and a sufficient amount of cg complexes, stable clustering patterns are obtained .
in figure 7c and 7g, components segregate and form clustering domains , reproducing drm/raft-like structures.
the improvement made in our simulation from those described in previous studies  <cit>  is demonstration of oligomerization-induced clustering, whereby the presence of a small number of plasma membrane protein complexes bound with relatively strong coupling correlates with the production of stabilized rafts  that are rich in c and g  <cit> .
in fact, the production of tl complexes is facilitated by the corralling of t and l in the clustering domains, while the slow motion of the tl complexes within the clusters in turn appears to hinder breakup of the clusters.
to make comparisons with spt as well as sfvi experiments, in particular, we displayed typical sample trajectories obtained from this simulation .
we observed two types of diffusional movements in the presence of clustering: the slow movement observed mainly for t and l , and the relatively free motion exhibited by c and g  but with a slower diffusion rate than that in the absence of clustering.
interestingly, these trajectories, especially that of g in figure 8d, appear to exhibit lumps with a similar size , indicating "hop diffusion" across clusters.
this size approximately equals the diameter of the clustering domain in figure  <dig>  for t and l, this hop diffusion rarely occurs, because most t's and l's rapidly form tl complexes in the clusters, and thus the energy cost of escaping from them increases.
hence, in figures 8a and 8b, t manages to escape from a cluster and enters another one wherein l remains to be trapped.
experimentally, images of trajectories of single or small groups of dope  molecules on the membrane of an nrk fibroblastic cell recorded at time resolutions of  <dig> μs were tracked, and found to exhibit hop diffusion with a typical scale of around  <dig>  μm  <cit> .
in the dynamics of membrane proteins, the rate of diffusion of lck  was reduced with decreasing distance from the point of stimulation in a tcr cluster  <cit> .
these observations are consistent with our simulation results.
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  <cit>  based on gillespie's algorithm  <cit> .
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  <cit> .
we are now performing systematic analyses of quantification of intrinsic noise by examining its dependence on kinetic constants such as km and vmax.
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"  <cit>  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  <cit> , with these proteins insulated from others with different signals.
using simulation to verify this will require  automatic identification of clustering domains or rafts, and  analyses of the net flow transferred between these domains.
these topics will be addressed in the near future.

 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 and single fluorophore video imaging 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 and mrnas in the nucleus. these techniques have also enabled the sizes of microdomain structures 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. 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 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 termed "rafts". furthermore, there we obtained important observations implying molecular trajectories of raft constructs give rise to characteristic types of diffusion associated with "hop diffusion" the production of membrane protein complexes are facilitated by entering a clustering domain the escaping rate of protein complexes from the clustering domain seems to be less than that of unbound molecules 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, where nx, ny, and nz 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 → this type of random walk can be considered a wiener process with a time-dependent gaussian distribution with the following diffusion coefficient: the fastest rate of diffusion here is equivalent to d = 1/ meaning that the particle must take the random walk step. the binding process suppose that a particle of chemical species s has just entered the interaction range of another particle of t species through the movement trial described above, and that these particles may bind to each other. in this binding process, whether the particles can form a st complex must be determined. specifically, we need to obtain both the candidates for complexes in a predefined table and the complexes that have already formed. figure shows a typical example of the case of a binary complex. for simplicity, we consider only the dynamics projected onto a particular plane. here, s, t, and u denote the chemical species of particles. the area enclosed by the dash-dot circle indicates the interaction range for particle s, which is defined as a sphere of radius √ in this case, the binding process consists of the series of steps listed below. 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. 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 ), particle s is able to bind with particle t. a uniform random number ξ is generated and compared with probability p = exp ≡ p where Δe1/kbt is a dimensionless activation energy. when ξ <p the movement of particle s is accepted; when ξ ≥ p it is rejected. here we assume that the rate of the production of st complex is related to the conditional probability p when the movement is accepted, we assign t to the bond variable of particle s, and vice versa. conservation of stoichiometry it is not necessarily true that all pairs of s and t particles within the interaction range can bind to each other. figure shows this exceptional case, in which molecule t has already bound to a molecule of species u and therefore cannot bind to molecule s. in this case, these particles undergo the series of steps listed below. the movement trial drives particle s upward, thereby causing particle t to enter the interaction range. by referring to the predefined table, it is found that molecule s can bind with molecule t. 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. the dissociation process in this process, the species s and t in the bond variables are assigned to particles t and s, respectively, and are cleared upon the acceptance of unbinding. a uniform random number ξ is generated and compared with probability p = exp ≡ p where Δe2/kbt is a dimensionless activation energy. when ξ <p the movement of particle s is accepted; when ξ ≥ p it is rejected. 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 Δe2/kbt. when the movement is accepted, t in the bond variable of particle s is cleared, and vice versa. the modification process this process allows each particle to undergo a trial in which its chemical species is replaced by a different one. the steps of the procedure are as follows: a uniform random number ξ is generated and compared with probability p = exp ≡ p where Δe3/kbt is a dimensionless activation energy. when ξ <p the modification of particle s to v is accepted; when ξ ≥ p it is rejected. when the modification is accepted, s is replaced by v. time and length scales we define unit time τ of a simulation as a cycle in which every particle has undergone the movement trial only once. also the trial of the modification process is performed only once for every particle per single unit time. we then relate the unit time and unit length to real ones. specifically, our choices of scale conversions are: when we are interested in relatively fast dynamics within a small volume, sec ≡ × 105τ and μm ≡ l, or when we study the long-time behavior of reactions within a relatively large volume, sec ≡ × 103τ and μm ≡ l. these two combinations of scale conversions are selected so that both give d ≈ when d = 1/ i.e, the fastest diffusion. this is readily checked by the following calculation. since the diffusion coefficient is estimated using d ≡ 3l2d as an approximation of equation, and by noting that l = × 10- and d = × 105/ for case, and l = × 10- and d = × 103/ for case, we obtain d = 3l2d ≈ equally for and. this value approximates experimental values for globular proteins in the cytoplasm. kinetic coefficients  dynamics of bound molecules to show that our algorithm and its implementation – especially for the binding process – correctly function as designed in a simulation, figure displays a trajectory of an se complex with p =. within this time range, molecules s and e remain continually bound to each other, while undergoing random movements in space. 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. in our monte carlo simulations, we estimated average time evolutions over samples, starting with a random initial distribution of particles. figure plots the average of time evolution of for various initial concentrations and two sets of activation energies listed in table and. we confirmed that the ratio / of steady-state values correctly reproduces exp with plotting symbol ▽ in figure 6a and 6b), where Δg = Δef1-Δer1-Δef2+Δer2+Δef3-Δer = kbt. to make comparisons with rate equations, we applied the scale conversion rule in table to the activation energies in table and thereby obtained the corresponding kinetic parameters in table and respectively. the curved lines in figure are then drawn by numerically solving the following rate equations: with the following constraints: where, x, y, and z denote,, and, 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 phase transition into a liquid-ordered one. in the intermediate level of this addition, the lo phase and liquid-disordered phase may coexist. this may give rise to lateral heterogeneity in the plasma membrane, thereby segregating cholesterols into cholesterol-rich domains. 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. taking these factors into account, we incorporated key components and their mutual interaction parameters other than those related to constituents in the bulk domain. we considered the interactions among cholesterols, glycosphingolipids, tcr, and lat molecules: a lower magnitude of coupling is assumed for cg, gt, and gl complexes, and a higher affinity is assumed for the tl complex. with this setup and a sufficient amount of cg complexes, stable clustering patterns are obtained. in figure 7c and 7g, components segregate and form clustering domains, reproducing drm/raft-like structures. the improvement made in our simulation from those described in previous studies is demonstration of oligomerization-induced clustering, whereby the presence of a small number of plasma membrane protein complexes bound with relatively strong coupling correlates with the production of stabilized rafts that are rich in c and g. in fact, the production of tl complexes is facilitated by the corralling of t and l in the clustering domains, while the slow motion of the tl complexes within the clusters in turn appears to hinder breakup of the clusters. to make comparisons with spt as well as sfvi experiments, in particular, we displayed typical sample trajectories obtained from this simulation. we observed two types of diffusional movements in the presence of clustering: the slow movement observed mainly for t and l, and the relatively free motion exhibited by c and g but with a slower diffusion rate than that in the absence of clustering. interestingly, these trajectories, especially that of g in figure 8d, appear to exhibit lumps with a similar size, indicating "hop diffusion" across clusters. this size approximately equals the diameter of the clustering domain in figure for t and l, this hop diffusion rarely occurs, because most t's and l's rapidly form tl complexes in the clusters, and thus the energy cost of escaping from them increases. hence, in figures 8a and 8b, t manages to escape from a cluster and enters another one wherein l remains to be trapped. experimentally, images of trajectories of single or small groups of dope molecules on the membrane of an nrk fibroblastic cell recorded at time resolutions of μs were tracked, and found to exhibit hop diffusion with a typical scale of around μm. in the dynamics of membrane proteins, the rate of diffusion of lck was reduced with decreasing distance from the point of stimulation in a tcr cluster. these observations are consistent with our simulation results.  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 based on gillespie's algorithm. 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. we are now performing systematic analyses of quantification of intrinsic noise by examining its dependence on kinetic constants such as km and vmax. 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" 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, with these proteins insulated from others with different signals. using simulation to verify this will require automatic identification of clustering domains or rafts, and analyses of the net flow transferred between these domains. these topics will be addressed in the near future. 
these domains corralled the diffusive motion of membrane proteins. 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. we propose an algorithm designed for three-dimensional simulation of the reaction-diffusion dynamics of molecules, based on a particle model. even single-particle tracking and single fluorophore video imaging are now being applied to observation of molecular-level dynamics. sample trajectories of raft constructs exhibited "hop diffusion". the simulation algorithm was tested for a reversible enzyme reaction model and its validity was confirmed.
 spatio-temporal dynamics within cells can now be visualized at appropriate resolution, due to the advances in molecular imaging technologies. even single-particle tracking and single fluorophore video imaging 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.  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 associated with rafts. sample trajectories of raft constructs exhibited "hop diffusion". these domains corralled the diffusive motion of membrane proteins.  these findings demonstrate that our approach is promising for modelling the localization properties of biological phenomena. june 20– symposium of computations in bioinformatics and bioscience in conjunction with the international multi-symposiums on computer and computational sciences hangzhou, china 