Automatic Simulation of the Chemical Langevin Equation

Biochemical systems have important practical applications, in particular to understanding critical intra-cellular processes. Often biochemical kinetic models represent cellular processes as systems of chemical reactions, traditionally modeled by the deterministic reaction rate equations. In the cellular environment, many biological processes are inherently stochastic. The stochastic fluctuations due to the presence of some low molecular populations may have a great impact on the biochemical system behavior. Then, stochastic models are required for an accurate description of the system dynamics. An important stochastic model of biochemical kinetics is the Chemical Langevin Equation. In this work, we provide a numerical method for approximating the solution of the Chemical Langevin Equation, namely the derivative-free Milstein scheme. The method is compared with the widely used strategy for this class of problems, the Milstein method. As opposed to the Milstein scheme, the proposed strategy has the advantage that it does not require the calculation of exact derivatives, while having the same strong order of accuracy as the Milstein scheme. Therefore it may be used for an automatic simulation of the numerical solution of the Chemical Langevin Equation. The tests on several models of practical interest show that our method performs very well.


Introduction
A fundamental problem in the post-genomic biology is to describe and analyze the complex dynamical interactions which take place at the level of a single cell.Recent experimental techniques made it possible to study gene regulatory networks in living cells [1] as well as to generate synthetic gene networks [2].
There are currently several levels of refinement used for modeling the cellular dynamics.Often chemical kinetic models represent cellular processes as systems of chemical reactions.Traditionally, these processes were modeled as continuous deterministic systems, by reaction rate equations.However, the random fluctuations which are captured by the experiments [3][4][5][6] are neglected by such models.These fluctuations are due to low molecular numbers of some biochemical species.Then, stochastic models are required for an accurate description of the system dynamics.A stochastic model of the well-stirred biochemical systems is the Chemical Master Equation [7].Various algorithms proposed for the exact simulation of the solution of the Chemical Master Equation [8,9] are computationally very expensive for most practical applications.Approximate algorithms were designed and analyzed in the literature to speed-up the simulation for biochemical systems modeled with the Chemical Master Equation [10][11][12][13][14]. Nonetheless, more sophisticated techniques are necessary for dealing with systems which manifest stiffness.Stiffness is due to the presence of the multiple time-scales in the system, as some reactions are much faster than others.
As an intermediate model between the Chemical Master Equation and the reaction rate equations, the Chemical Langevin Equation (CLE) [15] is considered a very attractive choice in modeling many important biological processes.CLE consists of a system of stochastic differential equations, nonlinear and with non-commutative multiplicative noise.Most biochemical systems of interest typically involve many components interconnected in a complex manner.Thus, it is important to have efficient and accurate algorithms for simulating their mathematical models and in particular if they are stiff.However, the construction of algorithms to simulate and approximate the solution to these mathematical models is a challenging task, and research in this field is only at the initial stages [16][17][18].One of the widely used numerical methods to simulate the Chemical Langevin Equation is the Milstein scheme [19,20].This scheme has strong order of accuracy one, however it necessitates the calculation of some exact derivatives.This is a drawback of the Milstein strategy.
This paper provides a derivative-free numerical method for the strong approximation of the solution of the Chemical Langevin Equation.To our knowledge, the derivative-free Milstein scheme was not utilized before in the simulation of stochastic models of biochemical kinetics.The advantages of this method include: it is of strong order of accuracy one and it does not entail the calculation of exact derivatives.The derivative-free Milstein strategy estimates the derivatives by using finite differences.The proposed method may therefore be used for designing automatic simulation algorithms for generic models of well-stirred biochemical systems, in the Langevin regime.
The paper is organized as follows.Section 2 gives an introduction to the strong numerical solution of Itô stochastic differential equations.In Section 3 we discuss a stochastic continuous model of well-stirred biochemical kinetics, namely the Chemical Langevin Equation.Section 4 presents our proposed numerical strategy for the Chemical Langevin Equation.Numerical tests on several models of practical interest, showing the accuracy of the method provided, are given in Section 5. Finally, we summarize our conclusions in Section 6.

Background
A brief introduction to the numerical solution of Itô stochastic differential equations (SDE), which are essential in stochastic biochemical kinetic modelling, is presented below.The Itô formulation of an SDE system is where X is an N-dimensional stochastic process.Here are N-dimensional and represent the drift and the diffusion coefficients, respectively, while denotes an M-dimensional Wiener process with independent components.

 , M W
We are interested in the strong numerical solution of SDE.Strong numerical approximations are employed when an accurate approximation of the solution of an SDE on individual trajectories is desired, while weak numerical approximations are utilized when the approximation of the moments of the exact solution is sufficient.
Let X L be the numerical approximation on [0, T], after L steps with stepsize h T L  , of the exact solution of (1) and let   be a constant.

Strong Convergence
The approximation of is said to have strong order of convergence γ if there exists a constant , independent of h and for any

Weak Convergence
The approximation   of is said to have weak order of convergence γ if, for any polynomial P there exists a constant , independent of h and for any Here   E  denotes the expectation of a random variable and  a norm of an N-dimensional vector.
The focus of this work is on SDE with non-commutative noise [20], as the Chemical Langevin Equation has multiplicative non-commutative noise.For this class of problems, to obtain numerical methods of strong order of accuracy 1 on each interval   , t t h  , in addition to the computation of the Wiener increments , the simulation of the stochastic double Itô integrals ij I is necessary or, equivalently, of the Levy areas.The double Itô integrals ij I are defined as  .The double Itô integrals are estimated in terms of their Fourier series expansion truncated after p terms (see also [20]): , where, for any 1 , 2π are independent normally distributed with mean 0 and variance 1, for any 1  and any . Numerical experiments in the literature indicate that is sufficient for an accurate approximation , 5 p i j I of the double Itô integrals ij I .In our simulations we choose p = 5.

Milstein Method
The classical strong order 1 numerical method due to Milstein is used in the literature for approximating the exact solution of the Chemical Langevin Equation [19].
The Milstein scheme on the time interval   , , , where the Wiener increments are denoted by The differential operator j L is defined as

Derivative-Free Milstein Method
The strong order 1 Milstein strategy has the disadvantage that it requires derivative calculations, an issue for generating automatic simulation algorithms.The derivativefree Milstein schemes [21] overcome this difficulty.The derivative-free Milstein strategy for the general SDE (1), driven by M independent Wiener processes can be obtained from the Milstein method by replacing the derivatives by finite differences.Note that these differences require intermediate approximations at other points.The derivative-free Milstein scheme can be written as where the intermediate values are for . 1 i M  

Stochastic Continuous Models of Biochemical Kinetics
It has recently been acknowledged that stochastic models are more accurate than their deterministic counterparts for representing cellular dynamics.Biological processes at the single cell level are often modeled as systems of biochemical reactions.Below we discuss a key stochastic model of well-stirred biochemical kinetics.This model is valid for isothermal biochemically reacting systems with relatively large molecular numbers, in a constant volume.
Assume that N biochemical species 1 , , N S S  undergo M reaction channels 1 , , M R R  . The well-stirred assumption leads to a simplification of the molecular dynamics model.Under the above assumptions, the system state can be represented by a stochastic process . The components of the dynamical state vector are , the number of molecules of the species present in the system at time t, for any .
is the stoichiometric matrix of the biochemical system.The propensity   j a x of the reaction j R is defined as  d j a x t is the probability that a single reaction j R occurs in the time-interval   , d t t t  , given the state x at time t.The existence of the propensity function is a consequence of the kinetic theory.
A unimolecular reaction  and by the propensity  (the reaction being called dimerization).Assume that the system has a macroscopic time-scale.More precisely, we assume that a time step h exists satisfying simultaneously the conditions.1) h is small enough such that no propensity varies significantly in the interval   for t s t h    and each 1 .j M   2) and h is large enough such that each reaction j R occurs many times in the time-interval   The conditions 1) and 2) are satisfied for biochemical systems with abundant molecular numbers.Then, the dynamical state of the system may be approximated by a continuous Markov process satisfying Here  

Derivative-Free Simulation of the Chemical Langevin Equation
In this paper we propose to utilize the derivative-free Milstein strategy for simulating the Chemical Langevin Equation.The derivative-free Milstein scheme has strong order of accuracy 1 as the Milstein method, but achieves it without making use of exact derivatives.The computation of the exact derivative required by the Milstein technique constitutes a difficulty for designing automatic simulating algorithms for the CLE, as it necessitates the user's input of the expression of the exact derivative.The method we propose for the CLE avoids this problem.The Chemical Langevin Equation ( 10) is a particular case of the SDE (1), with the drift coefficient and the diffusion coefficients for . 1 j M   The derivative-free Milstein method applied to the CLE ( 10) is derived by substituting the drift and diffusion coefficients ( 12) and ( 13), respectively, in the scheme (6) to get on the time-interval   , n n t t h  .The approximations at the intermediate points are for . 1 i M   In the next section we test this method by comparing it to the Milstein scheme for the CLE as well as with Gillespie's algorithm [8,9], on some models of biochemical systems of interest in applications.

Numerical Experiments
Below are presented numerical tests of our proposed method for approximating the solution of the Chemical Langevin Equation of a generic model of biochemical systems.The simulations are performed in Matlab [22].The numerical strategy proposed is tested on several models of biochemically reaction systems arising in applications.We compare our method with the Milstein method, which is typically employed for simulating the CLE.However, there is no known biochemical system for which the Chemical Langevin Equation model has a closed form solution.To validate the accuracy of our numerical strategy we compare the histogram obtained with our method with the one computed with Gillespie's algorithm [8,9].Gillespie's algorithm is a Monte Carlo simulation strategy which generates trajectories in exact accordance with the probability distribution of the Chemical Master Equation.The Chemical Langevin Equation model is an approximation of the Chemical Master Equation model, valid in the regime of large molecular population numbers.While there is also a modeling error when comparing the histograms generated with our method for the CLE and with Gillespie's algorithm for the Chemical Master Equation, we note that the agreement of the numerical results is excellent, thus our scheme is shown to be very accurate.

Michaelis-Menten Model
Consider the Michaelis-Menten model [23], which deals with a very important mechanism of enzymatic catalysis.Four molecular species are involved in three reactions

S S S S S S S S S
The species S 1 is a substrate, S 2 is an enzyme, S 3 represents an enzyme-substrate complex, while S 4 is a product.The biochemical model shows how the enzyme transforms the substrate into a product.The reaction rate parameters are and 3 , while the propensity functions associated with the reactions (16) are The solution of the system ( 16) is subject to the initial conditions The simulation is performed on the timeinterval [0, 30].Finally, the state-change vectors are the columns of the following stoichiometric matrix The simulations with the derivative-free Milstein method with stepsizes h 1 = 10 −1 , h 2 = 5 × 10 −2 and h 3 = 10 −1 , with the Milstein scheme for the step h 1 = 10 −1 and with Gillespie's algorithm are presented in Figure 1.Each integration is performed over 10,000 trajectories.We note that our method has a similar computational cost with the Milstein scheme.We computed the ratio of the execution time of the Milstein method and that of our derivative-free Milstein technique.The value of this ratio for the sequence of steps above is between 0.98 and 1.01, showing that the two methods have almost the same computational cost on this model.Figure 1 presents the histograms at t = 30 for the species S 1 , S 2 , S 3 , and S 4 , respectively.The accuracy of our method is excellent on this model.

Stiff Biochemical Model
Below we illustrate our derivative-free scheme on a more complex system, a stiff non-linear biochemical model, consisting of the following reversible reactions

S S S S S S S S S S S S S S S S S
The reaction rate constants are and c 6 = 2.The system is integrated on the interval  with initial conditions .The reactions above are characterized by the propensities Also, the state-change vectors are the corresponding columns of the stoichiometric matrix We ran the simulations for the proposed derivativefree Milstein, the Milstein and the Gillespie algorithms over 10,000 trajectories.The system is stiff as several orders of magnitude separate some of the propensities.As in the case of deterministic problems, in stochastic systems stiffness poses challenges to the numerical simulation.Our proposed method can integrate the system efficiently and accurately.Figure 2 shows the histograms computed at time with the proposed derivative-free Milstein, Milstein and the Gillespie algorithms, respectively.We illustrate the behavior of our method applied with stepsize and compare it with the behavior of the Milstein strategy for the same step.The agreement between the results of the two numerical integrators is very good.We also study the be-   havior of our method for a sequence of time-steps h 1 = 10 −5 , h 2 = 5 × 10 −6 and h 3 = 10 −6 .As expected, the accuracy of our algorithm improves when the step is reduced.In addition, we computed the ratio of the execution times of the Milstein scheme and of the proposed tive-free Milstein method for the steps above and found that the ratio ranges between 1.003 and 1.031, showing a very similar computational cost of the two methods.Finally, for the given sequence of time-steps the derivative free Milstein histogram matches closely that of the Gillespie's algorithm, which shows that our method is very accurate.
The advantage of our derivative-free Milstein algorithm over the existing Milstein strategy is that ours can generate the numerical solution automatically, and does not require the user's input regarding the computation of the exact derivatives, as the Milstein scheme does.

Conclusion
In this work we described the derivative-free Milstein method for approximating the solution of the Chemical Langevin Equation.Chemical Langevin Equation is a key model of well-stirred biochemical systems, with many important practical applications.Many models arising in practice are mathematically stiff and therefore their simulation may be quite challenging.The method we discuss in this paper achieves strong order of accuracy one as the Milstein scheme, which is currently the most widely used simulation technique for the Chemical Langevin Equation.Unlike the Milstein scheme, the method we provided above does not require the computation of exact derivatives, which is a drawback of the Milstein technique.The tests on key biochemical models arising in applications show the excellent accuracy of our method.
Each reaction j R is completely characterized by its propensity and its state-change vector.The state-change vector of the reaction , j j R ν , is an N-dimensional vector with the component ij  being the variation in the number of molecules of the species produced by the firi S ing of one reaction j R .The matrix