Accelerated Sensitivity Analysis in High-Dimensional Stochastic Reaction Networks

Existing sensitivity analysis approaches are not able to handle efficiently stochastic reaction networks with a large number of parameters and species, which are typical in the modeling and simulation of complex biochemical phenomena. In this paper, a two-step strategy for parametric sensitivity analysis for such systems is proposed, exploiting advantages and synergies between two recently proposed sensitivity analysis methodologies for stochastic dynamics. The first method performs sensitivity analysis of the stochastic dynamics by means of the Fisher Information Matrix on the underlying distribution of the trajectories; the second method is a reduced-variance, finite-difference, gradient-type sensitivity approach relying on stochastic coupling techniques for variance reduction. Here we demonstrate that these two methods can be combined and deployed together by means of a new sensitivity bound which incorporates the variance of the quantity of interest as well as the Fisher Information Matrix estimated from the first method. The first step of the proposed strategy labels sensitivities using the bound and screens out the insensitive parameters in a controlled manner. In the second step of the proposed strategy, a finite-difference method is applied only for the sensitivity estimation of the (potentially) sensitive parameters that have not been screened out in the first step. Results on an epidermal growth factor network with fifty parameters and on a protein homeostasis with eighty parameters demonstrate that the proposed strategy is able to quickly discover and discard the insensitive parameters and in the remaining potentially sensitive parameters it accurately estimates the sensitivities. The new sensitivity strategy can be several times faster than current state-of-the-art approaches that test all parameters, especially in “sloppy” systems. In particular, the computational acceleration is quantified by the ratio between the total number of parameters over the number of the sensitive parameters.


Variance of the estimator of the sensitivity bound
In order to get confidence intervals for an estimator that consists of simpler estimators, e.g., the estimator of the sensitivity bound ((4) in Main Text), the variance of the complex estimator must be computed. The Delta method gives a way to express the variance of the complex estimator using the variances of the simpler ones. In the following discussion we denote by N (0, Σ) the normal distribution with mean 0 and covariance matrix Σ and by d − → the convergence in distribution.
Theorem 1.1 (Delta method). Suppose {X n } ∞ n=0 is a sequence of random vectors in R K such that √ n (X n − θ) for some θ ∈ R K and Σ ∈ R K×K . Let g : R K → R and J(θ) = ∇g(x)| x=θ . Then, Proof. See Theorem 3.1 in [1].
The practical interpretation of the Delta method is that for large fixed N the variance of g(X N ) can be approximated by 1 N J(θ) T ΣJ(θ). After a Definition and an intermediate Lemma we will prove a Lemma that gives and approximation for the variance of the sensitivity bound ((4) in Main Text).
Definition 1.1 (moments). For a random variable X we define the moments, and the central moments Moreover, for two random variables, X and Y , we define the joint moments by, are independent and identically distributed random variables. Then the following are true: and holds that 2. an unbiased estimator for the variance, µ 2 , of distribution f is given by with variance equal to Proof. The proof of the first part is an application of the Central Limit Theorem and can be found in [1]. The second part can be proved using the Delta method and the details can be found in Section 3.1 of [1].
Let X be an estimator for Proof. For simplicity in the presentation we will work with the estimator of (8). The same result holds for the estimator (8), see the Example 3.2 in [1].
First, notice that X S 2 Y is equal to φ(X, Y , Y 2 ) for ϕ(x, y, z) = x(y 2 − z). Moreover, from the multivariate Central Limit Theorem, see [1], it holds that, where the covariance matrix Σ is equal to The Jacobian of ϕ evaluated at θ = (m 1 , n 1 , n 2 ) is equal to Finally, use Theorem 1.1 to conclude equation (10) with λ given by Equation (11) follows after basic algebraic manipulations.  (22). This result indicates that the variance of the bound estimator is much less than that of the coupling estimator. The ratio of variance for Parameter 11 does not appear in the graph because the variance of the coupling estimator is equal to zero.

Technical details regarding the computational cost
In this section, we compare the computational cost of (a) the proposed strategy, involving estimation of the upper bound in the inequality ((4) in Main Text) and the subsequent estimation of just the most significant sensitivity indices (SIs),S k, m (see (12) in Main Text), to (b) the estimation of the full sensitivity matrix, using the finite-difference coupling method for each one of the entries of the sensitivity matrix. The comparison will be done under the requirement that the relative confidence intervals of all estimators involved in both approaches will be of length equal to δ 1. The relative confidence interval is defined as the usual confidence interval, however the estimator of the standard deviation is replaced by a normalized standard deviationσ whereσ is the estimated standard deviation of the random variable X and x = 1 x i is the estimator of the mean of X. With this definition, random variables with large mean values have standard deviation comparable with random variables with small mean values. plays the role of a separator between large and small values. In our study we used = 0.01. Next, we set notation before proceeding to the comparison of the estimators. The time averaged observable will be that of equation (1) where V c k, is the estimated normalized standard deviation, as defined in (16), using M k, samples, and α = 1.96 is the inverse of the cumulative normal distribution evaluated at 0.975 leading to a 95% confidence interval, see [2]. The superscript c in (17) stands for the coupling method, see Section "Step 2: Finding the most sensitive parameters" in Main Text. Under the requirement that L c k, = δ we get the minimum number of samples for every pair (k, ) where x is the ceiling function defined by x = min{n ∈ Z : n ≥ x}. In the coupling method, as a variant of the SSA method [3], for every perturbed parameter we compute the SIs of all species. Hence, the minimum number of trajectories needed to be produced by this method is determined by the species with the highest variance. Thus, we define the quantity Let M c be the sum of M c k over all k, i.e. the total number of samples required by the coupling method such that all estimated SIs will have a relative confidence interval equal to δ. This number will be compared to the minimum samples required by the method proposed in this paper under the requirement that all estimators' relative confidence interval will be of the same length.
Working in the same way, we find the minimum number of trajectories required in order the relative confidence interval of the estimator of the sensitivity bound (SB) ((4) in Main Text) to be less or equal to δ, where the superscript b stands for "bound". The statistical estimators of the SB are presented in supporting information File S2 and the variance of the estimator is computed using Lemma 1.2. In one trajectory we collect samples from both FIM and IAT, i.e. the minimum number of iterations will be determined by the bound of the sensitivity of the -th species of the k-th parameter with the higher variance. Thus we define the quantity which is the minimum number of samples required for the upper bound estimator to have a relative confidence interval of less or equal to δ. In order to compare the straightforward SI estimator ((12) in Main Text) with the strategy proposed here, we have to compare the numbers in (19) and (21). Note that in our strategy, after the first screening using the SB is performed, we have to estimate accurately the SIs using the coupling method. Thus the ratio of the total computational cost of the proposed strategy to the computational cost of coupling in terms of number of samples is, where I and S are the sets of indices of insensitive and sensitive parameters, respectively. Note that in the above equation we dropped the ceiling functions in (19) and (21) making the calculation easier and introducing a negligible error. The interpretation of G 1 is "how many times is the coupling faster than Step 1 of the proposed strategy", i.e. the computation of the upper bound, and that of G 2 "how many times is the coupling applied to all sensitivity indices is faster than Step 2 of the proposed strategy", i.e. neglecting SIs with small values. The ratio 1/G is the total speed-up of the proposed two-step strategy.