Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Comput Biol Chem. Author manuscript; available in PMC 2010 June 1.
Published in final edited form as:
PMCID: PMC2693469

Hybrid stochastic simulations of intracellular reaction-diffusion systems


With the observation that stochasticity is important in biological systems, chemical kinetics have begun to receive wider interest. While the use of Monte Carlo discrete event simulations most accurately capture the variability of molecular species, they become computationally costly for complex reaction-diffusion systems with large populations of molecules. On the other hand, continuous time models are computationally efficient but they fail to capture any variability in the molecular species. In this study a novel hybrid stochastic approach is introduced for simulating reaction-diffusion systems. We developed a dynamic partitioning strategy using fractional propensities. In that way processes with high frequency are simulated mostly with deterministic rate-based equations, and those with low frequency mostly with the exact stochastic algorithm of Gillespie. In this way we preserve the stochastic behavior of cellular pathways while being able to apply it to large populations of molecules. In this article we describe this hybrid algorithmic approach, and we demonstrate its accuracy and efficiency compared with the Gillespie algorithm for two different systems. First, a model of intracellular viral kinetics with two steady states and second, a compartmental model of the postsynaptic spine head for studying the dynamics of Ca+2 and NMDA receptors.

Keywords: Biochemical Networks, Stochastic Simulation, Hybrid Algorithm, Chemical Master Equation, Calcium Dynamics

1 Introduction

With the observation that stochasticity is important in biological systems, chemical kinetics are now receiving wider interest. For example, the dendritic spine of a Purkinje cell has a volume of 10−19 (m3) and so 1 μM of a chemical corresponds to 60 molecules and 100 nM corresponds to six molecules ([Morishita and Aihara, 2004]). Because of this small number of molecules, stochasticity can considerably influence intracellular signal transduction ([Rao et al., 2002]). In 1976 Gillespie introduced two versions of a stochastic simulation algorithm (SSA) for simulating spatially homogeneous chemical systems, the Direct Method and the First Reaction method ([Gillespie, 1976]). Kinetic Monte Carlo algorithms ([Frenkel and Smith, 1996]) have been used extensively to study stochastic cell signaling processes, however, these methods are computationally expensive. This means as models become more and more complex the SSA becomes computationally slow, and therefore a rate limiting step in the study of cell signaling. Many efforts have tried to optimize or accelerate the SSA ([Gibson and Bruck, 2000, Cao et al., 2004, McColluma et al., 2005]), and approximate accelerated stochastic algorithms have been developed ([Rathinam et al., 2003, Goutsias, 2005, Cao et al., 2005, Chatterjee et al., 2005]). Gillespie introduced the explicit tau leap method, which approximates continuous-time discrete-state Markov process of fast-reactions as Poisson distributed ([Gillespie, 2001]), and hybrid algorithms that partition the reactions into Fast and Slow subsets have attempted to address this issue. By partitioning the reactions into these two subsets and simulating the Slow reactions with the SSA and the Fast reactions with an accelerated ([Puchalka and Kierzek, 2004, Rao et al., 2002]) or deterministic rate-based method ([Haseltine and Rawlings, 2002, Salis and Kaznessis, 2005]), the complexity of the system is reduced significantly and consequently the simulation running time.

In this study we have developed a novel hybrid stochastic approach for simulating reaction-diffusion systems. We assign each reaction channel and diffusion Markov jump event based on their propensities to two subsets, the Fast and the Slow. We treat the Fast set in the continuous time regime with deterministic rate based equations, and the Slow set in the discrete time regime with the SSA. In contrast with other similar methods ([Rudiger et al., 2007, Kiehl et al., 2004]), we dynamically partition each reaction, as well as each diffusion Markov event, to these two subsets during the whole simulated time. The key feature of our algorithm is the scaling of the propensity of each reaction/diffusion jump for each subset by a fractioning parameter. In this way, each reaction channel and diffusion Markov event is simultaneously simulated stochatically as well as deterministically. The level, at which each reaction is simulated stochastically is defined by a fractioning parameter derived by using a smooth, S-shaped function. In that way we define a fractional propensity for each reaction by scaling its original propensity function value. By using such a representation we avoid rapid changes of the propensity functions and generate smoother transitions of the reaction channels between discrete and continuous time regimes. In addition, by using a smooth function for partitioning and calculating the fractional propensities for the reactions, we capture their variability in a dynamic fashion during the whole simulation time. The portion of the captured variability is defined again from the fractional propensities. Finally, with this hybrid algorithm we achieve accelerated simulations while producing accurate results compared with the SSA. Large scale biochemical networks ([Manninen et al., 2006, Arkin et al., 1998]) have been used in the past as test beds for evaluating hybrid algorithms ([Wylie et al., 2006, Takahashi et al., 2004]). On the other hand simple reaction systems that are still computationally demanding and stiff in their nature, have been used extensively in many studies as well ([Alfonsi et al., 2005, Cao and Petzold, 2006, Chiam et al., 2006, Chatterjee et al., 2005, Rathinam et al., 2003]). The advantage of the small systems is that they provide us with a comprehendible insight that we often cannot achieve with large scale systems. In the present paper we demonstrate the application, accuracy and performance of our algorithm in the modeling of two different systems important in computational biology : (i) intracellular viral kinetics and (ii) calcium and N-methyl-D-aspartic acid receptor (NMDAR) dynamics in the postsynaptic spine head.

Mathematical models of viral dynamics have proved useful in describing viral progression upon infection and viral decline after drug treatment. The various kinetic models used cover situations ranging from generic temporal three-variable cases, which include uninfected cells, infected cells and virions, to complex spatio-temporal schemes taking into account the specifics of certain diseases ([Wodarz and Nowak, 2002, Perelson, 2002, Tuckwell and Wan, 2004]). Experimental data describing the details of intracellular interactions are now extensive ([Knipe and Howley, 2001, Poranen et al., 2002, Smith and Helenius, 2001, Pornillos et al., 2002]), and modelling the interplay of ensembles of virions and cells has yielded quantitative results that have greatly improved our understanding of immunological phenomena ([Kesmir and De Boer, 1999]). In this study we have used a nonlinear model of a generic non-lytic virus developed by Rsrivastava (2002) to test our algorithm. In this model, the viral nucleic acids were classified as either genomic or template, and two distinct properties of this system made it attractive for our study. First, the different time scales of the reactions provide us the opportunity to apply our algorithm by simulating fast and slow reactions with deterministic and stochastic algorithms, respectively. We were interested in testing the accuracy of our results for this model, and one way for doing so is to compare the first and second moment of the concentration of the different species. The second attractive property of the system is the existence of two steady states for the number of the template species. For these systems Cao, (2006) have introduced the concept of distribution distance for the evaluation of errors in stochastic simulation methods for chemically reacting systems. One of those is the Kolmogorov distance, and this is the criterion we have used to test our results. In doing so we demonstrate how our method, by using fractional propensities for the chemical reactions, can accurately reproduce results generated by the SSA.

The second application of our algorithm is the modeling of calcium and NMDAR dynamics in the postsynaptic spine head. Modeling studies often use compartmental models in spherical, cylindrical or Cartesian coordinates to simulate spine geometry ([Holmes, 1990, Volfobsky et al., 1999, Gold and Bear, 1993, Naoki et al., 2005]). In compartmental models the spatial localization of reactions in a cell is represented by independent compartments that interact only through diffusion of a small subset of reactants across compartmental boundaries. Freeware software ([Ichikawa, 2004, Franks et al., 2001]) which use stochastic or deterministic algorithms have been used widely. In reaction diffusion models localization is represented by spatial variables in the partial differential equations (PDEs). Key assumptions of such models are violated in subcellular domains when the number of molecules is statistically small or when many of the interacting molecules are immobilized. The dendritic spine is a good example of such a model. As Ca+2 signaling in spines is initiated by NMDARs it flows through the receptor channel and encounters a meshwork of binding effector molecules that are immobilized or diffuse only very slowly. This situation can be simulated by a stochastic model in which binding and reaction rates are represented as probabilities that can be calculated from kinetic rates measured in solution. Inside the dendritic spine, Ca+2 ions trigger signal transduction pathways that regulate a variety of neurobiological processes including synaptic plasticity, one of the major mechanisms underlying learning and memory. This is achieved by Ca+2’s role as a regulator of many postsynaptic enzymes, which trigger modification of synaptic strength or structure. Most forms of long-term synaptic plasticity, such as long-term potentiation (LTP) and long-term depression (LTD) require postsynaptic Ca+2 signals ([Zucker, 1999]), and we have shown how stochasticity can affect plasticity curves ([Shouval and Kalantzis, 2005]). Recently, Rudiger (2007) have employed a hybrid algorithm to describe calcium dynamics in the cytosol ([Rudiger et al., 2007]). In their algorithm they model the channels stochastically and the concentration fields of calcium deterministically. In our model of the dendritic spine head we partition dynamically each reaction and jump event during the whole simulation time. We show that our method is not only faster, but also consistent with the SSA with respect to the mean and the variance of calcium dynamics and NMDAR transitions. For simplicity we have used a cylindrical compartmental model for the spine head and the neck. The end of the spine neck behaves as a sink for the diffused molecules, and synaptic activity was modeled with a step pulse function of release of Glutamate with amplitude 1 mM and duration 1 ms. The NMDARs were on the upper cross sectional surface of the spine head, and they were simulated with a five-state Markov model ([Franks et al., 2002]). The only source of calcium was the NMDARs, and the postsynaptic membrane voltage was assumed to be constant, which implies that the rate of influx of calcium ions was fixed. The calcium reacts with calcium pumps and a generic buffer which was simulated with a two-state Markov model. For the numerical solution of ordinary differential equations (ODEs) or PDEs in continuous time there are several numerical methods. A forward Euler for the deterministic was sufficient and also efficient compared to higher order methods. For the SSA we have used the Direct method (DM) version.

2 Exact stochastic algorithm

2.1 Chemical reaction master equation

The Gillespie algorithm is based on the assumption that our system is a well-stirred mixture of N≥1 molecular species {S1,…, SN} that chemically interact, inside some fixed volume Ω through M ≥1 reactions {R1,…, RN}. We specify the dynamical state of this system by X(t) [equivalent] (X1(t),…, XN (t)), where Xi(t) [equivalent] the number of Si molecules in the system at time t, i = 1, …, N. For each reaction Rj there is a propensity function aj(x) which describes the probability, given X(t) = x, that one Rj reaction will occur somewhere inside the Ω in the next time interval [t, t + dt]. The propensity function aj(x) is equal to the product of a rate constant cj and the number of possible combinations of reactant molecules involved in reaction Rj. For the following reactions the propensity function is equal to: X → *, Xa + XbXc and Xa + XaXb, we get aj = cjXa(t), aj = cjXa(t)Xb(t) and aj = cjXa(t)(Xa(t) − 1)/2 respectively.

Once the reaction Rj occurred, the number of molecules for each species is updated according to the state vector or else the stoichiometric matrix νj which is equal to the change in the number of Si molecules. produced by one Rj reaction (j = 1,…, M ; i = 1,…, N).

The time evolution equation for the probability P(x, t|x0, t0) that X(t) will equal x given that X(t0) = x0 for tt0 is given by the chemical master equation (CME):


The Gillespie algorithm allows the numerical solution of the CME. Given the actual time t, the probability that the next stochastic event occurs in the infinitesimal time interval [t + τ, t + τ + dt] and is a Ri event, is given by


where a0 = Σj aj is the sum of all propensities. The probability P (τ, j) can be realized by drawing two random numbers r1 and r2 from a uniform distribution in the interval [0, 1] and choosing τ and i such that




In this way the next event to occur is Ri and it will occur after time τ. This is the Direct method introduced by Gillespie.

2.2 Reaction-Diffusion master equation

Applications of the Gillespie algorithm have been introduced for simulating spatially inhomogeneous reaction-diffusion systems in mesoscopic volumes. That method is based on the discretization of the volume in small cells and modeling the diffusion of molecules between neighboring cells as a Markov jump process ([Bernstein, 2005, Baras and Mansour, 1996]). Bernstein demonstrated that the algorithm is capable of simulating both pure diffusion and the non linear autocatalytic reaction-diffusion system under certain conditions for the size of each cell. For crowded systems where the probability of one molecule to jump to a neighboring cell is high, the algorithm can become again slow. Recent studies have developed also hybrid mesoscopic algorithms in order to reduce the running time for reaction diffusion systems ([Chiam et al., 2006, Wylie et al., 2006]).

The basic lines of the master equation formulation of reaction-diffusion systems can be summarized as follows. We subdivide the reaction volume into spatial cells ΔVr and consider as variables the numbers of particles xi;r(t) of species i = 1,…, M at time t at the cell r. We assume the jump process of one particle from cell r to an adjacent one as a Markov process. The mean jump frequency of species i, i, is related to Fick’s diffusion coefficient of the species with the following equation:


where d represents the space dimension and l is the characteristic length of the cell:


For the above system the resulting probability distribution P(x(r), t) obeys the so called multivariate master equation:



The sum k in eq. ?? runs over the first nearest neighbors of the cell r. We must note that the size of each compartment cannot be arbitrarily large. The linear dimension of the each cell must be typically of the order of the reactive mean free path which is defined as the average distance travelled by a particle before it undergoes a reactive collision. The local probability that the next j reaction will occur out of M distinct types at the volume V of the cell is obtained by modifying eq. 2:




and the time that the next reaction take place is given by:


2.3 Langevin Reaction Equation with Diffusion jumps

Assuming a chemical reaction system where fast reactions occur many times in an in-finitesimal time interval dt yet none of the propensity functions change appreciably, then these can be approximated as a continuous Markov processes. These continuous Markov processes are governed by Fokker-Planck equations which are more computationally tractable ([Risken, 1989]) than the master equations that govern jump Markov processes. These Markov processes for a system of chemical reactions is described as Chemical Langevin Equation (CLE) and can be interpreted as an Ito SDE (2) stated as:


where N(0,1) refers to the standard normal distribution. CLE have been used in the past for describing neuronal signal transduction pathways ([Manninen et al., 2006]). Many numerical methods have been developed for solving SDE. The simplest one is the Euler-Maruyama method which is the following:


where ΔWj is a normal Gaussian random number with a mean zero and a variance of Δt. The above equation can be modified for Reaction-Diffusion systems:


In the above equation xi;r is the number of molecules of type i in the compartment r which participate in M reactions.

3 Hybrid algorithm

In principle, for simulating cellular processes the ideal is the full, microscopic discrete event approach, as anything else is only an approximation. In this section we formulate a hybrid approach for simulating reaction or reaction-diffusion systems, which aims to preserve the stochastic behavior of the system while remaining efficient for simulating stiff systems with different time scales. This is achieved by combining the discrete-time SSA with the time-step integration of PDEs. In order to do this we partition the set of reaction channels and the diffusion Markov events into a ‘Fast’ or ‘Slow’ subset based on its propensity for the two groups. The Fast subset is implemented in the continuous time regime with PDEs, and the Slow subset is implemented in the discrete time regime with the SSA. Our approach is novel in that we apply fractional propensity functions for each reaction or diffusion jump event, and the calculation of these fractional propensities is attained using a smooth S-shaped function f(ai, t). One example of an S-shaped partitioning function is the sigmoid function:


We see that the sigmoid function (eq. 6) has two parameters. The first one, μ, defines the half maximum and the second one, λ, defines the slope of the curve (Fig. 1). Similarly with other approaches ([Griffith et al., 2006, Chiam et al., 2006, Alfonsi et al., 2005, Puchalka and Kierzek, 2004]) those two parameters are user-defined in our algorithm and there is not a systematic way to predefine them, rather the partitioning should be intuitive. Maintaining at least two orders of magnitude difference between the values of the partitioned reaction probabilities has been recommended in the past ([Haseltine and Rawlings, 2002]). That criterion seems to be empirically valid for proper separation of the time scales in a stiff system. Usually, preliminary results using the Gillespie algorithm help us to indicate possible bottlenecks in our system and set the partition parameters accordingly. The role of those two parameters will be discussed in the next section with the application of our algorithm to the viral kinetics.

Fig. 1
Partition sigmoid function f(x). Two parameters define the shape of the curve. The half maximum in both cases is the same (parameter μ) but the slope differs (parameter λ).

With this representation each reaction contributes to both the deterministic as well as the stochastic regime, and the degree of this contribution depends on the product of the propensity of each reaction or Markov jump event with the partitioning parameter f(ai, t). Those reactions that occur with high rates are simulated mostly on the continuous time regime where as reactions with low frequency mostly on the discrete time regime. In the absence of a hard cutoff value for classifying reaction channels into the Fast or Slow subsets our partitioning strategy leaves some reactions assigned in both continuous and discrete time regimes. In this case the molecular species that participate in the reaction are represented by a non integer number. The integer portion reflects changes due to the SSA algorithm, and the floating point value reflects changes due to the PDEs. However if a stochastic event occurs we may need to update the molecular populations even if it is not integer. In that case we either round off the concentrations or keep non integer numbers. Another approach is the “random-rounding” ([Bhalla, 2004]). In that case a non integer population is rounded randomly up or down to the nearest integer. In our simulations sets we maintained non integer numbers. The reason is that rounding off the populations whenever we have a stochastic event we may have accumulations of the rounding errors. Additionally we wanted to preserve the law of mass conservation.

For reaction-diffusion systems the modeling domain Ω is discretized into a Cartesian mesh of total NX ×NY ×NZ cells, each of dimensions NX1×NY1×NZ1. As was mentioned earlier molecular diffusion is modelled as a set of random jump processes from one cell to adjacent ones. Therefore, in each cell the set of M reactions and Markov jump events is carried out independently. The hybrid evolution equation of the species X(t) in one cell can be described with the following equation:


where the first term in the sum is the deterministic part, the second term is the stochastic and dNj(t) obeys the CME. The local probability that the next j reaction or jump event will occur out of M distinct types at the volume V of the cell r is obtained by:


where the sum of the scaled propensities is given by:


and the time interval until the next reaction of jump event:


In the continuous regime we solve numerically the differential equations of our system. However we need to keep in mind that numerical methods should satisfy conditions and criteria for accuracy and converge to solutions. Our algorithmic approach is summarized as follows:

  1. initialize and set t = 0
  2. while (t < T) do
  3. Calculate Propensities a(t) in Ω
  4. Calculate partitioning parameter f(ai;r, t)
  5. Drop a random number r1 from a uniform distribution
  6. Calculate τ(as(t), r1) for the fractional Slow propensities
  7. if (Δt) ≤ τ do
  8. t0 = t
  9. while (t + dtτ + t0) do
  10. Solve Deterministic for fractional Fast reactions: [Rd, X(t)]
  11. Update propensities a(t)
  12. Partition propensities to Fast [Rd(t), af (t)] and Slow [Rs(t), as(t)]
  13. Calculate τ(as(t), r1)
  14. t = t + Δt
  15. end while
  16. Find next reaction for (Rs(t), X(t), r2)
  17. else
  18. Find next reaction for (Rs(t), X(t), r2)
  19. t = t + τ
  20. end if
  21. end while

However during the time interval of two stochastic events the population of the different molecular species change and consequently the propensities are not constant. In that case we need to consider non homogeneous Poisson processes and integrate numerically the CME. In our approach, as we notice, we update the stochastic time step τ during the whole duration of simulation time and in combination with the sigmoid partition function we preserve smooth transitions of reactions from the continuous to the discrete time regime and vice versa. However, if the assumption of stationary processes results in large errors then we need to modify the evolution equation accordingly:


and the amount of time that must pass until the next reaction take place is given by the following equation:


The following changes could be applied to our algorithm:

  • (7) set ξ = 0
  • (8) while ξ < −ln(r1) do
  • (9) Solve Deterministic for fractional Fast reactions: [Rd, X(t)]
  • (10) Update propensities a(t)
  • (11) Partition propensities to Fast [Rd(t), af (t)] and Slow [Rs(t), as(t)]
  • (12) ξ=ξ+j=1M(1f(αj,t))αj(X(t))dt
  • (13) t = t + Δt
    (14) end while

It is worth noting that that integration of the CME may prove to be computationally expensive. Haseltine and Rawlings (2002) has proposed a modification of their hybrid algorithm to avoid the integration of the CME by scaling the stochastic time step t by artificially introducing a probability of no reaction into the system. by using fractional propensities we have slow changes of the propensities, and in most cases we can assume stationary Poisson processes for the reactions. Preliminary results from test simulations usually indicate us if we need to assume non homogeneous Poisson processes or not. Comparison of the mean and the variance could be used for evaluating our results. However as we shall see alternative metrics may be considered depending of the system we study. In Appendix we formulate two metric measurements for evaluating the accuracy of our approach.

4 Application of the algorithm for modelling intracellular viral kinetics

In the past, models have been developed for the simulating intracellular growth of bacteriophage T7 and HIV-1 for a better understanding of different anti-viral strategies. In this section we apply our algorithm to a nonlinear model for a non lytic virus, which was suggested by Srivastava et al (2002). The main components of this model are the viral genomic (gen) and template (tem) nucleic acids and the viral structural protein (struct). In this model the genome can either form a template, which takes part catalytically in the synthesis of the genome, or it may be packaged in structural proteins (Fig. 2). The model due to its stochastic behavior it provides us the opportunity to use it as a test bed for determining the accuracy as well as the performance of our hybrid algorithm.

Fig. 2
Diagram of the viral model. Notice two catalytic reactions k3 and k5

That model can be described with the following ODEs in a matrix form:


It is worth bearing in mind that the viral model has two catalytic reactions, k3 and k5. The parameters of eq. 8 adopted from Srivastava et al (2002).

As we can see from Table 1 the rate of gen production is 1 molecule−1day−1 where as the synthesis parameter of structural protein is 1000-fold greater. Because of this difference between the kinetic rates we can apply our algorithm for simulating the dynamics of the three species. In Figure 3 we see the histogram of the tem population after 200 days over 5,000 simulation runs and initial conditions (gem,struct,tem)=(0,0,1). In other words, it shows the probability distribution of the Template population at the end of the simulation time. The x-axis represents the template population and the y-axis the probability. Additionally it demonstrates the existence of two steady states as have been shown by Srivastava et al (2002). The first one, with the sharp peak, is unstable (or else “repulsive”) therefore if the system visit that state stays there or else is moving away. The second one is stable (or else “attractive”). Because of the stability of the second steady state we expect the existence of a basin of attraction. Due to that, we notice a greater spread of the probability distribution. However, because of the existence of two steady states low-order moments are not of much relevance for estimating and comparing our algorithm. ([Cao and Petzold, 2006]).

Fig. 3
Benchmarking simulations of the viral model for 5000 realizations using the SSA. From the histogram of the tem population after 200 days we can see the existence of the two steady states of the system.
Table 1
Parameters of the viral model

For our purpose we shall use the cumulative distribution function (cdf) FX, because the comparison of the cdf is more appropriate for estimating the errors of our simulations. Using data from simulations we estimate the cdf, and then we use the Kolmogorov distance to compare the cdf estimated from our hybrid algorithm’s results with the one derived using the Gillespie algorithm. In our case we do not know the exact cdf, but instead use the empirical version derived from our simulation data. One can see from Fig. 3 that the error is small, so we can accept the empirical obtained from the SSA as a good approximation of the theoretical cdf. In general, a distance measure quantifies the extent to which two states behave in the same way. While these distance measures are usually given by a certain mathematical expression, they often possess a simple operational meaning, i.e. they are related to the problem of distinguishing two systems. The Kolmogorov distance for two random variables X and Y with cdf FX and FY is defined as:




The Kolmogorov difference of the two distributions takes values between 0 and 1. A value close to zero implies that distributions are more alike while a value closer to 1 indicates that the distributions are more disjointed.

For the estimation of the fractional propensities we used a S-shaped, built-in membership function in Matlab. The initial conditions were (gem,struct,tem)=(0,0,1), and we processed results over 5000 trials. For demonstration purpose we have used three sets of the parameters μ and λ for the S-shaped partitioning function f(aj(x)). In the first two sets the main difference is the value of the parameter λ, which defines the slope of the partitioning function. For this first set of simulations, where the partitioning function is steep (μ=1 λ=1.5), the mean template population is 14.6273 and the variance 76.2997, but the cdf (Fig. 4, dashed line) deviates significantly from that obtained using the SSA algorithm (Fig. 4, solid line). By using a smoother partitioning function (μ=1.5 λ=25) with almost the same half maximum μ, the mean of the template is 13.9832 and the variance 80.4190, and a good match to the cdf obtained by the SSA is achieved (Fig. 4 dotted line). For the third set of parameters (Fig. 4 dashed-dotted line) μ=100 and λ=140, our results are almost identical with those obtained from the SSA algorithm, which is what one would expect.

Fig. 4
Mean results from 5000 realizations of the viral model. Comparison of the empirical cdf for the SSA (solid line), and our hybrid algorithm for three sets of the parameters μ and λ (dashed line, μ=1 λ=1.5, dotted line, ...

From our results it is clear that two steady states exist in the model. This prohibits us from using the first and second moment as a criterion for the estimation of errors for our hybrid simulation method. By using a hard threshold or a smooth partitioning function and applying our method we obtain similar results for the mean and variance of the template nucleic acid. However, only in the second case we do accurately reproduce their cdfs. Our method reliably links the two time regimes by allowing smooth transitions of the reaction channels between the two time regimes. Reactions with moderate propensities participate both in the continuous and discrete time regime, thereby contributing to the overall variability of the system. The degree of that participation depends on the fractional propensities obtained after scaling them with the fractioning parameter f(aj(X(t)). Quantitative results from measurements of the Kolmogorov distance are listed in Table 2.

Table 2
Kolmogorov distance between SSA and the hybrid method

Finally we tested the performance of our algorithm. We measured the average time spent executing one simulation run. Since the time interval, in the discrete regime, is proportional to the product of the kinetic rate and the number of species participating in each reaction we tested the performance of our algorithm for different initial conditions. As we can see in Table 3 the time of the SSA increase, as we increase the initial number of species. Even if the viral model is not computational demanding, the simulation time has been reduced significantly for each case when we use our hybrid approach. The reason for that is the nice time scale separation that we obtain for the particular system. Notice that the time required for the deterministic approach does not depend on the initial condition but only on the total number of time steps.

Table 3
Average simulation time

5 Calcium dynamics at the dendritic spine

5.1 Model of the spine head

In this section we describe the model we used for studyinge Calcium dynamics in the dendritic spine. Just as in the viral example described above we incorporate differential equations for the continuous time regime and the SSA for the discrete time regime. For simplicity, we shall only refer to the deterministic rate-based equations of the model, even though our simulation results for both the SSA and hybrid algorithm implement stochastic modeling during all reaction and diffusion events described below. The model has 13 compartments (Fig. 5a): six for the spine head and seven for the neck, and the most distant boundry of the spine neck acts as a trap for calcium ions. For compartment i, in the continuous time regime, the dynamics of calcium can be described with mass-action kinetics:


where DCa is the diffusion coefficient of the calcium ions, and L is the length of the i compartment. For the SSA the mean jump frequency of the Ca+2, Ca is given by:

Fig. 5
(a) Postsynaptic spine:(Left) Calcium enters the spine head through the NMDA receptors. Ions diffuse inside the spine, react with the calcium buffer B or leave to the extracellular space through the pumps.(Right) Compartmental model of the dendritic spine. ...

For the coupling of diffusion from the head to the neck, we need to scale the diffusion coefficient with the ratio of their cross sectional areas.

The ATP-driven pumps were modelled with the following equation:


where kpi is the rate constant computed from the velocity of the ATP-driven pump (1.4 × 10−4 cm/msec) multiplied by Ai/Vi, where Ai is the membrane area of compartment i and Vi is its volume.

In the dendritic spine there is a large family of calcium binding proteins, among them calmodulin, calbindin and calcineurin. In our model we have used a 2-states Markov model for a generic calcium buffer. The reactions of calcium with the buffer can be described by the following two equations:


where [B] and [CaB] is the concentrations of free buffer and bound buffer with calcium respectively. The parameters k1buffer and k1buffer are the association and dissociation rates.

Calcium enters the spine head from the outermost compartment of the head through the NMDARs, which are simulated as a 5-state Markov model (Fig. 5b). The complete set of the kinetic rates for the spine head model are listed in Table 4. The release of Glutamate was modeled as a step pulse of duration 1 ms and amplitude 1 mM. For the voltage dependence of the NMDA receptors we assumed voltage clamp conditions, which results in a fixed rate of ion influx through open receptors. For all the deterministic equations a forward Euler with time step 10−9 s was sufficient and also efficient compared with other methods of higher order (e.g. Runga-Kutta fourth order).

Table 4
Complete list of the parameters of the spine model

5.2 Application of the Hybrid Algorithm for single release of Glutamate

In this section we apply our hybrid algorithm to the model of the postsynaptic spine head during a single release of Glutamate. The duration of the release was 1 msec and the amplitude 1 mM. In Fig. 6a we report one sample of the total number of NMDA receptors in the open state for one simulation run using the SSA, and in Fig. 6b we see the calcium concentration for the same simulation run. Inspection of Fig. 6b shows that rapid changes of the calcium transient follows the changes in the total open NMDA receptors. If there is no change in the number of open receptors the calcium concentration decays relatively smoothly as we can see after the last NMDA transition. This observation implies that in our model the main source of variability in the calcium dynamics is the NMDA receptors. This is due to the slow transition of the NMDA receptors in contrast to the fast reactions and diffusion of the calcium ions. This range of timescales exhibited in the postsynaptic spine head provides us the opportunity to test our adaptive partitioning, and see how well it overcomes possible bottlenecks while still capturing the stochastic behavior of the model.

Fig. 6
Individual stochastic simulations using the SSA. (a) Number of open NMDA receptors for 400 msec. The total number of NMDA receptors was 10. (b) Calcium concentration for the same simulation run. Rapid changes of calcium follows after changes of the open ...

For this system the first and second moment is sufficient for evaluating results obtained from our hybrid algorithm. Fig 7a shows the average number of NMDA receptors in the open state using the SSA (solid line) and our hybrid algorithm for 5000 trials. We performed simulations using two sets of parameters for partitioning the reactions: (i) μ = 500, λ = 20 (dotted line) and (ii) μ = 250, λ = 20 (dashed line). What is clear for both sets of parameters is that our results are almost identical with those obtained form the SSA. Not only do the means agree, but our results are also consistent with the SSA in respect to the variance of NMDAR transitions to the open state (Fig. 7b).

Fig. 7
Results from 5000 runs for the SSA (solid line) and for our hybrid approach. Two sets of parameters were used for the partitioning function, μ = 500, λ = 20 (dotted line) and μ = 250, λ = 20 (dashed line). (a) Average number ...

Using the same set of simulation parameters we also compared our approach with the SSA by simulating the concentration of Calcium. At T=0 the concentration of calcium was set to zero, and the concentration of Buffer in the spine head was 143 μM and 80 μM in the neck, a concentration high enough to avoid saturation. The concentration of pumps was uniformly distributed on the spine head and modeled with a fixed extrusion rate, as described previously. Because of the finite element approach we were able to study the effect of the spine geometry on the spatial variability of calcium transients. For instance, it is important to know the difference in the concentration or variance of calcium at the top and bottom portion of the spine head as this could impact the activation of enzymes which are important for induction of synaptic plasticity. The correlation of the spine morphology with any spatial gradients of calcium ions and plasticity is certainly an interesting issue but this in not the goal of the current study; therefore for testing our results we have used the average concentration of calcium in the whole spine head. Fig. 7c illustrates the mean concentration of calcium and Fig 7d shows the variance. We see that for a relatively high half maximum of the partitioning function (dashed line) we achieve good agreement with the SSA (solid line), and given a lower half maximum (dotted line) both the mean and variance are almost identical with the SSA.

In a recent study, off-lattice Kinetic Monte Carlo simulations of the calcium dynamics in the postsynaptic spine head required up to 60 hours for 1 s of real time ([Keller et al., 2008]). From our simulations results we noticed that the SSA requires time steps of the order of 10−11 s whereas the numerical integration of the PDEs 10−9 s. Therefore we could say that stochastic simulation of calcium dynamics, even for a simplified model of the spine head, is indeed a computationally demanding task. To demonstrate the performance increase achieved by our hybrid algorithm we measured the average time spent executing one simulation run. However, we need to mention that the simulation time depends on many parameters such as the numerical method used for solving deterministic equations, the random number generator used since this may require up to a hundred CPU cycles and the optimization techniques used for implementing the algorithm as well as the compiler that is used. In our case, each of our simulations were implemented in Matlab, and it is likely that our performance increase could be improved even more by implementing our simulations in a more efficient programming language. The platform we used for running the simulations is a cluster of 15 nodes with 2 dual core AMD Opteron 2.4 GHz and 16 GB RAM per node, running Red Hat Enterprise 3 Linux. In Fig. 8 we see the average simulation time for different sets of the parameters μ and λ as well as the two extremes, the SSA and the deterministic. We notice that the simulation time has been reduced significantly due to the assignment of some reactions to the continuous time regime. For each type of reaction we also recorded the average number of occurrences (Table 5), and we see that only the diffusion jump events of calcium ions and their reaction with the calcium buffer have been reduced significantly, suggesting that the main bottleneck of the system is the large population of Ca+2 and their high diffusion coefficient.

Fig. 8
Comparison of the average simulation time required for the SSA and the hybrid approach
Table 5
Estimation of the number of reactions and diffusion events in the discrete regime

6 Conclusions

Simulating complex biological systems is an indispensable tool for advancing our understanding of highly dynamic intracellular events. Representing these chemical reaction systems with deterministic rate-based equations fails to capture the stochastic nature of the underlying events, so Gillespie introduced the SSA for stochastically simulating well stirred systems. In living cells however, the volume is not homogenous, and spatial variations in the concentration of reacting species need to be considered. While one would ideally utilize the SSA to simulate reaction-diffusion systems as realistically as possible, for complex systems where some reaction channels occur with high frequency the SSA algorithm becomes computationally costly and impractical. In addition, the coexistence of molecular species in small quantities and infrequent reactions do not satisfy the Langevin approach. One way to overcome this difficulty, which arises in stiff systems, is the separation of time scales. In the past several methods have been developed in order to accelerate stochastic simulations of biochemical networks by partitioning the reactions into these two subsets and simulating the Slow reactions with the SSA and the Fast reactions with an accelerated ([Puchalka and Kierzek, 2004, Rao et al., 2002]) or deterministic rate-based method ([Salis and Kaznessis, 2005, Haseltine and Rawlings, 2002, Alfonsi et al., 2005])

In this paper we have presented a new hybrid algorithm for simulating stochastic and deterministic reaction systems and reaction systems with diffusion events. We have discussed in detail our partitioning strategy and the coupling of continuous and stochastic time regimes using fractional propensity functions. Our hybrid method is an extension of the exact Gillespie method for the multivariate chemical master equation, but it is more computationally efficient when some reactions or diffusion jump events occur with high frequency. The applicability and accuracy of our approach for efficiently simulating stochastic events in stiff reaction or reaction-diffusion systems was demonstrated using two different models, both important in computational biology. We need to notice that, while the applicability of similar approaches have been demonstrated for several systems, still these techniques have not matured to the point where we can know when to apply these methods. Therefore we have focused on improving the performance of SSA.

The first system was a model of bacteriophage T7, which incorporates only chemical reactions and not diffusion events. The particular model incorporates both fast and slow reactions while achieving two steady states for the population of the templates. We demonstrated the existence of two steady states, and we quantified the accuracy of our results by comparing the Kologorov distance of the cdf for the template populations obtained from our approach with the SSA. We showed how our adaptive partitioning of the propensities using a smooth S-shaped function generates results that agree with those of the SSA. The second test system was the postsynaptic spine head which has previously served as a biological setting for modeling calcium dynamics ([Holcman and Schuss, 2005, Bhalla, 2004]). In our study the spine head and neck were implemented as a simplified compartmental model, and the stiffness of this model arises from the different time scales of slow reactions (e.g. NMDAR transitions) versus the fast reaction or jump diffusion events of calcium ions. We demonstrate that our hybrid algorithm can accurately describe calcium and NMDAR dynamics much more efficiently than the SSA.

One advantage of a hybrid stochastic approach is that it helps to determine the contribution of different model components to its stochastic behavior. In a previous study ([Shouval and Kalantzis, 2005]) we demonstrated how stochasticity affects plasticity curves under the assumption that the only source of variability in calcium transients was due to the NMDARs. Now, results from the present study support the idea that the main source of variability in calcium dynamics is the NMDARs, and this observation is an indicator of the validity of our previous assumption. A clear understanding of the main sources of stochasticity in a more realistic biophysical model of the spine head is essential for the derivation of stochastic plasticity learning models. Other possible sources of stochasticity depend on the presynaptic component of the synapse. However, further mathematical analysis and simulations are required to explicitly test this issue.

We want to stress that in our present model of the postsynaptic spine head, we studied only the mean and variance of the average concentration of the calcium dynamics. However for different geometries there may have spatial variability of the calcium with respect to the mean and the variance. In this study we do not address the question of how spine morphology affects the spatiotemporal calcium dynamics, but the model is capable of being expanded for the study of spatial gradients of not only calcium but of a variety of enzymes in the spine head as well. For instance the only source of calcium in our model is the NMDARs, but in the future, the effect of internal stores of calcium (i.e. IP3) on calcium dynamics and plasticity will be tested. Another area in which to expand our model is the NMDARs themselves. NMDARs are heterotetramers composed of two NR1 subunits and two NR2 subunits. There are four different subtypes of NR2 subunits (A, B, C and D), and the NR2A and B subunits have been shown to play an important in regulating synaptic plasticity. More specifically, during early stages of development the majority of NMDARs contain NR2B subunits and are localized at the spine head. Later in development, NR2B subunits are replaced by NR2A in an activity dependent manner. In this study we used a 5-state Markov model for the NMDAR, but one possible extension of our model would be the replacement of the Markov model of the NMDARs with one which takes into account the different kinetics of receptors composed of different subunits ([Erreger et al., 2005, Kampa et al., 2005]).

Finally, in the present study we clamped the postsynaptic membrane voltage, which results in a fixed rate of calcium influx through the NMDARs. Our algorithm can also be used for studying calcium dynamics under different protocols for inducing plasticity where the voltage is not fixed. Vargas-Caballero and Robinson ([Vargas-Caballero and Robinson, 2004]) have proposed a Markov model for NMDARs that takes into account the voltage dependence of the receptor’s conductance. Previous studies ([Haseltine and Rawlings, 2002]) have used the CLE for simulating the Fast subset of reactions. In our code we simulate stochastically the chemical reactions as well as the diffusion of the calcium ions. An extension of our method is the replacement of the mass-action rate equations with modified CLE which take into account diffusion jump events.


This work was funded by grant 2 P01- NS038310-06A2 from the National Institute of Health. The author would like to thank Harel Shouval, Yoshihisa Kubota for helpful discussions and Matt Swulius who provided many useful suggestions for improving the paper.

7 Appendix

Evaluation of the accuracy of a hybrid approach could be based on the usage of different metrics. One of them is the Kolmogorov distance, as discussed previously. Here we formulate two different metrics. The first one is the relative error of the mean of the ith species concentration with respect to the results obtained from the fully stochastic approach,


The Σ denotes the sum over all the voxels r and the brackets [left angle bracket][right angle bracket] denote the average over a large enough set of realizations. The second metric is the relative error of the second moment, the variance,


where the variance function is given by:


Those metrics can be used also as a guide for estimating the accuracy of simulation results for different set of the parameters μ and λ. However depending on the system different metrics (i.e. histogram distance) may be used.


Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.


  • Alfonsi A, Cances E, Turinici G, Ventura BD, Huisinga W. Exact simulation of hybrid stochastic and deterministic models for biochemical systems. ESAIM Proc. 2005;14:1–13.
  • Arkin A, Ross J, McAdams HH. Stochastic kinetic analysis of developmental pathway bifurcation in phage λ-infected escherichia coli cells. Genetics. 1998;149:16331648. [PubMed]
  • Baras F, Mansour MM. Reaction-diffusion master equation:a comparison with microscopic simulations. Phys Rev E. 1996;54:6139. [PubMed]
  • Bernstein D. Simulating mesoscopic reaction-diffusion systems using the gillespie algorithm. Phys Rev E. 2005;71:041103. [PubMed]
  • Bhalla US. Signaling in small subcellular volume ii. stochastic and diffusion effects on synaptic network properties. Biophys J. 2004;87:745–753. [PubMed]
  • Cao Y, Gillespie D, Petzold L. The slow scale stochastic simulation algorithm. J Chem Phys. 2005;122:014116. [PubMed]
  • Cao Y, Li H, Petzold LR. Efficient formulation of the stochastic simulation algorithm for chemically reacting systems. J Phys Chem. 2004;105:1876–1889. [PubMed]
  • Cao Y, Petzold L. Accuracy limitations and the measurement of errors in the stochastic simulation of chemically reacting systems. J Comp Physics. 2006;212:6–24.
  • Chatterjee A, Katsoulakis M, Vlachos D. Binomial distribution based tau -leap accelerated stochastic simulation. J Chem Phys. 2005;122:024112. [PubMed]
  • Chiam KH, Tan CM, Bhargava V, Rajagopal G. Hybrid simulations of stochastic reaction-diffusion processes for modeling intracellular signaling pathways. Phys Rev E. 2006;74:051910. [PubMed]
  • Erreger K, Shashank MD, Banke T, Wyllie D, Traynelis S. Subunit-specific gating controls rat nr1/nr2a and nr1/nr2b nmda channel kinetics and synaptic signaling profiles. J Physiol. 2005;563:345–358. [PubMed]
  • Franks KM, Bartol TM, Sejnowski TJ. A monte carlo model reveals independent signaling at central glutamatergic synapse neuromuscular junction. Biophys J. 2001;83:2333–2348. [PubMed]
  • Franks KM, Bartol TM, Sejnowski TJ. A monte carlo model reveals independent signaling at central glutamatergic synapses. Biophys J. 2002;83:2333–48. [PubMed]
  • Frenkel D, Smith B. Understanding molecular simulation:from algorithms to applications. Elsevier; San Diego: 1996.
  • Gibson MA, Bruck J. An efficient formulation of the stochastic simulation of chemical systems with many species and many channels. J Phys Chem. 2000;104:1876–1889.
  • Gillespie DT. A general method for numerically simulating the stochastic evolution of coupled chemical reactions. J Comput Phys. 1976;22:403–434.
  • Gillespie DT. Approximate accelerated stochastic simulation of chemically reacting systems. J Chem Phys. 2001;115:1716–1733.
  • Gold J, Bear M. A model of dendritic spine calcium concentration exploring possible bases for a sliding synaptic modification threshold. Proc Natl Acad Sci U S A. 1993;91:3941–3945. [PubMed]
  • Goutsias J. Quasiequilibrium approximation of fast reaction kinetics in stochastic biochemical systems. J Chem Phys. 2005;122:184102. [PubMed]
  • Griffith M, Courtney T, Peccoud J, Sanders WH. Dynamic partitioning for hybrid simulation of the bistable hiv-1 transactivation network. Bioinofrmatics. 2006;22:2782–2789. [PubMed]
  • Haseltine EL, Rawlings JB. Approximate simulation of coupled fast and slow reactions for stochastic chemical kinetics. J Chem Phys. 2002;117:6959–6969.
  • Holcman D, Schuss Z. Modelling calcium dynamics in dendritic spines. SIAM J Appl Math. 2005;65:1006–1026.
  • Holmes W. Is the function of dendritic spines to concentrate calcium? Brain Research. 1990;519:338–342. [PubMed]
  • Ichikawa K. Localization of activated ca2+/calmodulin-dependent protein kinase ii within a spine: modeling and computer simulation. Neurocomputing. 2004;58:443–448.
  • Kampa BM, Clements J, Jonas P, Stuart GJ. Kinetics of mg+2 unblock of nmda receptors implications for spike-timing dependent synaptic plasticity. ESAIM Proc. 2005;14:1–13. [PubMed]
  • Keller D, Franks K, Bartol T, Sejnowski T. Calmodulin activation by calcium transients in the postsynaptic density of dendritic spines. PLOS One. 2008;3(4):e2045. [PMC free article] [PubMed]
  • Kesmir C, De Boer RJ. A mathematical model on germinal center kinetics and termination. J Immunol. 1999;163:2463–2469. [PubMed]
  • Kiehl T, Mattheyses R, Simmons M. Hybrid simulation of cellular behavior. Bioinofrmatics. 2004;20:316–322. [PubMed]
  • Knipe DM, Howley M. Fields virology. Lippincott Williams and Wilkins; Philadelphia: 2001.
  • Manninen T, Linne M, Ruohonen K. Developing ito stochastic differential equation models for neuronal signal transduction pathways. Comp Biol Chem. 2006;30:280–291. [PubMed]
  • McColluma M, Peterson GD, Cox CD, Simpson ML, Samatova NF. The sorting direct method for stochastic simulation of biochemical systems with varying reaction execution behavior. J Comput Biol Chem. 2005;30:39–49. [PubMed]
  • Morishita Y, Aihara K. Noise-reduction through interaction in gene expression and biochemical reaction processes. J Theor Biol. 2004;228:315–325. [PubMed]
  • Naoki H, Sakumura Y, Ishii S. Local signaling with molecular diffusion as a decoder of calcium signals in synaptic plasticity. Mol Syst Biol. 2005;1:0027. [PMC free article] [PubMed]
  • Perelson AS. Modelling viral and immune system dynamics. Nature Rev Immun. 2002;2:28–36. [PubMed]
  • Poranen MM, Daugelavicius R, Bamford DH. Common principles in viral entry. Ann Rev Microbiol. 2002;56:521–538. [PubMed]
  • Pornillos O, Garrus JE, Sunquist WI. Mechanisms of enveloped rna virus budding. Trends Cell Biol. 2002;12:569–579. [PubMed]
  • Puchalka J, Kierzek AM. Bridging the gap between stochastic and deterministic regimes in the kinetic simulations of the biochemical reaction networks. Biophys J. 2004;86:1357–1372. [PubMed]
  • Rao CV, Wolf DM, Arkin AP. Control, exploitation and tolerance of intracellular noise. Nature. 2002;420:321–237. [PubMed]
  • Rathinam M, Petzold LR, Cao Y, Gillespie DT. Stiffness in stochastic chemically reacting systems:the implicit tau-leaping method. J Chem Phys. 2003;119:12784–12794. [PubMed]
  • Risken H. The fokker-planck equation:methods of solution and applications. Springer; Berlin: 1989.
  • Rudiger S, Shuai JW, Huisinga W, Nagaiah C, Wamecke G, Parker I, Falcke M. Hybrid stochastic and deterministic simulations of calcium blips. Biophys J. 2007;93:1847–1857. [PubMed]
  • Salis H, Kaznessis Y. Accurate hybrid stochastic simulation of a system of coupled chemical or biochemical reactions. J Chem Phys. 2005;122:54103. [PubMed]
  • Shouval HZ, Kalantzis G. Stochastic properties of synaptic transmission affect the shape of spike time dependent plasticity curves. J Neurophys. 2005;93:1069–1073. [PubMed]
  • Smith AE, Helenius A. How viruses enter animal cells. Science. 2001;304:237–242. [PubMed]
  • Takahashi K, Kaizu K, Hu B, Tomita M. A multi-algorithm multi-timescale method for cell simulation. Bioinformatics. 2004;20(4):538546. [PubMed]
  • Tuckwell HC, Wan FYM. On the behaviour of solutions in viral dynamical models. BioSystems. 2004;73:157–161. [PubMed]
  • Vargas-Caballero M, Robinson HP. Fast and slow voltage-dependent dynamics of magnesium block in the nmda receptor: the asymmetric trapping block model. J Neurosci. 2004;24:6171–6180. [PubMed]
  • Volfobsky N, Parnas H, Segal M, Korkotian E. Geometry of dendritic spines affects calcium dynamics in hippocampal neurons: Theory and experiments. J Neurophys. 1999;82:450–462. [PubMed]
  • Wodarz D, Nowak MA. Mathematical models of hiv pathogenesis and treatment. BioEssays. 2002;24:1178–1187. [PubMed]
  • Wylie D, Hori Y, Dinner A, Chakraborty AK. A hybrid deterministic-stochastic algorithm for modelling cell signaling dynamics in spatially ingomogeneous environment and under the influence of external fields. J Phys Chem. 2006;110:12749–12765. [PubMed]
  • Zucker RS. Calcium and activity-dependent synaptic plasticity. Curr Opin Neurobiol. 1999;9:305–313. [PubMed]