# Adaptive Time-Stepping Using Control Theory for the Chemical Langevin Equation

Authors: Silvana Ilie; Monjur Morshed
Journal: Journal of Applied Mathematics
 (2015)

---

## Abstract

Stochastic modeling of biochemical systems has been the subject of intense research in recent years due to the large
number of important applications of these systems. A critical stochastic model of well-stirred biochemical systems in
the regime of relatively large molecular numbers, far from the thermodynamic limit, is the chemical Langevin equation.
This model is represented as a system of stochastic differential equations, with multiplicative and noncommutative
noise. Often biochemical systems in applications evolve on multiple time-scales; examples include slow transcription
and fast dimerization reactions. The existence of multiple time-scales leads to mathematical stiffness, which is a major
challenge for the numerical simulation. Consequently, there is a demand for efficient and accurate numerical methods to
approximate the solution of these models. In this paper, we design an adaptive time-stepping method, based on control
theory, for the numerical solution of the chemical Langevin equation. The underlying approximation method is the
Milstein scheme. The adaptive strategy is tested on several models of interest and is shown to have improved efficiency
and accuracy compared with the existing variable and constant-step methods.


---

## Body

## 1. Introduction

Stochastic modelling is essential for studying key biological processes, such as signaling chemical pathways in a cell, when some molecular species are in low numbers. The random fluctuations due to low amounts of certain biochemically reacting species have been observed experimentally [1–3]. Mathematically, the behaviour of such biochemical systems is accurately described in terms of Markov processes. For systems which may be assumed to be well-stirred, the dynamics of the system is governed by the chemical master equation [4]. While exact simulation algorithms for the chemical master equation exist in the literature [5, 6] they are typically quite intensive computationally, and therefore approximate simulation strategies were proposed [7–10]. However, the level of detail provided by the chemical master equation is often not necessary. In particular, when all reacting species are present in relatively large numbers, an approximate model, the chemical Langevin equation [11], is more computationally attractive to simulate that the chemical master equation.The chemical Langevin equation (CLE) is a stochastic differential equation (SDE) of dimension equal to the number of reacting biochemical species in the system. It has noncommutative multiplicative noise. Moreover, the biochemical systems arising in applications usually evolve on several time scales, meaning their models are mathematically stiff. Stiffness is a serious challenge for the numerical solution of deterministic models, and even more so of stochastic ones. Indeed, in the regions where the problem is stiff, the explicit numerical integrators are forced to drastically reduce the time-step to satisfy the accuracy criteria. This leads to a highly inefficient simulation when fixed step size algorithms are employed. By contrast, adaptive time-stepping strategies reduce the step size in the regions where stiffness is present such that the error is maintained below the required tolerance, but relax the time-step as soon as the integration exists these regions. This significantly reduces the computational cost of the simulation, while preserving the desired accuracy of the numerical solution. In the literature, the numerical solution of the chemical Langevin equation is generally computed using the Milstein scheme [12, 13]. The Milstein scheme is a strong order one numerical method [14]. Gaines and Lyons [15] proved that numerical techniques of at least strong order of accuracy one are required to guarantee that the numerical solution computed with adaptive time-stepping strategies converges to the exact solution for noncommutative SDE and in particular for the CLE.Techniques for adapting the time-step for the weak solution of stochastic differential equations were proposed by Szepessy et al. [16]. For the strong numerical solution of stochastic differential equations with multidimensional Wiener processes, most of the existing adaptive time-stepping strategies were designed for systems with commutative noise (see [17–19]). Gaines and Lyons [15] developed an adaptive scheme for the strong solution of SDE, but the scheme is quite restrictive. It is based on a Brownian tree structure and the only variations in the step size allowed are doubling or halving. We note that the chemical Langevin equation belongs to the class of noncommutative SDEs, which are more challenging to solve numerically. A variable time-stepping method for the chemical Langevin equation using the Brownian tree structure of Gaines and Lyons [15] was introduced by Sotiropoulos and Kaznessis [20] and showed to be expensive. Ilie and Teslya [21] developed an adaptive technique for the mean-square numerical solution of the chemical Langevin equation with small noise.This paper proposes a variable time-stepping strategy for the strong numerical solution of the chemical Langevin equation, based on proportional integral- (PI-) control. Similar to the existing adaptive methods [13, 20, 21], we employ the Milstein scheme to advance the integration. The technique uses estimates of the local error, based on the work by Sotiropoulos and Kaznessis [20]. This variable time-stepping method may be applied to any biochemical system which can be modelled with the chemical Langevin equation, having an arbitrary magnitude of the random fluctuations. To the best of our knowledge, this isthe first PI-controller for noncommutative stochastic differential equations. This approach extends to a class of noncommutative Itô SDE, the chemical Langevin equation model, and the work on PI-controllers for ODE by Söderlind [22, 23] and for commutative Stratonovich SDE by Burrage et al. [17]. In addition, the PI-control of the step size outperforms the integral- (I-) control considered by the previous works for the chemical Langevin equation [13, 21]. Finally, our strategy allows rejection of the time-step when the error is above the tolerance, while guaranteeing that the statistics of the numerical solution is not biased.An outline of the paper is given below. Section2 presents a stochastic continuous model of well-stirred biochemical kinetics, the chemical Langevin equation. Section 3 gives a brief description of the strong numerical solution of Itô stochastic differential equations. In Section 4, we design an adaptive numerical technique for the chemical Langevin equation. The advantages of the proposed variable step size method over the existing methods are illustrated on several models of biochemical systems of interest in applications, in Section 5. The conclusions are given in Section 6.


## 2. Chemical Langevin Equation

AssumeN biochemical species S1,…,SN participate in M reaction channels R1,…,RM. The biochemical system, held at constant temperature, is homogeneous. The system state at time t≥0 is denoted by X(t)=X1t,…,XNtT, where Xi(t) represents the number of Si molecules at time t, for any i=1,…,N. Mathematically, the system state is modelled as a Markov process. When one reaction Rj fires, the state of the system is updated using a vector νj called the state-change vector corresponding to the reaction Rj. This is an N-dimensional vector with entries νij, denoting the variation in the Si molecular species caused by one reaction Rj. Then, V={νij}1≤i≤N,1≤j≤M is the stoichiometric matrix of the biochemical system. To each reaction Rj it corresponds a propensity aj(x) defined as follows: for an infinitesimal time increment dt, aj(x)dt is the probability that, given the state x at time t, one reaction Rj happens during [t,t+dt).For a unimolecular reaction(1)Sk⟶cj⋯
the propensity function is aj(X)=cjXk, while for a bimolecular reaction
(2)Si+Sk⟶cj⋯
the propensity is of the form aj(X)=cjXiXk when i≠k and aj(X)=cjXi(Xi-1)/2 when i=k.The system state may be represented as(3)Xt=X0+∑j=1M‍νjPj∫0t‍ajXsds,
where Pj are independent unit rate Poisson processes [24].Assume that the leap condition is satisfied: there exists a time steph small enough such that
(4)ajXs≈constantintheintervalt,t+h,
for each 1≤j≤M. Under this assumption, the representation (3) may be approximated by
(5)Xt+h=Xt+∑j=1M‍νjPjajXt·h.
Equation (5), due to Gillespie, is known as the tau-leaping formula [8]. A further reduction is possible if, in addition, h may be chosen large enough such that each reaction Rj fires many times in the interval [t,t+h]. More precisely, h exists such that
(6)ajXth≫1
for any 1≤j≤M.A step sizeh exists for which conditions (4) and (6) are simultaneously satisfied when all molecular species are present in large amounts in the biochemical system. In this case, the Poisson variables in (5) may be approximated by normal random variables with the same mean and variance, since (6) holds. Consequently, the tau-leaping formula (5) yields
(7)Xt+h=Xt+∑j=1M‍νjajXth+ajXthNj0,1,
where Nj(0,1) are normal random variables with mean 0 and variance 1. Regarding the step-size h as an infinitesimal dt in (7) leads to the following equation:
(8)dXt=∑j=1M‍νjajXtdt+∑j=1M‍νjajXtdWjt,
where (Wj)1≤j≤M are independent Wiener processes. Equation (8) is called the chemical Langevin equation (CLE) [11]. It is a noncommutative Itô SDE with multiplicative noise. Note that the dynamical state of the system X(t) is represented in (8) as a continuous Markov process.


## 3. Numerical Methods for SDE

Consider a general system of Itô stochastic differential equations driven by anM-dimensional Wiener process, W=(W1,W2,…,WM), in the form
(9)dXt=ft,Xtdt+∑j=1M‍gjt,XtdWjt.
In the differential equation (9), X denotes an N-dimensional stochastic process and f(t,X(t)) and gi(t,X(t)) are N-dimensional drift and diffusion coefficients, respectively. The initial condition is
(10)X0=x0
for t=0.The SDE is assumed to have noncommutative noise, as is the case of the chemical Langevin equation model. If the differential operatorLj is defined as
(11)Lj=∑k=1N‍gk,j∂∂Xk,j=1,…,M
then the SDE is callednoncommutative if for some j1,j2=1,…,M, with j1≠j2,
(12)Lj1gj2≠Lj2gj1.
It is called acommutative SDE otherwise (see [14, p. 348]).Below, we discuss briefly the strong numerical solution of SDE, with a focus on a numerical method of strong order of accuracy one. The strong numerical solution of an SDE is computed when approximations of the exact solution of individual paths are of interest. By contrast, when the approximation of the moments of the exact solution are desired, weak numerical approximations are employed.The numerical approximation on[0,T] of the exact solution X(t) of (9), after L steps with step-size h=T/L, is denoted by XL. This approximation is said to havestrong order of convergenceγ>0 if there exists a constant C>0, independent of h and δ>0, such that the following inequality holds for any h∈(0,δ):
(13)EXtL-XL≤Chγ,
where ∥·∥ is some norm of a vector of dimension N and E(·) is the expectation of a random variable.Numerical methods of strong order of accuracy1 for a noncommutative SDE require, on each subinterval on each interval [t,t+h], the simulation of the Wiener increments
(14)ΔWj=Wjt+h-Wjt,1≤j≤M
and of either the double Itô integrals Iij [19],
(15)Iij=∫tt+h‍∫tt+s‍dWisdWjr,1≤i,j≤M,
or the Levy areas.Approximations of the double Itô integralsIi,j for 1≤i,j≤M, using the truncation after p terms of their Karhunen-Loève or Fourier series expansion [14, p. 198–203] (see also [19]), are computed as
(16)Ii,jp=12hξiξj-12haj,oξi-ai,oξj+hAi,jp.
For 1≤i,j≤M, and 1≤r≤p,
(17)ξj=ΔWjh,ζj,r,ηj,r~N0,1
are independent normally distributed random variables with mean 0 and variance 1. Here the following notation is used:
(18)Ai,jp=12π∑r=1p‍1rζi,rηj,r-ηi,rζj,r,aj,0=-1π2h∑r=1p‍1rζj,r-2hρpμj,p,ρp=112-12π2∑r=1p‍1r2.
An accurate simulation of the double Itô integrals requires a minimum value p=5 (see also [19]).Now, let us introduce a numerical method of strong order of accuracy1, typically used in the literature to simulate the CLE, namely, the Milstein scheme [14]. TheMilstein method on the time interval [tn,tn+h] computes
(19)Xn+1=Xn+ftn,Xnh+∑j=1M‍gjtn,XnΔWjnXn+1=+∑j1,j2=1M‍Lj1gj2tn,XnIj1,j2p,
where Lj is defined in (11) and the Wiener increments are ΔWjn=ΔWj(tn+h)-ΔWj(tn).


## 4. Variable Step-Size Control in the Simulation of the Chemical Langevin Equation

An adaptive time-stepping technique for the strong (pathwise) numerical solution of the chemical Langevin equation is proposed below. The underlying numerical technique is the Milstein scheme. Since the Milstein method is of strong order of accuracy one, numerical solutions employing it on variable time-step meshes converge to the exact solution of the chemical Langevin equation, as the step size converges to zero [15]. Recall that the CLE is a noncommutative SDE; hence the adaptive schemes developed in [17, 19] for commutative SDE do not apply. While the Brownian tree approach to adaptivity due to Gaines and Lyons [15] may be utilized, it has been shown to be quite expensive by Sotiropoulos and Kaznessis [20].For the proposed method, the sequence of time steps depends on the particular trajectory and is obtained using proportional-integral- (PI-) control [22]. Step rejections are allowed when the error is above the prescribed tolerance. When a step is rejected, we apply a strategy that guarantees that the correct Brownian path is followed [13, 19]. This strategy ensures that the numerical solution is not biased.


### 4.1. Milstein Scheme for the Chemical Langevin Equation

We will apply below a numerical technique, due to Milstein, to the stochastic model of well-stirred biochemical kinetics considered above. Remark that the chemical Langevin equation (8) is an Itô SDE of the form (9), with the drift coefficient of the form
(20)f(X)=∑j=1M‍ajXνj
and the diffusion coefficients
(21)gjX=ajXνj
for 1≤j≤M.Substituting the drift (20) and diffusion coefficients (21) in the numerical scheme (19) for a generic SDE leads to the Milstein method for the stochastic continuous model (8). The Milstein strategy employed to the CLE (8), on a time-interval [tn,tn+1], becomes
(22)Xkn+1=Xkn+∑j=1M‍νk,jajXnhnXkn+1=+∑j=1M‍νk,jajXnΔWjnXkn+1=+12∑j1=1M‍∑j2=1M‍∑l=1N‍νl,j1νk,j2∂aj2∂Xl(Xn)aj1(Xn)aj2(Xn)I(j1,j2),
for any k=1,…,N.In this paper we adhere to the usual practice of controlling the local error. The accuracy criterion is the local error on each Brownian path that, over each step, should be below the user-prescribed tolerance. Hence, the numerical integration may generate different time-step sequences on different paths. To apply the above error criterion, accurate estimates of the (pathwise) local error need to be computed. Lete(Xn,hn,ΔWn) be the estimation of the local error for the time interval [tn,tn+1], on some Brownian trajectory. The local error may be approximated by the sum of a drift part and a diffusion part. The drift term of the local error is estimated by
(23)edriftXn,hn=12dfdxfXn∞·hn2
for a drift f(X) given by (20). According to Sotiropoulos and Kaznessis [20], one can estimate the diffusion part of the pathwise local error produced by the Milstein scheme as
(24)ediffusionXn,hn,ΔWn=ΔWn,36GXn∞ediffusionXn,hn,ΔWn=×JGGXn∞
with
(25)GX=∑j=1M‍aj(X)νj
and JGG being the Jacobian of G. Here ΔWn,(3) denotes a column vector with the jth entry being (ΔWjn)3.


### 4.2.  Adaptive PI-Control for CLE

The local error committed by the discretization method, the Milstein scheme (22), for the CLE is required to be smaller than a prescribed tolerance, ε, at each step and on each Brownian path; that is,
(26)eXn,hn,ΔWn≤ε.Before describing the proposed adaptive time-stepping strategy, let us discuss the important problem of step rejection in the numerical integration of an SDE, when the desired accuracy is not achieved. For SDE, a step rejection must be performed such that the statistics of the approximate solution is not biased. More precisely, when the numerical solution is advanced on a particular Brownian path and a step is rejected as it failed to satisfy the accuracy requirement, the subsequent steps must be chosen such that the same Brownian path is maintained. In this work, this is achieved using a Brownian bridge [18]. If the step h is tried and rejected on the current interval, [t,t+h], its corresponding Wiener increments,
(27)ΔWh=Wt+h-Wt~h·N0,1,
were sampled. Denote the sampled value ΔWh=wh. Then, a smaller time-step is tried, τ<h, and the corresponding Wiener increments on the subintervals [t,t+τ] and [t+τ,t+h] are created
(28)ΔWτ=Wt+τ-Wt,ΔWh-τ=Wt+h-Wt+τ.
To ensure that the integration preserves the already created Brownian path, a Brownian bridge is employed:
(29)ΔWτ=τhwh+τh-τh·η,ΔWh-τ=h-τhwh-τh-τh·η.
Here η denotes a new normally distributed random variable with mean 0 and variance 1.One step-size change method is the integral- (I-) controller for predicting the future step,hn+1, based on the current step hn and the current local error e(Xn,hn,ΔWn). TheI-controller is
(30)hn+1=hnθ·εeXn,hn,ΔWn1/K,
for some constant K. The safety factor θ≤1 is introduced to reduce the number of rejected steps. The values considered in our simulations are θ=0.8 and K=2. Nonetheless, some step rejects may occur since the local error depends on the particular Brownian path; hence the principal error function has random values.To reduce the number of step rejections, we propose a proportional integral- (PI-) controller for varying the step size in the numerical CLE solving. Based on control theory, these PI-controllers present several advantages over the classical I-controllers, including improved efficiency and computational stability of the numerical solution. For more details on the control theory approach to adaptive time-stepping, we refer the reader to [22].ThePI-controller is given by
(31)hn+1=hnθ·εeXn,hn,ΔWnkIhn+1=×eXn-1,hn-1,ΔWn-1eXn,hn,ΔWnkP.
For a further reduction of number of step rejections, the adaptive algorithm employs a more conservative approach to controlling the time-step:
(32)hn+1=hnmin⁡fmax⁡,×eXn-1,hn-1,ΔWn-1eXn,hn,ΔWnkPhn+1=hnmin⁡fmax⁡fmin⁡,θ·εeXn,hn,ΔWnkI×eXn-1,hn-1,ΔWn-1eXn,hn,ΔWnkPhn+1=hnminaxfmax⁡,×eXn-1,hn-1,ΔWn-1eXn,hn,ΔWnkP,
which is similar to an adaptive procedure for the numerical solution of ODE [25]. Note that the factors fmax⁡>1 and fmin⁡<1 ensure that the next step does not increase or decrease too much, so that fewer steps are rejected. A similar conservative technique is applied to the I-controller (30).The characteristic equation of the PI-controller (31) is
(33)q2-1-kkI-kkpq-kkp=0.
Here k is the order of the local error, e(Xn,hn,ΔWn)=ϕnhnk, where ϕn is the principal error function. For the PI-controller to be stable, the roots of its characteristic polynomial should be located inside the unit circle [22]. The parameters kI, referred to as the integral gain, and kP, known as the proportional gain, are chosen to obtain the desired properties of the controller.Numerical testing on many problems showed that the following parameter values of the PI-controller give excellent results for mildly stiff models of biochemical systems:(kkI,kkP)=(0.1,0.2). Remark that the controller is stable, since its roots are within the unit circle. The PI-controllers for ODE, which were designed to produce smooth sequences of time-steps, do not perform well on SDE, where the goals are to increase the accuracy and reduce the number of rejected steps.


## 4.1. Milstein Scheme for the Chemical Langevin Equation

We will apply below a numerical technique, due to Milstein, to the stochastic model of well-stirred biochemical kinetics considered above. Remark that the chemical Langevin equation (8) is an Itô SDE of the form (9), with the drift coefficient of the form
(20)f(X)=∑j=1M‍ajXνj
and the diffusion coefficients
(21)gjX=ajXνj
for 1≤j≤M.Substituting the drift (20) and diffusion coefficients (21) in the numerical scheme (19) for a generic SDE leads to the Milstein method for the stochastic continuous model (8). The Milstein strategy employed to the CLE (8), on a time-interval [tn,tn+1], becomes
(22)Xkn+1=Xkn+∑j=1M‍νk,jajXnhnXkn+1=+∑j=1M‍νk,jajXnΔWjnXkn+1=+12∑j1=1M‍∑j2=1M‍∑l=1N‍νl,j1νk,j2∂aj2∂Xl(Xn)aj1(Xn)aj2(Xn)I(j1,j2),
for any k=1,…,N.In this paper we adhere to the usual practice of controlling the local error. The accuracy criterion is the local error on each Brownian path that, over each step, should be below the user-prescribed tolerance. Hence, the numerical integration may generate different time-step sequences on different paths. To apply the above error criterion, accurate estimates of the (pathwise) local error need to be computed. Lete(Xn,hn,ΔWn) be the estimation of the local error for the time interval [tn,tn+1], on some Brownian trajectory. The local error may be approximated by the sum of a drift part and a diffusion part. The drift term of the local error is estimated by
(23)edriftXn,hn=12dfdxfXn∞·hn2
for a drift f(X) given by (20). According to Sotiropoulos and Kaznessis [20], one can estimate the diffusion part of the pathwise local error produced by the Milstein scheme as
(24)ediffusionXn,hn,ΔWn=ΔWn,36GXn∞ediffusionXn,hn,ΔWn=×JGGXn∞
with
(25)GX=∑j=1M‍aj(X)νj
and JGG being the Jacobian of G. Here ΔWn,(3) denotes a column vector with the jth entry being (ΔWjn)3.


## 4.2.  Adaptive PI-Control for CLE

The local error committed by the discretization method, the Milstein scheme (22), for the CLE is required to be smaller than a prescribed tolerance, ε, at each step and on each Brownian path; that is,
(26)eXn,hn,ΔWn≤ε.Before describing the proposed adaptive time-stepping strategy, let us discuss the important problem of step rejection in the numerical integration of an SDE, when the desired accuracy is not achieved. For SDE, a step rejection must be performed such that the statistics of the approximate solution is not biased. More precisely, when the numerical solution is advanced on a particular Brownian path and a step is rejected as it failed to satisfy the accuracy requirement, the subsequent steps must be chosen such that the same Brownian path is maintained. In this work, this is achieved using a Brownian bridge [18]. If the step h is tried and rejected on the current interval, [t,t+h], its corresponding Wiener increments,
(27)ΔWh=Wt+h-Wt~h·N0,1,
were sampled. Denote the sampled value ΔWh=wh. Then, a smaller time-step is tried, τ<h, and the corresponding Wiener increments on the subintervals [t,t+τ] and [t+τ,t+h] are created
(28)ΔWτ=Wt+τ-Wt,ΔWh-τ=Wt+h-Wt+τ.
To ensure that the integration preserves the already created Brownian path, a Brownian bridge is employed:
(29)ΔWτ=τhwh+τh-τh·η,ΔWh-τ=h-τhwh-τh-τh·η.
Here η denotes a new normally distributed random variable with mean 0 and variance 1.One step-size change method is the integral- (I-) controller for predicting the future step,hn+1, based on the current step hn and the current local error e(Xn,hn,ΔWn). TheI-controller is
(30)hn+1=hnθ·εeXn,hn,ΔWn1/K,
for some constant K. The safety factor θ≤1 is introduced to reduce the number of rejected steps. The values considered in our simulations are θ=0.8 and K=2. Nonetheless, some step rejects may occur since the local error depends on the particular Brownian path; hence the principal error function has random values.To reduce the number of step rejections, we propose a proportional integral- (PI-) controller for varying the step size in the numerical CLE solving. Based on control theory, these PI-controllers present several advantages over the classical I-controllers, including improved efficiency and computational stability of the numerical solution. For more details on the control theory approach to adaptive time-stepping, we refer the reader to [22].ThePI-controller is given by
(31)hn+1=hnθ·εeXn,hn,ΔWnkIhn+1=×eXn-1,hn-1,ΔWn-1eXn,hn,ΔWnkP.
For a further reduction of number of step rejections, the adaptive algorithm employs a more conservative approach to controlling the time-step:
(32)hn+1=hnmin⁡fmax⁡,×eXn-1,hn-1,ΔWn-1eXn,hn,ΔWnkPhn+1=hnmin⁡fmax⁡fmin⁡,θ·εeXn,hn,ΔWnkI×eXn-1,hn-1,ΔWn-1eXn,hn,ΔWnkPhn+1=hnminaxfmax⁡,×eXn-1,hn-1,ΔWn-1eXn,hn,ΔWnkP,
which is similar to an adaptive procedure for the numerical solution of ODE [25]. Note that the factors fmax⁡>1 and fmin⁡<1 ensure that the next step does not increase or decrease too much, so that fewer steps are rejected. A similar conservative technique is applied to the I-controller (30).The characteristic equation of the PI-controller (31) is
(33)q2-1-kkI-kkpq-kkp=0.
Here k is the order of the local error, e(Xn,hn,ΔWn)=ϕnhnk, where ϕn is the principal error function. For the PI-controller to be stable, the roots of its characteristic polynomial should be located inside the unit circle [22]. The parameters kI, referred to as the integral gain, and kP, known as the proportional gain, are chosen to obtain the desired properties of the controller.Numerical testing on many problems showed that the following parameter values of the PI-controller give excellent results for mildly stiff models of biochemical systems:(kkI,kkP)=(0.1,0.2). Remark that the controller is stable, since its roots are within the unit circle. The PI-controllers for ODE, which were designed to produce smooth sequences of time-steps, do not perform well on SDE, where the goals are to increase the accuracy and reduce the number of rejected steps.


## 5. Numerical Experiments

We test the adaptive time-stepping technique, based on control theory, developed above for solving numerically the CLE on three interesting models of biochemical systems. The performance of the PI-controllers is compared to that of the I-controller and of the existing constant-step methods. The comparison is done as follows: for each of the PI- or I-controller the algorithm is run with a certain tolerance and the number of attempted (i.e., of accepted and rejected) steps is stored. Then the fixed step size method is applied with the largest number of attempted steps, between the I- and PI-adaptive schemes, and its local error is measured. In our experiments, we choose the following values for the parameters:fmin⁡=0.5 and fmax⁡=2.In order to validate the accuracy of our variable step size strategy for the numerical solution of the CLE, we compare the histogram generated with our scheme for the CLE with that computed with the exact algorithm for the chemical master equation due to Gillespie [5, 6]. Although chemical Langevin equation approximates the chemical master equation (CME) model, and thus some error in modelling exists, we obtain a very good agreement between the histograms using Gillespie’s algorithm for the CME and our adaptive strategy for the CLE. This shows the excellent accuracy of the proposed method.


### 5.1. Stiff Biochemical Reaction Model

Consider the following stiff and nonlinear biochemical reaction system(34)S1+S2⟶c1S3,S3⟶c2S1+S2,S1+S3⟶c3S2,S2⟶c4S1+S3,S3+S2⟶c5S1,S1⟶c6S3+S2.
The propensities characterizing these reactions are
(35)a1X=c1X1X2,a2X=c2X3,a3X=c3X1X3,a4X=c4X2,a5X=c5X2X3,a6X=c6X1.
The stoichiometric matrix
(36)V=-11-111-1-111-1-111-1-11-11
has as columns the state-change vectors of the reactions above. The system parameters are as follows: c1=100, c2=104, c3=10-4, c4=10-2, c5=10, and c6=103 for the reaction rate constants, and X(0)=100,1000,100′ for the initial conditions. The interval of integration is [0,2×10-2].We simulated 10,000 trajectories with the integral (I), the proposed proportional integral- (PI-) adaptive and the fixed step size algorithms of the CLE for the sequence of tolerances1,10-1, and 10-2. The total number of attempted, accepted, and rejected steps as well as the corresponding errors is stored in Table 1. Also, we present the error generated with the fixed time-step method with the same total number of steps as the most expensive between I- and the PI-strategies, for the same tolerance. Remark that the PI-controller takes fewer attempted steps compared with the I-controller, for the same ε, being thus more efficient. The ratio of rejected to attempted steps for the PI- adaptive scheme ranges from 5% to 8% and it decreases as the tolerance decreases. For the I-adaptive technique the same ratio of rejected to attempted steps is much larger, between 20% and 30% for the tolerances tried. The fixed step size method produces a local error of up to 100 times larger than the adaptive methods, for the same computational cost.Table 1
The stiff model: the number of steps taken by the I-adaptive, PI-adaptive, and fixed step-size methods on the time-interval [0, 2 × 10−2].
ε
Method
Tried steps
Accepted steps
Failed steps
Error/ε
1
I-adapt.
1943
1357
585
0.9985
PI-adapt.
1800
1652
148
0.9971
Fixed
1943
—
—
56.7290
10−1
I-adapt.
5298
3910
1388
0.9995
PI-adapt.
4926
4587
339
0.9990
Fixed
5298
—
—
76.2772
10−2
I-adapt.
14590
11543
3046
0.9998
PI-adapt.
13816
13061
755
0.9997
Fixed
14590
—
—
100.6531To verify the accuracy of the proposed PI-controller for adapting the time-step when simulating the CLE, we compare the histograms at timet=2×10-2 obtained with the Milstein PI-adaptive scheme and the Gillespie algorithms, respectively. The histograms, plotted in Figure 1, show an excellent agreement between our method and the exact algorithm; consequently our scheme is very accurate.Figure 1
The stiff model: the histograms atT=2×10-2 for the species S1,S2, and S3 obtained with the SSA (green) and the PI-adaptive method for the CLE with tolerances ε=10-1 (blue) and ε=10-2 (red). The simulation uses 10,000 paths.


### 5.2. Schlögl Model

The Schlögl model [26] is well-known for its bistable behavior. While for the reaction rate equation model the trajectories converge to one of the two stable states, for the chemical Langevin equation, the trajectories may switch between the two stable states due to the noise in the system.The Schlögl model consists of the following reactions:(37)A+2X⟶c13X,3X⟶c2A+2X,B⟶c3X,X⟶c4B.
The reaction rate constants take the values c1=3×10-7, c2=10-4, c3=10-3, and c4=3.5. The state change vectors of the reactions in the Schlögl model are the corresponding columns of the stoichiometric matrix:
(38)V=1-11-1-110000-11.
The reaction propensities may be computed as
(39)a1X=c1AXX-12,a2X=c2XX-1X-26,a3X=c3B,a4X=c4X.
The numbers of molecules of the species A and B have constant values, A=105 and B=2×105. The initial condition for the molecular number of the species X is X(0)=250.The system is integrated on the time-interval[0,40]. The simulations of our PI-adaptive, the existing I-adaptive, and the constant-step schemes used for the CLE and of the SSA are done for 10,000 trajectories. The total number of attempted, accepted, and rejected steps and the measured error for the proposed PI-adaptive and the I-adaptive algorithms, as well as the error committed by the constant-step scheme, are depicted in Table 2. The simulations are performed for the tolerances ε=10-1, ε=10-2, and ε=10-3. The proposed variable time-stepping algorithm produces an error of up to 19 smaller than that of the fixed-step-size method. Furthermore, the ratio of the number of rejected steps to the total number of steps tried for our PI-controlled time-step method is very small, between 3% and 5%, while for the I-adaptive scheme ranges from 11% to 23%.Table 2
The Schlögl model: the number of steps taken by the I-adaptive, PI-adaptive, and fixed stepsize methods on the time-interval[0,40].
ε
Method
Tried steps
Accepted steps
Failed steps
Error/ε
10−1
I-adapt.
1235
950
285
0.9976
PI-adapt.
1198
1123
75
0.9942
Fixed
1235
—
—
10.7990
10−2
I-adapt.
3472
2869
603
0.9992
PI-adapt.
3334
3179
155
0.9981
Fixed
3472
—
—
14.6088
10−3
I-adapt.
9724
8617
1107
0.9997
PI-adapt.
9587
9275
312
0.9993
Fixed
9724
—
—
18.7245In Figure2 we show the evolution in time of the error scaled by the tolerance, on a sample path, for the tolerance 10-3 for both our variable time-step and the existing fixed-step schemes. Moreover, we examined the accuracy of the numerical solution obtained with our adaptive method in Figure 3. We compared the probability distribution at time t=40 for the species X computed with the PI-adaptive method for the CLE with ε=10-2 and ε=10-3 and with the SSA. The histograms agree very accurately.Figure 2
The Schlögl model: the local error/tol versus time for the adaptive and fixed step-size algorithms on a sample trajectory, withtol=10-3 on [0,40].Figure 3
The Schlögl model: the histogram at timeT=40 of the species X obtained with the SSA (green) and the PI-adaptive method for the CLE with tolerances ε=10-2 (blue) and ε=10-3 (red). The simulation uses 10,000 paths.


### 5.3. Decay-Dimerization Model

The final system used to investigate the behaviour of our variable step-size algorithm is the decay-dimerization model [8]. This model consists of three molecular species involved in four reactions:
(40)S1⟶c10,S1+S1⟶c2S2,S2⟶c3S1+S1,S2⟶c4S3.
The propensity functions corresponding to these reactions are
(41)a1X=c1X1,a2X=c2X1X1-12,a3X=c3X2,a4X=c4X2,
with the reaction rate parameter values c1=1, c2=0.002, c3=0.5, and c4=0.04. The solution of this biochemical reaction model is subject to the initial conditions X1(0)=X2(0)=104 and X3(0)=100. The integration is performed on the time-interval [0,5]. In addition, the stoichiometric matrix for this system is
(42)V=-1-22001-1-10001.To analyze the efficiency of the proposed method, we simulated our PI-adaptive, the I-adaptive, and the existing constant-step algorithms on 10,000 paths. Table3 reports the numerical results: the number of attempted, accepted, and rejected steps of each of the variable step-size strategies, for tolerances ε=10-1,10-2, and 10-3, respectively. In addition, Table 3 records the local errors for each of the above schemes, for a given tolerance. Moreover, it shows the local error obtained using the fixed time-step method with the same total number of attempted time-steps as the most expensive technique, the integral- (I-) controller one (being thus of similar computational cost with it). Observe that the I-controller attempts more steps than our PI-controller, for the same tolerance; therefore it is more expensive. Likewise, the I-adaptive technique has a ratio of rejected to attempted steps between 6% and 14%, while the same step rejection ratio for our PI-adaptive scheme is below 3% and is decreasing with the decrease of tolerance. It is important to remark that the existing fixed step-size method yields a local error which is up to 780 times larger than the proposed variable time-stepping method. Consequently, our variable step-size algorithm is much more accurate than the constant time-step method, for a similar computational cost.Table 3
The decay-dimerization model: the number of steps taken by the I-adaptive, PI-adaptive, and fixed stepsize methods on the time-interval[0,5].
ε
Method
Tried steps
Accepted steps
Failed steps
Error/ε
10−1
I-adapt.
1344
1150
194
0.9980
PI-adapt.
1321
1275
46
0.9952
Fixed
1344
—
—
573.9135
10−2
I-adapt.
3888
3505
383
0.9992
PI-adapt.
3852
3752
100
0.9984
Fixed
3888
—
—
685.7963
10−3
I-adapt.
11509
10831
678
0.9996
PI-adapt.
11496
11298
198
0.9994
Fixed
11509
—
—
782.6674To assess the accuracy of the proposed variable step size method, we perform simulations for 10,000 trajectories with the PI-adaptive Milstein method for the CLE with tolerances10-2 and 10-3 and with Gillespie’s algorithm. The histograms generated with these algorithms at time t=5 for each of the species are plotted in Figure 4. Again, the accuracy of our variable step size strategy is excellent.Figure 4
The decay-dimerization model: the histograms atT=5 for the species S1,S2, and S3 obtained with the SSA (green) and the PI-adaptive method for the CLE with tolerances ε=10-2 (blue) and ε=10-3 (red). The simulation uses 10,000 paths.


## 5.1. Stiff Biochemical Reaction Model

Consider the following stiff and nonlinear biochemical reaction system(34)S1+S2⟶c1S3,S3⟶c2S1+S2,S1+S3⟶c3S2,S2⟶c4S1+S3,S3+S2⟶c5S1,S1⟶c6S3+S2.
The propensities characterizing these reactions are
(35)a1X=c1X1X2,a2X=c2X3,a3X=c3X1X3,a4X=c4X2,a5X=c5X2X3,a6X=c6X1.
The stoichiometric matrix
(36)V=-11-111-1-111-1-111-1-11-11
has as columns the state-change vectors of the reactions above. The system parameters are as follows: c1=100, c2=104, c3=10-4, c4=10-2, c5=10, and c6=103 for the reaction rate constants, and X(0)=100,1000,100′ for the initial conditions. The interval of integration is [0,2×10-2].We simulated 10,000 trajectories with the integral (I), the proposed proportional integral- (PI-) adaptive and the fixed step size algorithms of the CLE for the sequence of tolerances1,10-1, and 10-2. The total number of attempted, accepted, and rejected steps as well as the corresponding errors is stored in Table 1. Also, we present the error generated with the fixed time-step method with the same total number of steps as the most expensive between I- and the PI-strategies, for the same tolerance. Remark that the PI-controller takes fewer attempted steps compared with the I-controller, for the same ε, being thus more efficient. The ratio of rejected to attempted steps for the PI- adaptive scheme ranges from 5% to 8% and it decreases as the tolerance decreases. For the I-adaptive technique the same ratio of rejected to attempted steps is much larger, between 20% and 30% for the tolerances tried. The fixed step size method produces a local error of up to 100 times larger than the adaptive methods, for the same computational cost.Table 1
The stiff model: the number of steps taken by the I-adaptive, PI-adaptive, and fixed step-size methods on the time-interval [0, 2 × 10−2].
ε
Method
Tried steps
Accepted steps
Failed steps
Error/ε
1
I-adapt.
1943
1357
585
0.9985
PI-adapt.
1800
1652
148
0.9971
Fixed
1943
—
—
56.7290
10−1
I-adapt.
5298
3910
1388
0.9995
PI-adapt.
4926
4587
339
0.9990
Fixed
5298
—
—
76.2772
10−2
I-adapt.
14590
11543
3046
0.9998
PI-adapt.
13816
13061
755
0.9997
Fixed
14590
—
—
100.6531To verify the accuracy of the proposed PI-controller for adapting the time-step when simulating the CLE, we compare the histograms at timet=2×10-2 obtained with the Milstein PI-adaptive scheme and the Gillespie algorithms, respectively. The histograms, plotted in Figure 1, show an excellent agreement between our method and the exact algorithm; consequently our scheme is very accurate.Figure 1
The stiff model: the histograms atT=2×10-2 for the species S1,S2, and S3 obtained with the SSA (green) and the PI-adaptive method for the CLE with tolerances ε=10-1 (blue) and ε=10-2 (red). The simulation uses 10,000 paths.


## 5.2. Schlögl Model

The Schlögl model [26] is well-known for its bistable behavior. While for the reaction rate equation model the trajectories converge to one of the two stable states, for the chemical Langevin equation, the trajectories may switch between the two stable states due to the noise in the system.The Schlögl model consists of the following reactions:(37)A+2X⟶c13X,3X⟶c2A+2X,B⟶c3X,X⟶c4B.
The reaction rate constants take the values c1=3×10-7, c2=10-4, c3=10-3, and c4=3.5. The state change vectors of the reactions in the Schlögl model are the corresponding columns of the stoichiometric matrix:
(38)V=1-11-1-110000-11.
The reaction propensities may be computed as
(39)a1X=c1AXX-12,a2X=c2XX-1X-26,a3X=c3B,a4X=c4X.
The numbers of molecules of the species A and B have constant values, A=105 and B=2×105. The initial condition for the molecular number of the species X is X(0)=250.The system is integrated on the time-interval[0,40]. The simulations of our PI-adaptive, the existing I-adaptive, and the constant-step schemes used for the CLE and of the SSA are done for 10,000 trajectories. The total number of attempted, accepted, and rejected steps and the measured error for the proposed PI-adaptive and the I-adaptive algorithms, as well as the error committed by the constant-step scheme, are depicted in Table 2. The simulations are performed for the tolerances ε=10-1, ε=10-2, and ε=10-3. The proposed variable time-stepping algorithm produces an error of up to 19 smaller than that of the fixed-step-size method. Furthermore, the ratio of the number of rejected steps to the total number of steps tried for our PI-controlled time-step method is very small, between 3% and 5%, while for the I-adaptive scheme ranges from 11% to 23%.Table 2
The Schlögl model: the number of steps taken by the I-adaptive, PI-adaptive, and fixed stepsize methods on the time-interval[0,40].
ε
Method
Tried steps
Accepted steps
Failed steps
Error/ε
10−1
I-adapt.
1235
950
285
0.9976
PI-adapt.
1198
1123
75
0.9942
Fixed
1235
—
—
10.7990
10−2
I-adapt.
3472
2869
603
0.9992
PI-adapt.
3334
3179
155
0.9981
Fixed
3472
—
—
14.6088
10−3
I-adapt.
9724
8617
1107
0.9997
PI-adapt.
9587
9275
312
0.9993
Fixed
9724
—
—
18.7245In Figure2 we show the evolution in time of the error scaled by the tolerance, on a sample path, for the tolerance 10-3 for both our variable time-step and the existing fixed-step schemes. Moreover, we examined the accuracy of the numerical solution obtained with our adaptive method in Figure 3. We compared the probability distribution at time t=40 for the species X computed with the PI-adaptive method for the CLE with ε=10-2 and ε=10-3 and with the SSA. The histograms agree very accurately.Figure 2
The Schlögl model: the local error/tol versus time for the adaptive and fixed step-size algorithms on a sample trajectory, withtol=10-3 on [0,40].Figure 3
The Schlögl model: the histogram at timeT=40 of the species X obtained with the SSA (green) and the PI-adaptive method for the CLE with tolerances ε=10-2 (blue) and ε=10-3 (red). The simulation uses 10,000 paths.


## 5.3. Decay-Dimerization Model

The final system used to investigate the behaviour of our variable step-size algorithm is the decay-dimerization model [8]. This model consists of three molecular species involved in four reactions:
(40)S1⟶c10,S1+S1⟶c2S2,S2⟶c3S1+S1,S2⟶c4S3.
The propensity functions corresponding to these reactions are
(41)a1X=c1X1,a2X=c2X1X1-12,a3X=c3X2,a4X=c4X2,
with the reaction rate parameter values c1=1, c2=0.002, c3=0.5, and c4=0.04. The solution of this biochemical reaction model is subject to the initial conditions X1(0)=X2(0)=104 and X3(0)=100. The integration is performed on the time-interval [0,5]. In addition, the stoichiometric matrix for this system is
(42)V=-1-22001-1-10001.To analyze the efficiency of the proposed method, we simulated our PI-adaptive, the I-adaptive, and the existing constant-step algorithms on 10,000 paths. Table3 reports the numerical results: the number of attempted, accepted, and rejected steps of each of the variable step-size strategies, for tolerances ε=10-1,10-2, and 10-3, respectively. In addition, Table 3 records the local errors for each of the above schemes, for a given tolerance. Moreover, it shows the local error obtained using the fixed time-step method with the same total number of attempted time-steps as the most expensive technique, the integral- (I-) controller one (being thus of similar computational cost with it). Observe that the I-controller attempts more steps than our PI-controller, for the same tolerance; therefore it is more expensive. Likewise, the I-adaptive technique has a ratio of rejected to attempted steps between 6% and 14%, while the same step rejection ratio for our PI-adaptive scheme is below 3% and is decreasing with the decrease of tolerance. It is important to remark that the existing fixed step-size method yields a local error which is up to 780 times larger than the proposed variable time-stepping method. Consequently, our variable step-size algorithm is much more accurate than the constant time-step method, for a similar computational cost.Table 3
The decay-dimerization model: the number of steps taken by the I-adaptive, PI-adaptive, and fixed stepsize methods on the time-interval[0,5].
ε
Method
Tried steps
Accepted steps
Failed steps
Error/ε
10−1
I-adapt.
1344
1150
194
0.9980
PI-adapt.
1321
1275
46
0.9952
Fixed
1344
—
—
573.9135
10−2
I-adapt.
3888
3505
383
0.9992
PI-adapt.
3852
3752
100
0.9984
Fixed
3888
—
—
685.7963
10−3
I-adapt.
11509
10831
678
0.9996
PI-adapt.
11496
11298
198
0.9994
Fixed
11509
—
—
782.6674To assess the accuracy of the proposed variable step size method, we perform simulations for 10,000 trajectories with the PI-adaptive Milstein method for the CLE with tolerances10-2 and 10-3 and with Gillespie’s algorithm. The histograms generated with these algorithms at time t=5 for each of the species are plotted in Figure 4. Again, the accuracy of our variable step size strategy is excellent.Figure 4
The decay-dimerization model: the histograms atT=5 for the species S1,S2, and S3 obtained with the SSA (green) and the PI-adaptive method for the CLE with tolerances ε=10-2 (blue) and ε=10-3 (red). The simulation uses 10,000 paths.


## 6. Conclusion

Variable step size control is critical for efficient approximation of the solution of stochastic differential equations and in particular for models which exhibit stiffness. The objective is to minimize the computational effort of the simulation, while maintaining the desired accuracy of the numerical solution. When the numerical integration enters a region of stiffness, an adaptive algorithm lowers the step size to satisfy the accuracy requirement. However, outside these regions of stiffness the algorithm relaxes the step size, thus gaining in efficiency.This paper developed a variable time-stepping strategy, based on control theory, for the strong numerical solution of a stochastic continuous model of biochemical kinetics. More precisely, we designedthe firstPI-controller for adapting the step size for a class of stochastic differential equations with noncommutative multiplicative noise, known as the chemical Langevin equation. The underlying numerical technique is a higher order of accuracy method due to Milstein. The controller, which uses low cost estimates of the local error generated by Milstein scheme, is recommended for models of biochemical systems which are mildly stiff. Unlike other variable time-stepping strategies [20] which alter the step only by halving or doubling, our technique allows for a flexible variation of the time-step. Furthermore, when a step is rejected as it does not satisfy the error criteria, the method guarantees that the statistics of the numerical solution are not biased.Numerical experiments show that the proposed PI-controller for adapting the step size issignificantly more accurate than the existing constant-step algorithms, for the same computational cost. Additionally, it has improved computational stability. Compared to the variable time-stepping method using integral control, our PI-adaptive strategy reduces greatly the number of step rejections, having a reduced computational effort for the same prescribed tolerance.Our future work will investigate how to efficiently control the time-step in the weak numerical solution of the chemical Langevin equation.



---

*Source: 567275-2015-02-09.xml*