Effective implicit finite-difference method for sensitivity analysis of stiff stochastic discrete biochemical systems.

Simulation of cellular processes is achieved through a range of mathematical modelling approaches. Deterministic differential equation models are a commonly used first strategy. However, because many biochemical processes are inherently probabilistic, stochastic models are often called for to capture the random fluctuations observed in these systems. In that context, the Chemical Master Equation (CME) is a widely used stochastic model of biochemical kinetics. Use of these models relies on estimates of kinetic parameters, which are often poorly constrained by experimental observations. Consequently, sensitivity analysis, which quantifies the dependence of systems dynamics on model parameters, is a valuable tool for model analysis and assessment. A number of approaches to sensitivity analysis of biochemical models have been developed. In this study, the authors present a novel method for estimation of sensitivity coefficients for CME models of biochemical reaction systems that span a wide range of time-scales. They make use of finite-difference approximations and adaptive implicit tau-leaping strategies to estimate sensitivities for these stiff models, resulting in significant computational efficiencies in comparison with previously published approaches of similar accuracy, as evidenced by illustrative applications.

1 Introduction

The dynamics of reaction networks in living organisms have been extensively studied in Systems Biology [1, 2]. Dynamics in such complex networks is driven by random thermal agitation. Consequently, stochastic modelling is a preferred approach for studying the behaviour of these biochemical reaction networks [3-7].

The most widely used stochastic model for describing spatially homogeneous biochemical reaction dynamics is the Chemical Master Equation (CME). State trajectories for the CME can be simulated via Gillespie's [8, 9] stochastic simulation algorithm (SSA), a Monte Carlo method. When describing systems with large molecular populations or widely varying time‐scales, Gillespie's algorithm is computationally expensive. To reduce computational costs, Gillespie [10] introduced a variant of his algorithm method called tau‐leaping, in which time steps are selected dynamically to avoid exhaustive calculations, with tolerable loss of accuracy.

The time steps for the explicit tau‐leaping strategy are limited to the fastest mode. Consequently, this strategy is not suitable for stiff biochemical systems. Specifically, the explicit tau‐leaping technique generalises the explicit Euler method for ordinary differential equations to discrete stochastic systems. Similar to Euler's method, when taking large time steps, the explicit tau‐leaping scheme shows instability for very stiff biochemical systems. Rathinam et al. [11] proposed the implicit tau‐leaping technique to overcome the stability issues of the explicit strategy, thus allowing larger time steps. The implicit tau‐leaping scheme employs larger stepsizes than the explicit tau‐leaping strategy for stiff discrete stochastic systems while producing a solution of similar accuracy as the explicit tau‐leaping scheme for the slow manifold and for the mean of the fast variables on the slow manifold [11]. The implicit tau‐leaping may damp the noise for some systems, as illustrated in [11] for systems reaching steady state.

Sensitivity analysis, which describes how model parameters related to system dynamics, is a key tool for model development and analysis [12]. Sensitivity analyses can be characterised as global (sampling over a broad region of the parameter space) or local (describing behaviour in the neighbourhood of a nominal operating point). Global approaches are computationally demanding because they require many samples from a (typically high‐dimensional) parameter space [13]. In contrast, local sensitivity analysis involves only small perturbations about a nominal parameterisation. Sensitivity analysis is a powerful tool for investigating the dynamic properties of a biochemical system. A small sensitivity reflects robustness, while large sensitivities can indicate parameters that have significant effects on outputs of interest. For biochemical reaction systems, the parameters of interest include the initial conditions, kinetic rate constants of the reactions, and environmental parameters (such as temperature). For deterministic models, local sensitivity analysis is rarely computationally demanding [14]. However, the corresponding analysis of stochastic models requires large ensembles of simulated sample paths.

Local sensitivity coefficients are usually determined as finite differences. The sensitivity of the expected value of an output f of a stochastic system to a parameter c can be described as , where represents the state of the system at time t and parameter c. The perturbation h is small compared to the nominal value of parameter c. To determine the sensitivity coefficient, two ensembles of sample paths are generated, corresponding to the parameter values c and c + h. A number of approaches to the calculation of such finite‐difference estimators have appeared recently in the literature, including the common random number (CRN) and common reaction path (CRP) approaches [15], the coupled finite‐difference (CFD) scheme [16], and the coupled tau‐leaping (CTL) algorithm [17].

For each of these sensitivity methods, sample paths of the nominal () and perturbed systems are simulated, using a common random seed to reduce variance. Anderson's CFD method, in which trajectories of the two systems are strongly coupled, was found to produce the smallest estimator variance of these approaches [18] and can be used effectively for non‐stiff models.

For stiff models, exact stochastic simulation is computationally expensive. The CTL [17] scheme was introduced to address this issue by applying an explicit tau‐leaping in the context of trajectory coupling. Because it uses the explicit tau‐leaping strategy, this method works well for moderately stiff problems.

In this paper, we propose a novel finite‐difference method for estimating sensitivities of stiff stochastic models of biochemical reaction networks. This approach, described below as the Coupled implicit‐Tau (CIT) algorithm, makes use of an implicit‐tau leaping scheme to efficiently generate sample trajectories of systems in which some reactions are in effective equilibrium. This adaptive time leaping bears similarities to an algorithm proposed by Cao et al. [19]. Similar to the CFD method mentioned earlier, an important property of the ‐leaping scheme for perturbation analysis is that the nominal and perturbed systems are strongly coupled. A similar method was developed by Anderson and Higham [20] for continuous time Markov Chains. This coupling reduces the variance in the finite‐difference estimator, allowing for a more precise measure of sensitivity.

The remaining of this paper is organised as follows. Section 2 gives an introduction to stochastic modelling and simulation of well‐stirred biochemical systems. In Section 3, we discuss the finite‐difference approximation of sensitivity for stochastic models of biochemical networks. The new CIT algorithm is presented in Section 4. Numerical tests comparing the CIT and the CFD methods are given in Section 5. Finally, we summarise our conclusions in Section 6.

2 Stochastic biochemical kinetics

We consider a well‐stirred system of biochemical reactions kept in a constant volume, at a constant temperature. The chemical species are involved in the biochemical reactions . The state vector indicates the species abundance at each time t.

Each reaction is characterised by a propensity function and a state change (stoichiometry) vector . The propensity is defined as follows: is the probability that a single reaction fire in the interval , provided that . The stoichiometry vector describes the change in molecule abundance as a consequence of reaction occurring: if the system is in the state , and then a reaction occurs, the system state becomes . The matrix is the stoichiometry matrix.

The conditional probability of the system being in state at time t, , provided that is denoted by . This probability obeys: which is known as the Chemical Master Equation.

2.1 Stochastic simulation algorithm

Gillespie proposed a Monte Carlo method, the SSA [9], for simulating sample paths with a probability in exact agreement with the solution of the CME. This algorithm is summarised below:

Initialise the system state at .

while t < T

Calculate the propensities and set .

Sample and from a uniform distribution on [0, 1], denoted by U(0, 1).

Set j, the smallest integer satisfying , be the index of the next reaction.

Set as the time for the next reaction.

Update and .

end while

2.2 Tau‐leaping methods

Exact Monte Carlo simulation algorithms [8, 9, 21] for the CMEs are often computationally expensive for problems of practical interest. An approximate strategy that reduces the computational cost of solving the CME is the tau‐leaping method, proposed by Gillespie [10]. In the tau‐leaping scheme, the system is advanced over many reactions with a preassigned stepsize . The step must satisfy the leap condition, which states that the propensities remain approximately constant over . When the leap condition is met, the number of reactions occurring between can be approximated by a Poisson random variable, , with mean and variance , when . Then the system state in may be calculated as given that . Here , with , are independent Poisson random variables. Formula (2) is called the (explicit) tau‐leaping method and was introduced by Gillespie [10].

Many biochemical systems arising in applications are stiff, displaying both slow and fast dynamics, with the fast modes being stable. The explicit tau‐leaping strategy is impractical for stiff systems, as its time step is limited to the fastest mode. To deal with this challenge, Rathinam et al. [11] proposed the implicit tau‐leaping method. The implicit tau‐leaping technique overcomes the stability issue of the explicit strategy, allowing larger steps in time. Consequently, for stiff stochastic biochemical systems, it is more efficient than the explicit method while maintaining a similar accuracy. In fact, the scheme is semi‐implicit, being implicit only in the mean part of each term , i.e. . If , the implicit tau‐leaping method updates the system state as

2.3 Stepsize selection for implicit tau‐leaping

A reversible reaction can come to equilibrium between reactants and products. When this occurs for some reversible reactions while the rest of the system is still undergoing significant variation, the system is said to be in partial equilibrium. Partial equilibrium occurs when the forward and backward propensities of a reversible reaction are approximately equal, i.e. their difference is much smaller than the propensities themselves. More precisely, if the propensities of the reversible reactions are denoted by and , the partial equilibrium condition is for some small quantity . (In the implementations below we used .)

We make use of the stepsize selection strategy introduced by Cao et al. [19]. For those reactions that are not in partial equilibrium, we demand that the mean and variance of each reactant population should satisfy where is the given tolerance. The scalar represents the highest order in which species reacts (see [22] for further details)

(A) If , take

(B) If , take , unless the left hand side of the reaction is , in which case take

(C) If , take , unless the left hand side of the reaction is , in which case take or the reaction is , in which case take

Following [22], we arrive at an efficient implementation of this leap condition by classifying all reactions that are not in partial equilibrium as critical or non‐critical, as follows. We begin by specifying the value of a control parameter, . Typically , see [22]. If some molecular amounts approach zero during the integration, then there is a trade‐off between maintaining these population numbers positive and the efficiency of the algorithm. A larger decreases the chance of negative population numbers while reducing the efficiency of the method. If a reactant is within firings of producing a zero population, it is called a critical reaction. Let us denote by , , and the set of indices of critical, non‐critical and not in partial equilibrium reactions, respectively. Let be the index set of the reaction channels that are non‐critical and not in partial equilibrium. Also, denotes the set of indices of species that are the reactants of non‐critical reactions.

Non‐critical, not in partial equilibrium reactions: The implicit tau‐leaping method is applied to the leap condition (5) implemented using the time‐step : with the auxiliary quantities Critical reactions: These reactions are implemented on the nominal or perturbed trajectory using the SSA, one critical reaction at a time. For advancing a critical reaction on the nominal trajectory two samples of the uniform random variable on the unit interval, U(0,1) are computed, and , along with the sum of all propensities of the critical reactions . Then, the index of the next critical reaction of the nominal trajectory is the smallest integer satisfying and the time to the next reaction is Similarly, for finding the next critical reaction on the perturbed trajectory, two samples of the uniform random variable on the unit interval are computed, and , and the sum of all propensities of the critical reactions . The next reaction on the perturbed trajectory has index given by the smallest integer obeying and occurs after a time step

2.4 Finite‐difference methods for sensitivity analysis of stochastic biochemical systems

Each of the techniques described above for approximating parametric sensitivities for stochastic discrete models of biochemical kinetics involves the forward finite‐difference estimator , where h represents a perturbation, c is the parameter of interest, X is the state of the chemical reaction system, and f is the output of interest. This finite‐difference estimator approximates the local sensitivity of the expected value of the quantity with respect to a parameter c, given a polynomial function f. (Note that higher‐order moments can be determined by appropriate combinations of expected sensitivities.)

Among the existing finite‐difference strategies [15, 16] for stochastic discrete models of biochemical kinetics, the CFD method due to Anderson in [16] was shown to produce the smallest estimator variance. It simulates the coupled trajectories with the next reaction method. This sensitivity estimator is based on the following tight coupling between the nominal process, , and the perturbed process, , with . Here and are independent unit rate Poisson processes.

3 CIT method

In this section, we propose a new technique for approximating local sensitivities for stochastic discrete models of biochemical systems. This method is effective and accurate for stiff models (involving multiple scales in time). Stiff systems are often encountered in applications, as biochemical systems regularly involve both fast and slow reactions. In contrast with the CRN, CRP and CFD finite‐difference schemes, which use exact SSAs to generate the nominal and perturbed trajectories, our strategy computes coupled paths using the (approximate) implicit tau‐leaping strategy. The coupling we employ is related to [23], which is used in the CFD method [16]. This coupling shares similarities to the coupling in [24] and is applied in [20] for designing multi‐level Monte Carlo methods for well‐stirred stochastic biochemical systems. The CTL method for sensitivity [17] uses finite‐differences to estimate the sensitivities and the (approximate) explicit tau‐leaping strategy to generate the coupled trajectories. However, the CTL was designed for biochemical networks that are at most moderately stiff. As opposed to these approaches, the novel CIT technique involves solving implicit equations. For stiff to very stiff models, the proposed CIT strategy allows much larger time‐steps than the previous methods. Consequently, the CIT algorithm is expected to be significantly more effective that the existing finite‐difference estimators for such systems.

In the CIT algorithm, the coupled (i.e. nominal and perturbed) implicit tau‐leaping trajectories are generated as follows with and . The Poisson random variables , and are independent. Here . The contribution of the shared term, , is expected to be significant, thus leading to a strong coupling. A consequence of this strong coupling is the tendency for reduced variance observed for this method (as illustrated in the next section). Once the Poisson terms are generated, Newton's method is applied to solve numerically both implicit equations: (14) for and (15) for , respectively.

For advancing the numerical solution, the CIT utilises an extension of the adaptive time‐stepping strategy introduced by Cao et al. [19], for the implicit tau‐leaping method, as outlined in the previous section. A candidate leap is computed for the critical and non‐critical reactions, independently, on each of the nominal and perturbed trajectories, and then the smallest leap size is chosen as the next step.

CIT algorithm

Set simulation parameters: the tolerance for tau‐leaping , the tolerance for Newton's method, TOL, the critical threshold , the final time T and the partial equilibrium parameter .

Initialise sample paths: the time and the states and .

While t < T

(a) Compute propensity functions: and .

(b) Partial equilibrium condition: for each set of reversible reactions in both systems, use (4) to check if they are in partial equilibrium.

(c) Find set of critical reactions for nominal and perturbed paths: for each non‐partial equilibrium reaction in the two systems, with or , compute where is the ‘greatest integer in’ and set .

(d) Compute candidate stepsizes, and , for non‐critical and not in partial equilibrium reactions: If no non‐critical reactions occur, . Else, determine . For on each of nominal and perturbed paths find:

(I) , the highest order at which appears in a non‐critical reaction.

(II) , the highest order at which species reacts.

(III) , using (7) and (8) and , employing (6).

(e) Compute candidate stepsizes, and , for critical reactions: calculate , , sample , from , and find , , using (10) and (12).

(f) Determine next stepsize and critical reaction index: Let and . Set for all critical reactions.

(I) If and , set .

(II) else if , a sample from . Choose as smallest integer satisfying (9). Take , .

(III) else if , a sample from . Choose as smallest integer satisfying (11). Take , .

(IV) else sample from . Choose as smallest integer satisfying (9). Take , .

(g) Step over non‐critical reactions: For each ,

(I) Generate samples from Poisson distributions

(II) Apply Newton's method to solve each of the systems

(III) Update , .

(g) Implement the step: update time and system states

(h) Approximate sensitivity on the sample path: at the current time.

4 Numerical results

In this section, we compare the CIT method with the CFD strategy on some examples of stiff biochemical systems. Recall that, of the published finite‐difference techniques for estimating the sensitivities, the CFD technique provides estimates with the lowest variance [16].

In our comparisons, we use ensembles of 10,000 paths of the CFD and of the new CIT methods, respectively. We apply the CIT algorithm as described above with tolerance , TOL = 0.01, and . We show that the CIT method produces smaller variances than the CFD strategy for the first two models and similar variances for the third model. The CIT estimator is found to be significantly faster than the CFD. The efficiency is measured by where the CPU time to simulate 10,000 trajectories is considered in each case.

4.1 Decay‐dimerisation model

The decay‐dimerisation model of [11] consists of three molecular species involved in four chemical reactions (Fig. 1). The reactions and propensities are given in Table 1, along with a set of nominal values for the rate constants.

Fig. 1

Decay‐dimerisation model reaction chain

Table 1

Decay‐dimerisation model

	Reaction	Propensity	Nominal rate constant
R_1	S_1⟶C_1⊘	a_1=C_1X_1	C_1=0.05
R_2	S_1+S_1⟶C_2S_2	a_2=C_2X_1(X_1−1)/2	C_2=50
R_3	S_2⟶C_3S_1+S_1	a_3=C_3X_2	C_3=1,000,000
R_4	S_2⟶C_4S_3	a_4=C_4X_2	C_4=0.05

The system was simulated on the time‐interval [0, 1], with initial conditions and the parameter . The mean of the state variable (i.e. the number of molecules), for the adaptive implicit tau‐leaping algorithm and for the next reaction method (used for the CFD), are plotted in Fig. 2. Fig. 2 shows the standard deviation of this state variable. The estimated sensitivity of with respect to the parameter and its standard deviation are shown in Figs. 2 and d. The perturbation parameter is h = 0.05 (i.e. 0.1% of the nominal parameter value). Fig. 2 demonstrates that the variance of the CIT estimator is small compared to that of the CFD, demonstrating accuracy. Moreover, the speed‐up of the CIT scheme over the CFD technique for estimating sensitivities for this particular simulation is

Fig. 2

Decay‐dimerisation model: 10,000 trajectories were generated on the time‐interval [0,1], with initial condition and parameters in Table 1

, Mean and standard deviation of the number of molecules for species were calculated by the next reaction method and the adaptive implicit tau‐leaping algorithm, , Finite‐difference estimates of the sensitivity of the abundance of with respect to , and the standard deviation of the estimators, for the CFD and CIT

4.2 Genetic positive feedback loop

We next consider a simple model of positive feedback in gene expression (Fig. 3), as presented in [25]. Referring to Table 2, x represents a monomeric protein, y the protein dimer, the unoccupied regulatory site on the gene coding for x, the dimer‐occupied site, and m, the mRNA transcript. The reactions, propensities and a set of nominal parameter values are included in the table.

Fig. 3

Schematic diagram of genetic positive feedback loop model

Table 2

Genetic positive feedback loop model

	Reaction	Propensity	Nominal rate constant
R_1	x+x⟶C_1y	a_1=C_1X(X−1)/2	C_1=5000
R_2	y⟶C_2x+x	a_2=C_2Y	C_2=10^6
R_3	y+d_0⟶C_3d_r	a_3=C_3YD_0	C_3=5000
R_4	d_r⟶C_4y+d_0	a_4=C_4D_r	C_4=10^6
R_5	d_0⟶C_5d_0+m	a_5=C_5d_0	C_5=10
R_6	d_r⟶C_6d_r+m	a_6=C_6D_r	C_6=20
R_7	m⟶C_7m+x	a_7=C_7M	C_7=1
R_8	x⟶C_8⊘	a_8=C_8X	C_8=0.8
R_9	m⟶C_9⊘	a_9=C_9M	C_9=7

We ran simulations from initial molecular amounts of over the time‐interval [0,2], with .

Fig. 4 presents the evolution of the mean amount of the x molecules over 10,000 paths, generated with the CIT algorithm and the next reaction method, respectively. The standard deviation of the molecular count of x as a function of time, for each of the two algorithms, is shown in Fig. 4. The behaviours of the estimated sensitivity of the x molecular numbers with respect to the parameter , using the CIT and the CFD methods are presented in Fig. 4, whereas the corresponding standard deviations of the CIT and CFD estimators are given in Fig. 4. The simulations are performed with a perturbation h = 0.5 (i.e. 0.01% of the nominal parameter value). From Fig. 4, we observe that the CIT estimator variance is low compared to the variance of the CFD estimator, therefore the sensitivity estimation of the new CIT method is more accurate. This low variance is confirmed as shown in Fig. 4. For the set of parameters used, the speed‐up, on time interval [0, 2], of the CIT over the CFD is significant

Fig. 4

Genetic positive feedback loop model. 10,000 sample paths with initial condition and parameters as given in Table 2 were generated on the time‐interval [0,2]

, Mean and standard deviation of the number of molecules for species x were calculated by the next reaction method and the adaptive Implicit tau‐leaping algorithm, , Mean and standard deviation of the finite‐difference estimators determined via the CFD and implicit tau‐leaping methods, of the sensitivity of the abundance of x to the parameter

4.3 Collins toggle switch model

The Collins toggle switch [26] is a gene regulatory network that exhibits bistability: it consists of two genes, each encoding a repressor of the other. Referring to Fig. 5, the species and are the gene's protein products, while and denote the corresponding mRNA transcripts. The parameters and denote the maximal transcription rates. Furthermore, and are the degrees of non‐linearity in the repression mechanisms. The system exhibits bistable‐like behaviour when and the maximal expression rates are adequately large. We introduce a scaling factor k on the propensity of mRNA transcription and degradation to allow tuning of the model stiffness. For Figs. 6 and 7, we used k = 1000. Table 3 lists the reactions, their propensities and nominal rate constants. This illustrative parameterisation was chosen to achieve regular transitions between well‐separated quasi‐steady states and significant stiffness.

Fig. 5

Collins toggle switch model reaction scheme diagram

Fig. 6

Collin's toggle switch model: a sample path of species and , with initial condition and the parameters in Table 3 generated with the implicit tau‐leaping method on the time‐interval [0, 8000]

Fig. 7

Collins toggle switch model. 10,000 sample paths with initial condition and parameters as given in Table 3 were generated on the time‐interval [0, 2000]

, Mean and standard deviation of the number of molecules for species were calculated by the next reaction method and the adaptive Implicit tau‐leaping algorithm, , Mean and standard deviation of the finite‐difference estimators determined via the CFD and implicit tau‐leaping methods, of the sensitivity of the abundance of to the parameter

Table 3

Collins toggle switch model

	Reaction	Propensity	Nominal rate constant
R_1	⊘⟶C_1m_1	a_1=kα_11+(X_2)^β	α_1=C_1=28.98β=4
R_2	m_1⟶C_2⊘	a_2=kC_2X_3	C_2=0.23
R_3	m_1⟶C_3p_1+m_1	a_3=C_3X_3	C_3=0.23
R_4	p_1⟶C_4⊘	a_4=C_4X_1	C_4=0.23
R_5	⊘⟶C_5m_2	a_5=kα_21+(X_1)^γ	α_2=C_5=28.98γ=4
R_6	m_2⟶C_6⊘	a_6=kC_6X_4	C_6=0.23
R_7	m_2⟶C_7p_2+m_2	a_7=C_7X_4	C_7=0.23
R_8	p_2⟶C_8⊘	a_8=C_8X_2	C_8=0.23

This system was integrated on the time‐interval [0, 2000], with initial conditions and . Sample trajectories for the species and , simulated with the underlying implicit tau‐leaping method are shown in Figs. 6 and b. Bistable behaviour is also observed for the species and (not shown). The mean and standard deviation number of the number of molecules for the implicit tau‐leaping algorithm and for the next reaction method are plotted in Figs. 7 and b, respectively. We remark that, for this bistable model, the variance of the implicit tau‐leaping technique matches well the variance of the (exact) next reaction method, in contrast with the behaviour of the implicit tau‐leaping scheme observed by Rathinam et al. [11] for systems reaching a steady state, where the implicit scheme reduced the variance of the numerical solution. Figs. 7 and d present the finite‐difference estimation of the sensitivity of the molecule count with respect to the parameter and the estimator's standard deviation for each of the CIT and CFD algorithms. In these simulations, the perturbation parameter is h = 0.05 (i.e. 0.2% of the nominal parameter value). The estimation of the sensitivity is similar for the CIT and the CFD methods, while the standard deviation of the CIT estimator is slightly larger than that of the CFD estimator. However, for the set of parameters in Table 3, the speed‐up of the CIT over the CFD is 74‐fold.

To gauge the performance of the method as stiffness increases, the performance of the CIT and CFD methods was studied for various values of scaling parameter k. The chosen range produced stiff to very stiff model formulations, resulting in speed‐up of the new CIT strategy compared to the existing CFD method of up to 468 times, as reported in Table 4. For non‐stiff models, the CIT algorithm will perform no better than the CFD method. For this model, a similar computational time for the two algorithms is obtained when the propensities of the fastest and slowest reactions are separated by about two orders of magnitude.

Table 4

Collins toggle switch model: the speed‐up of the CIT compared to the CFD for estimating the sensitivity of with respect to for h = 0.05 on time interval [0, 2000]

Method	Stiffness parameter k	Speed‐up
CIT	300	10.16
CIT	1000	74.71
CIT	3000	468.43
CFD	—	1

As shown in panel (d) of Figs. 1, 4 and 7, the variance of the CIT estimator is not always comparable to that of the CFD estimator (smaller in the first two examples, larger in the third). For first two models, we observe that our CIT method is more accurate and far more efficient than the existing CFD strategy.

For the third model, when the value of stiffness parameter k grows, our CIT method becomes increasingly more efficient than the CFD scheme. On the other hand, for the Collins toggle switch model, the variance of the CIT estimator is slightly larger than that of the CFD. The implicit tau‐leaping scheme damps the noise for systems reaching a steady state [11]. However, for the toggle switch model, the implicit tau‐leaping scheme does not cause noise reduction. Trajectories frequently switch between two states, the model exhibiting bistable behaviour. This behaviour restricts the noise damping property of the implicit tau‐leaping scheme and leads to a slightly larger variance of the CIT algorithm than that of the CFD, unlike for the previous two models. According to our numerical experiments, we conclude that our CIT method is expected to be more accurate and significantly more efficient than the CFD technique when the stiff system reaches a steady state. A theoretical study of the properties of the finite‐difference CIT sensitivity estimator and appropriate values of the perturbation parameter h leading to optimal convergence rates will be considered in our future work.

5 Conclusion

Noise plays an important role in the behaviour of well‐stirred biochemical systems when some species have low molecular counts. Usually, such systems are represented using the discrete stochastic CME model. For the class of biochemical systems that are mathematically stiff, implicit ‐leaping schemes are preferred over exact Monte Carlo SSAs; implicit methods are considerably more efficient for accurately determining the slow stochastic variables of the system, and for determining the mean behaviour of the fast variables. Indeed, the implicit tau‐leaping strategy allows large stepsizes while maintaining the solution close to the slow manifold. By contrast, exact stochastic simulation techniques are forced to take very small time‐steps when stiffness is encountered.

This work developed a new finite‐difference estimator of sensitivities for stochastic discrete models of biochemical kinetics. The new sensitivity estimator, the CIT or CIT, employs an adaptive implicit ‐leaping method to simulate the nominal and perturbed paths. To enhance the accuracy of the estimation, the CIT applies a strong coupling between the nominal and perturbed trajectories. The new sensitivity estimator is superior to Anderson's CFD estimator [16] for models of biochemical systems that are stiff to very stiff. We demonstrated that the CIT has a considerably reduced computational cost compared to the CFD while retaining an equivalent accuracy. Thus, for the class of stiff to very stiff biochemical reaction networks, the coupled implicit ‐leaping method is a preferred choice for sensitivity analysis.