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.
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 is the 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. Section 2 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
Assume biochemical species participate in reaction channels . The biochemical system, held at constant temperature, is homogeneous. The system state at time is denoted by , where represents the number of molecules at time , for any . Mathematically, the system state is modelled as a Markov process. When one reaction fires, the state of the system is updated using a vector called the state-change vector corresponding to the reaction . This is an -dimensional vector with entries , denoting the variation in the molecular species caused by one reaction . Then, is the stoichiometric matrix of the biochemical system. To each reaction it corresponds a propensity defined as follows: for an infinitesimal time increment , is the probability that, given the state at time , one reaction happens during .
For a unimolecular reaction the propensity function is , while for a bimolecular reaction the propensity is of the form when and when .
The system state may be represented as where are independent unit rate Poisson processes [24].
Assume that the leap condition is satisfied: there exists a time step small enough such that for each . Under this assumption, the representation (3) may be approximated by Equation (5), due to Gillespie, is known as the tau-leaping formula [8]. A further reduction is possible if, in addition, may be chosen large enough such that each reaction fires many times in the interval . More precisely, exists such that for any .
A step size 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 where are normal random variables with mean and variance . Regarding the step-size as an infinitesimal in (7) leads to the following equation: where 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 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 an -dimensional Wiener process, , in the form In the differential equation (9), denotes an -dimensional stochastic process and and are -dimensional drift and diffusion coefficients, respectively. The initial condition is for .
The SDE is assumed to have noncommutative noise, as is the case of the chemical Langevin equation model. If the differential operator is defined as then the SDE is called noncommutative if for some , with , It is called a commutative 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 of the exact solution of (9), after steps with step-size , is denoted by . This approximation is said to have strong order of convergence if there exists a constant , independent of and , such that the following inequality holds for any : where is some norm of a vector of dimension and is the expectation of a random variable.
Numerical methods of strong order of accuracy for a noncommutative SDE require, on each subinterval on each interval , the simulation of the Wiener increments and of either the double Itô integrals [19], or the Levy areas.
Approximations of the double Itô integrals for , using the truncation after terms of their Karhunen-Loève or Fourier series expansion [14, p. 198–203] (see also [19]), are computed as For , and , are independent normally distributed random variables with mean 0 and variance 1. Here the following notation is used: An accurate simulation of the double Itô integrals requires a minimum value (see also [19]).
Now, let us introduce a numerical method of strong order of accuracy , typically used in the literature to simulate the CLE, namely, the Milstein scheme [14]. The Milstein method on the time interval computes where is defined in (11) and the Wiener increments are .
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 and the diffusion coefficients for .
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 , becomes for any .
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. Let be the estimation of the local error for the time interval , 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 for a drift 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 with and being the Jacobian of . Here denotes a column vector with the th entry being .
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,
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 is tried and rejected on the current interval, , its corresponding Wiener increments, were sampled. Denote the sampled value . Then, a smaller time-step is tried, , and the corresponding Wiener increments on the subintervals and are created To ensure that the integration preserves the already created Brownian path, a Brownian bridge is employed: Here denotes a new normally distributed random variable with mean and variance .
One step-size change method is the integral- (I-) controller for predicting the future step, , based on the current step and the current local error . The I-controller is for some constant . The safety factor is introduced to reduce the number of rejected steps. The values considered in our simulations are and . 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].
The PI-controller is given by For a further reduction of number of step rejections, the adaptive algorithm employs a more conservative approach to controlling the time-step: which is similar to an adaptive procedure for the numerical solution of ODE [25]. Note that the factors and 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 Here is the order of the local error, , where 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 , referred to as the integral gain, and , 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: . 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: and .
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 The propensities characterizing these reactions are The stoichiometric matrix has as columns the state-change vectors of the reactions above. The system parameters are as follows: , , , , , and for the reaction rate constants, and for the initial conditions. The interval of integration is .
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 tolerances , and . 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 to and it decreases as the tolerance decreases. For the I-adaptive technique the same ratio of rejected to attempted steps is much larger, between and for the tolerances tried. The fixed step size method produces a local error of up to times larger than the adaptive methods, for the same computational cost.
To verify the accuracy of the proposed PI-controller for adapting the time-step when simulating the CLE, we compare the histograms at time 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.
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: The reaction rate constants take the values , , , and . The state change vectors of the reactions in the Schlögl model are the corresponding columns of the stoichiometric matrix: The reaction propensities may be computed as The numbers of molecules of the species and have constant values, and . The initial condition for the molecular number of the species is .
The system is integrated on the time-interval . 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 , , and . 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 and , while for the I-adaptive scheme ranges from to .
In Figure 2 we show the evolution in time of the error scaled by the tolerance, on a sample path, for the tolerance 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 for the species computed with the PI-adaptive method for the CLE with and and with the SSA. The histograms agree very accurately.
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: The propensity functions corresponding to these reactions are with the reaction rate parameter values , , , and . The solution of this biochemical reaction model is subject to the initial conditions and . The integration is performed on the time-interval . In addition, the stoichiometric matrix for this system is
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. Table 3 reports the numerical results: the number of attempted, accepted, and rejected steps of each of the variable step-size strategies, for tolerances , and , 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 and , while the same step rejection ratio for our PI-adaptive scheme is below 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 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.
To 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 tolerances and and with Gillespie’s algorithm. The histograms generated with these algorithms at time for each of the species are plotted in Figure 4. Again, the accuracy of our variable step size strategy is excellent.
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 designed the first PI-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 is significantly 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.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
This work was supported in part by a grant from the Natural Sciences and Engineering Research Council of Canada (NSERC).