Abstract
Background
Stochastic modeling and simulation provide powerful predictive methods for the intrinsic understanding of fundamental mechanisms in complex biochemical networks. Typically, such mathematical models involve networks of coupled jump stochastic processes with a large number of parameters that need to be suitably calibrated against experimental data. In this direction, the parameter sensitivity analysis of reaction networks is an essential mathematical and computational tool, yielding information regarding the robustness and the identifiability of model parameters. However, existing sensitivity analysis approaches such as variants of the finite difference method can have an overwhelming computational cost in models with a highdimensional parameter space.
Results
We develop a sensitivity analysis methodology suitable for complex stochastic reaction networks with a large number of parameters. The proposed approach is based on Information Theory methods and relies on the quantification of information loss due to parameter perturbations between timeseries distributions. For this reason, we need to work on pathspace, i.e., the set consisting of all stochastic trajectories, hence the proposed approach is referred to as “pathwise”. The pathwise sensitivity analysis method is realized by employing the rigorouslyderived Relative Entropy Rate, which is directly computable from the propensity functions. A key aspect of the method is that an associated pathwise Fisher Information Matrix (FIM) is defined, which in turn constitutes a gradientfree approach to quantifying parameter sensitivities. The structure of the FIM turns out to be blockdiagonal, revealing hidden parameter dependencies and sensitivities in reaction networks.
Conclusions
As a gradientfree method, the proposed sensitivity analysis provides a significant advantage when dealing with complex stochastic systems with a large number of parameters. In addition, the knowledge of the structure of the FIM can allow to efficiently address questions on parameter identifiability, estimation and robustness. The proposed method is tested and validated on three biochemical systems, namely: (a) a protein production/degradation model where explicit solutions are available, permitting a careful assessment of the method, (b) the p53 reaction network where quasisteady stochastic oscillations of the concentrations are observed, and for which continuum approximations (e.g. mean field, stochastic Langevin, etc.) break down due to persistent oscillations between high and low populations, and (c) an Epidermal Growth Factor Receptor model which is an example of a highdimensional stochastic reaction network with more than 200 reactions and a corresponding number of parameters.
Keywords:
Biochemical reaction networks; Sensitivity analysis; Relative entropy rate; Pathwise Fisher information matrix; p53 model; EGFR modelBackground
The need of an intrinsic understanding of the interplay between complexity and robustness of biological processes and their corresponding design principles is welldocumented, see for instance [15]. The concept of robustness can be described as “a property that allows a system to maintain its functions against internal and external perturbations” [3]. When referring to mathematical models of complex biological processes, one of the mathematical tools to describe the robustness of a system to perturbations is sensitivity analysis which attempts to determine which parameter directions (or their combinations) are the most/least sensitive to perturbations and uncertainty, or to errors resulting from experimental parameter estimation. Recently there has been significant progress in developing sensitivity analysis tools for lowdimensional stochastic processes, modeling wellmixed chemical reactions and biological networks. Some of the mathematical tools included loglikelihood methods and Girsanov transformations [68], polynomial chaos [9], finite difference methods and their variants [10,11] and pathwise sensitivity methods [12]. However, existing sensitivity analysis approaches can have an overwhelming computational cost, either due to high variance in the gradient estimators, or in models with a highdimensional parameter space, [13].
The aforementioned methods focus on the sensitivity of stochastic trajectories and corresponding averages. However, it is often the case that we are interested in the sensitivity of probability density functions (PDF), which are typically nonGaussian in nonlinear and/or discrete systems. In that latter direction, there is a broad recent literature relying on information theory tools, and where sensitivity is estimated by using the Relative Entropy and the Fisher Information Matrix between PDFs, providing a quantification of information loss along different parameter perturbations. We refer to [1418] for the case when the parametric PDF is explicitly known. For instance, in [16] the parametric PDF’s structure is known as it is obtained through an entropy maximization subject to constraints. Knowing the form of the PDF allows to carry out calculations such as estimating the relative entropy and identifying the most sensitive parameter combinations. Furthermore, the pathwise PDFs are also known in reaction networks when a Linear Noise Approximation (LNA) is employed and for this case the relative entropy can be explicitly computed allowing thus to carry out parametric sensitivity analysis, [18]. However, for complex stochastic dynamics of large reaction networks, spatial Kinetic Monte Carlo algorithms and molecular dynamics, such explicit formulas for the PDFs are not available in general.
In [19], we address such challenges by introducing a new methodology for complex stochastic dynamics based on the Relative Entropy Rate (RER) which provides a measure of the sensitivity of the entire timeseries distribution. Typically, the space of all such timeseries is referred in probability theory as the “path space”. RER measures the loss of information per unit time in path space after an arbitrary perturbation of parameter combinations. RER and the corresponding Fisher Information Matrix (FIM) become computationally feasible as they admit explicit formulas which depend only on the propensity functions (see (4) and (6), respectively). In fact, we showed in [19] that the proposed pathwise approach to sensitivity analysis has the following features: First, it is rigorously valid for the sensitivity of longtime, stationary dynamics in path space, including for example bistable, periodic and pulselike dynamics. Second, it is a gradientfree sensitivity analysis method suitable for highdimensional parameter spaces as the ones typically arising in complex biochemical networks. Third, the RER method does not require the explicit knowledge of the equilibrium PDFs, relying only on information for local dynamics and thus making it suitable for nonequilibrium steady state systems. In [19], we demonstrated these features by focusing on two classes of problems: Langevin particle systems with either reversible (gradient) or nonreversible (nongradient) forcing, highlighting the ability of the method to carry out sensitivity analysis in nonequilibrium systems; and spatially extended Kinetic Monte Carlo models, showing that the method can handle highdimensional problems.
In this paper, we extend and apply the pathwise sensitivity analysis method in [19] to biochemical reaction networks, and demonstrate the intrinsic sensitivity structure of the network. Such systems are typically modeled as jump Markov processes and they are simulated using either exact algorithms such as the Stochastic Simulation Algorithm (SSA), [2022] and the nextreaction method [23], or by employing approximations such as mean field ODEs, tauleap [24] and stochastic Langevin methods [25].
We show that the proposed pathwise method allows us to discover the intrinsic sensitivities of the reaction network by decomposing the FIM into diagonal blocks. The blockdiagonal structure of the proposed FIM reveals, in a straightforward way, the sensitivity interdependencies between the system parameters. For instance, if each propensity function depends only on one parameter –usually the reaction constant– then the FIM is a diagonal matrix (see (14)). The sparse representation of the FIM can be essential in optimal experimental design as well as in parameter identifiability and robustness where each subset of the parameters defined by a block of the FIM can be treated separately. Moreover, our earlier rigorous analysis [19] for the stationary regime suggests suitable extensions in the transient case which are here tested and validated. Finally, we present strategies for efficiently and reliably implementing the proposed method for highdimensional, complex stochastic systems using an array of existing accelerated versions of the SSA algorithm such as mean field, stochastic Langevin, τleap approximations and their variants, [21,2427].
We test the proposed set of methods and computational strategies in three examples of biochemical networks. First, we consider a prototypical protein production/degradation model, i.e, a singlespecies birth/death model, with explicitly known formulas for the stationary and the timedependent distribution. This model serves as a benchmark where the differences between the proposed pathwise FIM and the stationary FIM are highlighted. Second, we study the parameter sensitivities of a p53 gene model for cell cycle regulation and response to DNA damage, that incorporates the feedback between the tumor suppressor p53 gene and the oncogene Mdm2 [28]. This is a reaction network that exhibits random oscillations in its steady state, and for which continuum approximations of the SSA such as LNA break down due to persistent oscillations between high and low populations. Using the proposed method, we also study a far more complex network, the epidermal growth factor receptor (EGFR) model, describing signaling processes between mammalian cells [2931]. This is a highdimensional system both in the number of variables and parameters, including 94 species and 207 reactions. Having a gradientfree method for this example with parameter space of dimension 207 provides a significant advantage over gradient methods such as finite differencing, where the computation of a very high number of partial derivatives and/or directional derivatives is needed and with possibly significant variance that scales with the dimension, [11]. By contrast, the eigenvalue/eigenvector analysis of the proposed FIM identifies the order from least to most sensitive directions (determined by the eigenvectors of the FIM) by the corresponding eigenvalues.
In Methods, we present the derivation of the Relative Entropy Rate and its corresponding Fisher Information Matrix for continuoustime jump Markov processes as well as we reveal the blockdiagonal structure of the FIM for commonly encountered reaction networks, continued by the presentation of both unbiased and biased –but accelerated– statistical estimators for RER and FIM. Then, in the Results, we apply and validate the proposed pathwise sensitivity analysis methodology in three complex biological reaction networks.
Methods
We consider a wellmixed reaction network with N species, S = {S_{1}, …, S_{N}}, and M reactions, R = {R_{1}, …, R_{M}}. The state of the system at any time t ≥ 0 is denoted by an Ndimensional vector X(t) = [X_{1}(t), …, X_{N}(t)]^{T} where X_{i}(t) is the number of molecules of species S_{i} at time t. Let the Ndimensional vector ν_{j} correspond to the stoichiometry vector of jth reaction such that ν_{i,j} is the stoichiometric coefficient of species S_{i} in reaction R_{j}. Given that the reaction network at time t is in state X(t) = x, the propensity function, a_{j}(x), is defined so that the infinitesimal quantity a_{j}(x)dt gives the transition probability of the jth reaction to occur in the time interval [t, t + dt]. Propensities are typically dependent on the state of the system and the reaction conditions (i.e., external parameters) of the network such as temperature, pressure, etc. Mathematically, is a continuoustime, timehomogeneous, jump Markov process with countable state space . The transition rates of the Markov process are the propensity functions a_{j}(·), j = 1, …, M. The transition rates determine the clock of the updates (or jumps) from a current state x to a new (random) state x^{′} through the total rate while the transition probabilities of the process are determined by the ratio . We refer to Algorithm 1 for the details of the stochastic simulation.
Relative entropy
Assume that two probability distributions (or more generally probability measures) and have corresponding probability densities p = p(x) and . Then, the Relative Entropy or KullbackLeibler divergence of with respect to is defined as [32,33]
In a more general setting, relative entropy is defined as where is a function known as RadonNikodym derivative while the integration is performed with respect to the probability measure , [34]. A necessary condition for the relative entropy to be welldefined is that the RadonNikodym derivative exists which is satisfied when is absolutely continuous with respect to . Relative entropy has been utilized in a diverse range of scientific fields from statistical mechanics [34] to coding in telecommunications (information theory) [33] and finance [35], and it possesses the following three fundamental properties:
(i) it is always nonnegative,
(ii) it equals to zero if and only if almost everywhere, and,
(iii)
if and only if and are absolutely continuous with respect to each other.
From an information theory perspective, relative entropy quantifies the loss of information when is utilized instead of , [33]. In other words, relative entropy quantifies the inefficiency of assuming an incorrect or perturbed distribution instead of employing the true distribution . Therefore, even though not a metric, relative entropy has been used as a suitable quantity for the assessment of parametric sensitivities since the higher the relative entropy (i.e., the information loss) in some perturbed direction, the larger the sensitivity should be in this direction.
Pathwise relative entropy and relative entropy rate
Proceeding to the pathwise formulation of the relative entropy, we assume that the propensities depend on a parameter vector (i.e., ) while the continuoustime jump Markov process lies in the stationary regime. We denote by μ^{θ}(x) the steady state (or stationary) distribution of the stochastic process X(t). The stationary path distribution of the process in the interval [0, T] is denoted by . Notice that path distributions (i.e., timeseries distributions) are highdimensional complex objects; for instance, if we consider the simpler discretetime Markov chain, , defined by the transition probability density p(z, z^{′}), then, utilizing repeatedly the Markov property, the stationary path distribution of the timeseries (z_{0}, z_{1}, …, z_{T}) is given by
Proceeding, we consider another continuoustime jump Markov process defined by perturbing the propensity functions by a small vector . The corresponding steady state and path distributions of are denoted by μ^{θ+ε}(x) and , respectively. Let the two path distributions and be absolutely continuous with respect to each other which is satisfied when if and only if holds for all x ∈ E and j = 1, …, M. Then, the Relative Entropy of the path distribution with respect to is defined similarly to (1) as
where is the RadonNikodym derivative of with respect to . In fact, using the Girsanov’s formula, we can obtain an explicit expression for the RadonNikodym derivative in terms of the propensities, [34]. In the context of sensitivity analysis, the pathwise relative entropy is a measure of information loss due to an εperturbation of the model parameters, and consequently it is a natural measure of parametric sensitivity.
Moreover, in the stationary regime, relative entropy increases linearly in time, hence the Relative Entropy Rate (RER) which is the time average of the relative entropy,
is a welldefined quantity, [36]. As first proposed in [19], is a suitable timeindependent measure of sensitivity: it measures the rate of the loss of information due to an εperturbation of the model parameters, in the longtime, stationary dynamics regime of the stochastic process. Furthermore, RER admits an explicit formula given by (see Additional file 1 for a rigorous derivation)
Additional file 1. The detailed derivation of relative entropy rate and the associated Fisher information matrix.
Format: PDF Size: 203KB Download file
This file can be viewed with: Adobe Acrobat Reader
Thus, from a practical point of view, RER is an observable of the stochastic process which can be computed numerically as an ergodic average, requiring only the knowledge of the propensity functions and the stoichiometric matrix (ν)_{i,j}. Nevertheless, in order to carry out the sensitivity analysis in the parameter vector θ, the computation of RER for different ε’s is necessary which can be computationally challenging for highdimensional parameter spaces. Thus, a sensitivity analysis methodology which does not depend on ε’s –such methods are called “gradientfree”– is desirable and is developed next.
Pathwise Fisher information matrix
Even though not directly evident from (4), a Taylor series expansion of RER in terms of ε reveals that RER is locally a quadratic function of the parameter vector . Indeed, RER is nonnegative when ε ≠ 0 and equals to zero when ε = 0 thus the linear term in the Taylor expansion is zero. Therefore, RER is written –under smoothness assumptions on the propensity functions with respect to the parameter vector θ– as [19],
where is a K × K matrix that can be considered as a pathwise analogue for the steady state Fisher Information Matrix (FIM). Similarly to the steady state FIM for parametrized distributions [33], is the Hessian of the RER which geometrically corresponds to the curvature around the minimum value of the relative entropy rate. The pathwise FIM contains up to third order accuracy all the sensitivity information for the path distribution at point θ for any perturbation direction ε, therefore, the computation of the FIM is sufficient up to third order for the evaluation of all the local sensitivities of the path distribution around the parameter vector θ. Moreover, an explicit formula for the pathwise FIM is given by (see Additional file 1 for a derivation)
The implications of this explicit formula are twofold. First, it reveals that for many typical reaction networks the FIM has a special blockdiagonal structure which reflects the parameter interdependencies and it is discussed in detail below. Second, the FIM is based on the propensity functions as well as on their derivatives which are known – actually, they define the process– thus the FIM, similarly to RER, is numerically computable as an observable of the process. Subsequent sections present various strategies to numerically estimate both the RER and the FIM in an efficient fashion.
Furthermore, the pathwise FIM, , can be used not only for the estimation/approximation of RER via (5) but also to infer intrinsic knowledge for the system’s sensitivities [16,37]. In general, the spectral analysis of the FIM reveals the (local) comparatively most/least sensitive directions of the system around θ. Indeed, by ordering the eigenvalues of the FIM as
it can be inferred that the most sensitive direction corresponds to the eigenvector with eigenvalue while the least sensitive direction corresponds to the eigenvector with eigenvalue . Additionally, the FIM is one of the most useful tools for optimal experimental design. Many of the optimality criteria such as Doptimality where the determinant of the FIM is maximized or Aoptimality where the trace of the inverse of the FIM is minimized are based on FIM, [37].
In the same direction, robustness of the system to parameter perturbations or errors as well as parameter identifiability can be studied utilizing spectral analysis of the FIM. For instance, parameter identifiability is satisfied when all the eigenvalues of the FIM are above a given threshold, [18].
Sensitivity analysis at the logarithmic scale
In many biochemical reaction networks, the model parameters differ by orders of magnitude and a reasonable option for carrying out sensitivity analysis is to perform perturbations which are proportional to the parameter magnitude. This can be carried out by perturbing the logarithm of the model parameters instead of the parameters itself. Using the chain rule ∇_{logθ}f(θ) = ∇_{θ}f(θ).∇_{logθ}θ = θ.∇_{θ}f(θ) where ‘. ’ is defined as the element by element multiplication (i.e., (a.b)_{k} = a_{k}b_{k}, k = 1, …, K), we obtain the logarithmicallyscaled FIM:
where is given by (6). Similarly, the logarithmic perturbation for the RER is carried out using the perturbation vector θ.ε instead of ε. Notice that (5) continues to be valid for the logarithmic scale, i.e.,
Linking relative entropy and observables
As we discussed in the previous sections, relative entropy provides a mathematically elegant and computationally tractable methodology for the parameter sensitivity analysis of complex, stochastic dynamical systems. Such results focus on the sensitivity of the entire probability distribution, either at equilibrium or at the pathspace level, i.e., for the entire stationary timeseries. However, in many situations of chemical and biological networks, the interest is focused on observables such as mean populations, population correlations, population variance as well as pathspace observables such as time autocorrelations and extinction times. Therefore, it is plausible to attempt to connect the parameter sensitivities of observables to the relative entropy methods proposed here. Indeed, relative entropy can provide an upper bound for a large family of observable functions through the Pinsker (or CsiszarKullbackPinsker) inequality, [33].
More precisely, for any bounded observable function f, the Pinsker inequality states that
where ·_{∞} denotes the supremum (here, maximum) of f. An obvious outcome of this inequality is that if the (pseudo)distance between two distributions defined by is controlled, then the error between the two distributions is also controlled for any bounded observable. In the context of sensitivity analysis, inequality (9) states that if the relative entropy is small, i.e., insensitive in a particular parameter direction, then, any bounded observable f is also expected to be insensitive towards the same direction. In this sense, (9) can be viewed as a “conservative” –but not necessarily sharp– bound for the parametric sensitivity analysis of observables, including pathdependent observables such as longtime averages and autocorrelations.
From a practical perspective, (9) can be used as an indicator that suggests –even in the presence of a very highdimensional parameter space– which are the insensitive parameter directions for observables of stochastic dynamical systems. The leastsensitive directions can be verified computationally and we present in the sequel two examples of this practical strategy in the p53 and the EGFR models. More generally, determining insensitive directions in parameter space can by particularly useful in multiparameter models in systems biology characterized by “sloppiness”, i.e., when most model parameters allow for a vast range of perturbations without affecting the dynamics, [38]. As we concretely show in the EGFR example, our methodology can easily demonstrate and quantify such properties in stochastic dynamics through the use of the spectral analysis of the pathwise FIM, even if the models include a very large number of parameters.
Remark 1
We note that in order to carry out such an analysis in a mathematically rigorous manner which includes both sensitive and insensitive parameters, we need to require that the norm ·_{∞} in (9) is controlled. For instance, typical observables in biochemical reaction networks are the number of molecules for each species, hence f(x) = x. Thus, for reaction networks where the population size is large, the Pinsker inequality (9) will provide a bound that may not be sharp. In fact, it was recently shown that there are alternative bounds which are expected to be practically more useful than (9) in the sense that they provide sharp upper bounds for observables in terms of the relative entropy, [39]. These sharp bounds rely on the variational representation of the relative entropy and the existence of an explicit minimizer for the upper bound. Furthermore, [39] combined with our work on RER and pathwise FIM suggests new mathematical questions towards pursuing practical sharp upper bounds involving RER and pathwise FIM. In view of these comments, we conclude that (9), (a) constitutes a theoretical indicator that relative entropy is a reliable tool for sensitivity analysis and more generally for uncertainty quantification; (b) from a practical perspective, it is capable to rule out insensitive directions in parameter space, which in turn provides a significant advantage in the study of “sloppy” multiparameter models.
Blockdiagonal structure of the pathwise FIM
In chemical reaction networks, reactions typically depend only on a small subset of the parameter vector. Mathematically, this is described as
where while L_{j} ≪ K is the number of involved parameters in reaction R_{j}. Using (6), it can be shown that this parametric dependence of the propensities is directly reflected on the pathwise FIM. Indeed, after grouping the reactions into subsets in such a way that each subset contains the minimum number of reactions having common parameters, the pathwise FIM – upon rearrangement of the parameter vector– becomes a blockdiagonal matrix. The pathwise FIM is then written as
where are block matrices. The block matrices which are defined by the reaction subsets with the same parametric dependence are easily obtained by creating a graph whose nodes are the reactions and the parameters while the edges are their dependences. Then, the parameter nodes contained in a connected subgraph define a parameter subset which in turn corresponds to a block of the FIM. An illustration of this procedure is shown in Figure 1 where a reaction network with M = 9 reactions and K = 7 parameters is plotted. The parametric dependencies of the reactions are shown in the left panel where 4 subgroups of parameters are defined based on the graph connectivity. The resulting blockdiagonal structure of the FIM is shown on the right panel of Figure 1.
Figure 1. The graph representation (left panel) of the dependencies between the reactions (left column) and the model parameters (right column) as well as the corresponding blockdiagonal structure of the FIM (right panel). The grouping of the parameters is carried out by noting the connected parts of the graph. In this example K = 7 while the largest dimension of the blocks is L = 3.
Before proceeding with the theoretical computation of the FIM for various wellknown classes of biochemical reaction networks, we list some of the implications of this simplified structure of the FIM in sensitivity analysis and elsewhere.
(i) The sparsity of the FIM is proportional to the parametric decoupling between the reactions. Knowing a priori the zero elements of the FIM, there is no need to numerically compute them. It is clear that the computation cost for each sample drops from O(K^{2}) to O(KL) where L is the largest dimension of the block matrices.
(ii) The inverse of the FIM is also blockdiagonal and each block of the inverse FIM is the inverse of the respective block. This fact allows us in the parameter estimation problem to easily evaluate the lower bound of the variance, at least for the completedata case [40], i.e., obtain CramerRao bounds, [41,42] which are given by the diagonal elements of the inverse of the FIM.
(iii) Relation (11) implies that that optimality criteria in optimal experiment design, [37], are significantly simplified. For example, the determinant employed in the Doptimality test is given by the relation , while the trace of the inverse of the FIM utilized in the Aoptimality reduces to .
(iv) Given that parametric identifiability is characterized by the magnitude of the eigenvalues of the FIM, e.g. a zero eigenvalue corresponds to a nonidentifiable direction in parameter space, [18,43], then the blockdiagonal structure (11) can provide additional insights in parameter identification. For instance, the identifiability of the parameters of the group that corresponds to the ith block, A_{i}, will be increased if the smaller eigenvalues of the ith block can be increased. On the other hand, if the determinant of the ith block equals to zero then at least one of its eigenvalues is zero and thus the corresponding linear combinations of parameters are nonidentifiable. Similarly, the robustness of the system to perturbations of the parameters of the ith group will be increased if there is a way to decrease the larger eigenvalues of the ith block.
Overall, we note that extracting useful information regarding model parameters can be performed for each block of the pathwise FIM independently. Next, we discuss two specific examples of biochemical reaction networks where the explicit calculation of the blockdiagonal FIM is demonstrated.
Reactions with independent reaction constants (mass action kinetics)
An important class of wellmixed reaction networks take the general form “” where A_{j} and B_{j} are the reactant species while α_{j} and β_{j} are the respective number of molecules needed for the reaction. The reaction constant, θ_{j}, is the parameter of the jth reaction. The propensity function for the jth reaction is given as the product between a rate constant and a function of the current state x:
Typically, which stems from the law of mass action, however, it can take different forms depending on the modeling of the physical process. This reaction network has K = M parameters, while each propensity depends only on one parameter, i.e., L_{j} = 1 in (10) for j = 1, …, M. The (k, l)th element of the FIM in the logarithmic scale is explicitlyl given by
where μ^{θ} is the stationary distribution of the stochastic process. Furthermore, it holds that where δ(·) is the Dirac function, therefore the pathwise FIM is a diagonal matrix with elements given by
This result demonstrates that the sensitivity of a reaction constant is proportional to the equilibrium average of the respective propensity function. Moreover, due to the diagonal form of the FIM, it is straightforward to carry out the eigenvalue analysis and infer the most/least sensitive directions of the reaction network: the eigenvalues of the FIM are its diagonal elements while the eigenvectors are the standard basis vectors of . Hence, the most (respectively least) sensitive parameter is obtained from the largest (respectively smallest) diagonal element of the FIM. Furthermore, (14) demonstrates that the (local) robustness of the reaction network to a specific parameter is inversely proportional to the mean propensity of the corresponding reaction. Another observation stemming from the diagonal structure of the pathwise FIM is that each rate constant can be estimated from timeseries data independently from the other rate constants. This observation has been already pointed out and discussed in the context of maximum likelihood estimation for the completedata case ([40], Sec. 10.2).
Additionally, the mean firing rate of a reaction is equal to the mean propensity. Hence, it can be stated that the parameters that correspond to the faster reactions, i.e., to reactions with larger mean firing rate, are more sensitive in a pathwise entropy sense. It should be noted, however, that not all observables are sensitive to the parameters that correspond to the faster reactions and there are examples (see the protein productiondegradation model in the Results section) where steady state observables such as the equilibrium distribution remain insensitive to specific perturbation directions even though their mean propensity may be increased. Finally, we would like to remark that even though trivially holds true, the diagonal elements of the FIM are not linear functions of the corresponding reaction constants since the steady state distribution μ^{θ}, depends also on the parameter vector θ. In fact, high reaction constants do not necessarily imply large mean propensities and hence a more sensitive parametric dependence. This is specifically due to the mean value in (14) and as an illustrative example we refer to the simple protein productiondegradation model (e.g., compare (24) and (27)).
MichaelisMenten kinetics
Another class of reaction networks is the MichaelisMenten kinetics. In its simplest form (e.g., singlesubstrate reaction without intermediate), this chemical network contains a reaction between species A and B (i.e., A → B) with propensity function given by
This reaction subnetwork which is derived under a quasisteadystate assumption is one of the bestknown models of enzyme kinetics in biochemistry [44]. The parameter θ_{k} (usually denoted by V_{max}) represents the maximum rate achieved by the system, at maximum (saturating) substrate concentrations while the Michaelis constant (usually denoted by K_{m}) is the substrate concentration at which the reaction rate is half the maximum value. The propensities of this MichaelisMenten subnetwork depend on two parameters (L_{k} = 2 in (10)) thus the corresponding FIM block is a 2 × 2 matrix. The elements of the FIM matrix are given by
for the kth row while the k^{′}th row is given by
In general, biochemical reaction networks may have significantly more complex propensities, nevertheless, the computation of the FIM follows exactly the same calculation lines for any propensity function.
Strategies for the statistical estimation of RER and FIM
Previous sections introduced and justified RER and FIM as appropriate observables for measuring the sensitivity analysis of the reaction network’s parameters in longtime dynamics. This section presents strategies on how to efficiently estimate these quantities as ergodic averages of the underlying stochastic process.
Unbiased statistical estimators
Since the stationary distribution, μ^{θ}, is usually not known, both FIM and RER should be estimated numerically as ergodic averages. Indeed, the statistical ergodic estimator for RER is given by
where Δt_{i} is an exponential random variable with parameter given by the total rate, , while is the total simulation time. The sequence is the embedded Markov chain with transition probabilities from state x_{i} to state x_{i+1} is given by the ratio . The weight Δt_{i}, which is the waiting time at state x_{i}, is necessary for the unbiased estimation of the observable, [45]. Similarly, the unbiased estimator for the FIM is
Noticing that the computation of the local propensity functions for all j = 1, …, M is needed for the simulation of the jump Markov process when Monte Carlo methods such as SSA [45] is utilized, the computation of the perturbed transition rates is the only additional computational cost for the numerical RER while the additional cost for the estimation of the FIM is the computation of the derivatives of the propensities. Algorithm 1 summarizes the numerical computation of RER and FIM, employing the SSA for the simulation of the jump Markov process.
Accelerated statistical estimators
A typical feature of biochemical systems is that the modeled reaction network is large with hundreds or thousands of reactions and different time scales stemming from the orders of magnitude difference between the reaction rates and/or between the species concentrations, making the SSA extremely slow. A large number of multiscale approximations of the original SSA have been developed in order to handle such issues resulting to accelerated simulation algorithms. For example, meanfield approximation ignores the fluctuations of the stochastic process and yields a deterministic system of ordinary differential equations (ODE) for the mean population of the species [46,47]. Stochastic corrections to the meanfield model such as stochastic Langevin [25] and linear noise approximation [48] can be applied in order to improve the accuracy of the simulation. An alternative approximation of the jump Markov process is the tauleap method proposed by Gillespie [24] where a batch of events occurs at each timeincrement, τ. Several improvements of the basic tauleap algorithm on how to select adaptively the τ[49] or avoiding negative populations [27,50] have been proposed, however, their performance is heavily modeldependent.
Algorithm 1 SSAbased numerical computation of RER and FIM
In this subsection, we propose such approximations in order to efficiently compute the FIM and/or RER observables, while maintaining controlled bias in the statistical estimators. As an illustration, we present the wellknown meanfield approximation. The popularity of the meanfield modeling stems from their computational efficiency. To proceed, the stochastic process can be written as
where x(t) is the deterministic part (mean) of the process, ξ(t) is the stochastic zeromean part while η is the amplitude of the stochastic term. The amplitude of the stochastic term is proportional to the inverse square root of the reactant populations [25,48,51]. Thus, for large populations, the fluctuations of the timeevolving species populations become vanishingly small compared to the deterministic contributions. Consequently, the dominant part of the process is the deterministic term whose dynamics are governed by the ODE system
This ODE system is also known as reaction rate equations [25]. Restricted for simplicity to the special case with independent rate constants for each reaction, the diagonal elements of the FIM are approximated using (19) as
Typically, such ODE system is solved using an adaptive timestep numerical integrator up to final time . Thus, for large species populations (S_{i}≫1), the following numerical estimator for the FIM’s diagonal elements is obtained:
Relation (22) suggests an algorithm similar to Algorithm 1 for the numerical computation of the FIM where instead of SSA, an ODE solver is employed.
Remark 2
Multiscale approximations are usually valid for large populations and relatively simple systems which do not exhibit complex dynamics such as bistability or intermittency. Indeed, large deviation arguments [52] or even explicitly available formulas for escape times [53] demonstrate that stochastic approximations cannot always capture correctly exit times, rare events, strong intermittency, etc. even in relatively simple systems. However, in order to simulate large biochemical systems there is often no other alternative than such approximate models, which nevertheless need to be employed with the necessary caution.
Remark 3
In biochemical systems, we are interested not only in the steady state, i.e., the stationary distribution or timeseries, but also in the transient regime, e.g. signaling phenomena. The extension of the proposed sensitivity analysis method to the transient regime is justified by the fact that the timenormalized relative entropy can be also decomposed as a sum of simple integrals [33] which results to the fact that the statistical estimators (17) and (18) remain valid. In a subsequent section we present an example of a biochemical system (EGFR) which exhibits transient behavior, and where the proposed sensitivity analysis tools are tested and validated. The rigorous mathematical derivation of the relative entropy rate for the transient regime is out of the scope of this publication and a dedicated mathematical article on the timedependent relative entropy rate will follow.
Results and discussion
A simple protein production/degradation model
We first consider an elementary stochastic model for protein production and degradation, [54], which is also a component of more complex models for gene regulatory networks, [55]. In this simplified model, the protein is produced at a constant rate k_{1}, while it is degraded with rate k_{2}, corresponding to the reactions
Accordingly, the corresponding propensity functions for the current state x = x are:
We consider this simple stochastic model due to the available analytic representations of the steady state (equilibrium) distribution, timedependent moments and autocorrelations, ([46], Sec. 7.1). Consequently, we can both illustrate the proposed pathwise sensitivity analysis, as well as compare it to the standard equilibrium FIM, revealing concretely differences between the two approaches.
The equilibrium distribution, μ^{θ}, of this simple network is a Poisson distribution with parameter . Therefore, the equilibrium FIM for the parameter vector θ = [k_{1}, k_{2}]^{T} is given in logarithmic scale by
On the other hand, the pathwise FIM is computed via (14):
where we used that
The complete calculations can be found in the Additional file 2. Some of the implications of the differences between these two FIMs are discussed next.
Additional file 2. The calculation of equilibrium and pathwise FIMs for the protein production/degradation model.
Format: PDF Size: 121KB Download file
This file can be viewed with: Adobe Acrobat Reader
First, we observe that the equilibrium FIM, (25), is singular, i.e., one of the eigenvalues is zero. We readily see that in the parameter direction defined by the corresponding eigenvector, i.e., when the parameter ratio, , remains constant, the system is expected to be insensitive, at least with respect to the equilibrium distribution. Clearly, this is a fact verified directly from the Poisson equilibrium distribution μ^{θ} which depends only on the ratio. On the other hand, the pathwise FIM, (26), is not singular and all the directions are equally sensitive. This fact suggests that observables for dynamic quantities may be sensitive not only to parameter ratio perturbations but also to other parameter perturbations. Indeed, one such example is the stationary autocorrelation function, which in the case of the simple protein production/degradation model is explicitly given by ([46], Sec. 7.1),
where < ·, · >_{s} denotes stationary averaging. Based on this formula, it is obvious that the autocorrelation function is also sensitive to k_{2}, in addiction to the ratio . This example demonstrates that in contrast to the pathwise FIM, the equilibrium FIM is inadequate to fully capture the dynamic properties of the process. Moreover, the pathwise FIM depends linearly only on k_{1}, which shows that the reaction rate constants and propensity functions in (24) alone, can be misleading in the assessment of parametric sensitivity. Contrary, the mathematically correct equilibrium averaging of the propensities, i.e., (14) can lead to a completely different outcome, as can be readily seen when we compare (24) and (27).
In terms of parameter identifiability, the fact that one of the eigenvalues of (25) is zero implies that that the twodimensional parameter vector of the system is nonidentifiable. Indeed, the asymptotic normality of the maximum likelihood estimators, [41,42], states that their variance (also a lower bound according to the Cramer Rao theorem), which determines parameter identifiability of k_{1} and k_{2}, is the reciprocal of the eigenvalues of (25). A straightforward calculation involving the eigenvectors of (25) shows that the only identifiable parameter is the ratio of the reaction constants appearing in (25). Therefore parametric inference for both parameters from equilibrium data is not possible. On the other hand, the pathwise FIM (26) is not singular, which readily implies that both parameters can be identified through (complete) timeseries data, provided that k_{1} ≠ 0. Summarizing, this birth/death model is an example where equilibrium sampling is not enough for the identifiability of all the parameters, however, if dynamics data are available and are taken into account then all the parameters become identifiable as pathwise FIM asserts.
The p53 gene model
The p53 gene plays a crucial role for effective tumor suppression in humans as its universal inactivation in cancer cells suggests [28,56,57]. The p53 gene is activated in response to DNA damage and gives rise to a negative feedback loop with the oncogene protein Mdm2. Models of negative feedback are capable of oscillatory behavior with a phase shift between the gene concentrations. Here, we perform sensitivity analysis to a simplified reaction network between three species, p53, Mdm2precursor and Mdm2 introduced in [28]. The model consists of five reactions and seven parameters provided in Table 1. The nonlinear feedback regulator of p53 through Mdm2 takes place in the second reaction while the remaining four reactions fall in the special class where each reaction depends on one parameter. Due to these mechanisms a nontrivial steady state regime exists and can be characterized by random oscillations, see for instance Figure 2. The proposed sensitivity methodology is directly applicable, and the corresponding pathwise FIM, see (13) and Figure 1, consists of 5 diagonal blocks with respective size 1 × 1, 3 × 3, 1 × 1, 1 × 1, 1 × 1. Furthermore, the sensitivity analysis of this model has been performed earlier in [18] based on a linear noise approximation. Here, we present a detailed comparison between the two sensitivity analysis methodologies, since the one proposed here does not involve any approximation of the stochastic network dynamics.
Table 1. The reaction table with x corresponding to p53,y_{0 }to Mdm2precursor while y corresponds to Mdm2
Figure 2. Molecule concentration of p53, Mdm2precursor and Mdm2. Concentration oscillations as well as time delays (phase shifts) between the species are present due to the negative feedback loop. Furthermore, the concentration of p53 periodically approaches zero and since negative concentrations are not allowed, the stochastic characteristics of p53 are far from Gaussian.
Figure 2 shows the timeseries of the species for the parameter values in Table 2. Evidently, oscillatory behavior is observed at this parameter regime, where persistent random oscillations occur, ranging between high and low populations. On the other hand, the frequency of the oscillations is less variable as it has been already reported both experimentally and numerically [28]. Another interesting observation is that the concentration of p53 species usually attains the lower bound of its admissible value (populations cannot be negative) which results in stochastic effects far away from Gaussianity, as can be readily seen also in Figure 2.
Table 2. Parameter values for the p53 model
Proceeding, we denote by θ = [b_{x}, a_{x}, a_{k}, k, b_{y}, a_{0}, a_{y}]^{T} the parameter vector. The numerical estimator for RER as well as for the pathwise FIM in the logarithmic scale are computed utilizing Algorithm 1. Logarithmic sensitivity analysis is preferred because the range of the parameters values varies by orders of magnitude as can be seen in Table 2. The upper plot in Figure 3 shows the RER as a function of time for various perturbations. Viewing RER as an observable, it is striking the speed of relaxation of the estimator. Within two or three oscillation periods, RER has been converged to its value even though the three species have significant oscillations and stochasticity, as Figure 2 shows. A primary reason for the fast relaxation is the numerical estimator of RER where the summation is over all reactions even though only one reaction takes places at each jump (see (18)). Having the important property of fast convergence, global sensitivity analysis, where not only a point of the parameter regime but also large subsets of the parameter space, can be efficiently performed, [15]. The lower panel of Figure 3 shows the RER when only one of the parameters are perturbed by +10% or by 10%. Additionally, the RER computed from the FIM, utilizing (5), is also provided. The FIM approximation of RER is a second order approximation in terms of ε, hence the computation of FIM is typically enough to fully resolve the local sensitivities of a model. Evidently, the most sensitive parameters here are b_{x} and a_{k} while the least sensitive parameters are a_{x} and k.
Figure 3. RER in time (upper panel) for the parameter perturbation of b_{x} (blue), k (green) and a_{y} (red) by +10% (i.e.,ε_{0} = 0.1) as well as RER for various perturbation directions (lower panel) computed either directly (blue and red bars) or based on FIM (green bars). Direction e_{k} corresponds to perturbation of parameter θ_{k}.
Comparison to the LNAbased sensitivity approach
In [18], the authors suggested a linear noise approximation (LNA) for the stochastic evolution around the nonlinear meanfield equation, and based on this approximation a system of ODEs is derived for the mean and the covariance matrix of the approximation process. Since the noise of LNA is Gaussian, the mean and the covariance matrix contains all necessary information regarding the approximate stochastic model. Then, the associated FIM is derived and based on it, the sensitivities for each parameter are computed. Although there are regimes where this approximation is applicable (short times, high populations, systems with a single steady state, etc.), for systems with nontrivial longtime dynamics, e.g. metastable, it is not correct as large deviation arguments [52] and explicit formulas for escape times [53] show. Similar issues with nongaussianity in the longtime dynamics arise in stochastic systems with strongly intermittent (pulselike) or random oscillatory behavior [58]. In the p53 model considered in [18] which had the same parameter values as here, Figure 2 reveals that the timeseries of the p53 populations persistently fluctuate between high and low values, thus the LNA approximation may not be accurate at least when the concentration of the species is very low.
At first pass, when the parameters are grouped into two classes depending on their sensitivities, the two sensitivity approaches produce qualitatively similar results. Indeed, by visual inspecting the lower plot of Figure 3 in the current publication and Figure three in [18], the (more) sensitive parameters in both methods are b_{x}, b_{y}, a_{k}, a_{0}, a_{y} while the practically insensitive parameters are a_{x}, k. However, upon closer inspection, the two methods produce different results. Figure 4 shows the proposed FIM (left) based on the exact (without any approximations) pathwise relative entropy theory, as well the FIM proposed in [18] which is derived from the LNA of the reaction system. The results are completely different and the proposed pathwise FIM is sparse as expected. A striking difference between the two sensitivity approaches is that the sensitivity of parameter b_{x} in our proposed method is relatively high compared to the other parameters while the sensitivity of b_{x} in the LNAbased method is at least one order lower compared to the other parameters, see Figure 4 (dark blue) and also compare Figure 3 of this publication and Figure three in [18].
As a means of comparison between the methods, we perturb b_{x} as well as b_{y} by the same amount and observe the Power Spectral Density (PSD), i.e., the square of the absolute value of the Fourier transform of each species’ timeseries. Given the sustained random oscillations observed in the p53 model, see Figure 2, the PSD is a suitable observable since it identifies the dominant periodicities and corresponding amplitudes in stationary timeseries, [41]. Using the Pinsker inequality (9) as a guideline, we expect that the observable will not be sensitive to the least sensitive directions of the FIM, therefore, we focus on the most sensitive directions of the FIM identified in Figure 4. Figure 5 shows the averaged PSD for the three species of the model for the unperturbed case (black lines), the perturbation of b_{x} only (blue lines) as well as the perturbation of b_{y} only (yellow lines). One hundred realizations were used for the averaging procedure while the perturbation strength was +20%. The L_{1}norm, i.e., the integral of the absolute difference, between the unperturbed PSD and the b_{x}perturbation is 8.56 · 10^{5} while the L_{1}norm between the unperturbed PSD and and the b_{y}perturbation is 4.32 · 10^{5} which is about the half value. Hence, the averaged PSD –primarily the amplitude of the oscillations– is more sensitive to perturbations of b_{x} rather than to perturbations of b_{y} as our sensitivity analysis method predicted while the LNAbased method suggested the reverse order of sensitivity. We note that the choice of the L_{1} norm for the PSDs is justified since it describes the total energy (power) of the timeseries, when viewed as a signal, [41]. An explanation of the performance of the LNAbased sensitivity analysis stems from the fact that the p53 species does not have Gaussian noise when the population is close to zero, and which can indeed occur frequently, see Figure 2 (blue line). Additionally, notice that both b_{x} and b_{y} affect the concentration of p53 explicitly or implicitly through the associated reactions thus their sensitivities are heavily biased due to the wrong statistical approximation of the p53 species. Moreover, we note that other observables, e.g., the arg max (location of the maximum) of the PSD, can be sensitive to perturbations of b_{y}, see Figure 5. This observation is not contradictory to the findings of the proposed sensitivity methodology for two reasons: first, even though the Pinsker inequality (9) points towards the right direction regarding sensitive parameters for observables, it is only an upper (possibly crude) bound. Second, even though b_{x} is more sensitive in absolute value than b_{y} in terms of the proposed pathwise FIM, both b_{x} and b_{y} sensitivities have the same order of magnitude (see Figure 4) therefore there should exist observables which are also sensitive to b_{y} perturbations. Finally, for the sake of completeness, we report the observable values for the insensitive parameters. The L_{1}norm between the unperturbed PSD and the a_{x}perturbation is 9.03 · 10^{3} which is approximately two orders of magnitude less than the L_{1}norm for the sensitive parameters. Similar results hold for the other insensitive parameter, k. Thus, the outcome of the Pinkser inequality (9), i.e., that small RER values imply relatively insensitive parameters for the observables, is also numerically verified.
Figure 5. Power Spectral Densities (PSD) of the timeseries of the species in the p53 model for the unperturbed parameter regime (black), when b_{x }is perturbed by +20% (blue) as well as when b_{y }is perturbed by +20% (yellow). The value and the position of the prominent peak of the PSD are related with the amplitude and the frequency of the species oscillations. The visual comparison between the averaged PSDs suggests that the spectral properties are more sensitive to b_{x} than to b_{y}.
Epidermal growth factor receptor model
The EGFR model is a wellstudied system describing signaling phenomena of (mammalian) cells [2931]. As its name suggests, EGFR regulates cell growth, survival, proliferation and differentiation and plays a complex and crucial role during embryonic development and in tumor progression [60,61]. In this paper, we study the reaction network model for the dynamics of EGFR developed by Schoeberl et al. [31] which consists of 94 species and 207 reactions. Figure 6 presents the EGFR reaction network in an abstract level. Initially, the extracellular binding of EGF with the EGF receptors induce receptor dimerization. Then, two principal pathways, Shcdependent and Schindependend, are initiated leading to activation of RasGTP. Subsequently phosphorylation of MEK kinase through the activation of Raf kinase occurs leading to the phosphorylation of ERK kinase which regulates several proteins and nuclear transcription factors inside the cell. The detailed graphical description of the reaction network can be found in the Figures one & two of supplementary information in [31]. For completeness, all the reactions along with their rates are provided in the Additional file 3 of this publication.
Figure 6. Building blocks of the EGFR reaction network. Each module communicates with the adjacent modules through few species only. Additionally, with the exception of the first module, all the others are double, one external (i.e., outside the cell surface) and one internal.
Additional file 3. This file contains in plain text the reactions and the reaction constants of the EGFR model.
Format: TXT Size: 13KB Download file
The propensity functions for the reactions R_{1}, …, R_{97}, R_{100}, …, R_{207} of the EGFR network are written in the general form
with the exception of reaction pair R_{98}, R_{99} where their propensity functions are governed by the MichaelisMenten kinetics
where x is the current state of the reaction system while A_{j} corresponds to the reacting species. The parameter vector contains all the reaction constants,
with all values provided in the Additional file 3. Due to the specific values of the reaction constants as well as the initial population of the species (see Table 3), the firing rates between reactions differ by many orders of magnitude giving rise to a highly stiff network. Therefore, even though there are some stochastic implementations, [27], here for the purposes of RER and FIM calculations, we adopt the meanfield approximation discussed in the accelerated estimators subsection. We solve the derived system of ODEs with Matlab’s routine ode15s and compute the FIM at the steady state regime which corresponds to the time interval [500,700]. The completion of the internalization process needs about 500 seconds. It should be noted here that even though the simulation of the EGFR is performed utilizing a deterministic approximation model, the computed pathwise FIM has been derived from the stochastic network, i.e., (13). This approximation is expected to be valid in the sense of (19) due to the large populations considered here. Overall, the computed FIM is a sparse matrix and measures efficiently the sensitivities of the stochastic model in a gradientfree manner.
Table 3. Initial population of the species for the EGFR network
The upper plot of Figure 7 shows the diagonal elements of the FIM in descending order computed at the steady state regime. We report our results in the format of Figure 7 in order to be able to accommodate the large number of parameters in the model. The kth diagonal element of the FIM corresponds to RER where the perturbation takes place only to the kth parameter (see (5)). Figure 7 (upper plot) in conjunction with Table S1 of the Additional file 4 fully describe the (local) sensitivities of the reaction network. Table S1 in Additional file 4 presents the reaction constants ordered from the most sensitive to the least sensitive parameter. Moreover, the FIM is diagonal –except a small 2×2 block associated with the MichaelisMenten reactions– therefore the diagonal elements correspond to the eigenvalues of the FIM. The sensitivity analysis depicted in Figure 7, demonstrates that most model parameters allow for a vast range of perturbations without affecting the dynamics. Furthermore, this robustness to variations in most parameters was also reported in the original, fully deterministic EGFR model in [31]. This is a feature shared by many multiparameter models in systems biology and which is known as “sloppiness”, [38]. Our methodology can easily demonstrate such properties in stochastic dynamics, as we can readily see in Figure 7, even if the models include a large number of parameters.
Figure 7. Diagonal elements of the FIM computed at the steady state regime (upper plot) and at the transient regime (lower plot). Note the changes in sensitivity and consequently the parameter identifiability. The parameter sensitivities differ by orders of magnitude.
Additional file 4. The ordering of the parameter sensitivities for the EGFR model.
Format: PDF Size: 65KB Download file
This file can be viewed with: Adobe Acrobat Reader
The previous discussion refers to the analysis of the EGFR model to the steady state regime. On the other hand, EGFR is a signaling model whose transient regime, in addition to the steady state, is of great interest. As discussed in Remark 3, we can justify the application of the RER and FIM sensitivity analysis in the transient regime. Therefore, we compute the proposed FIM at the time interval [0,10], using (22). The lower plot of Figure 7 shows the diagonal elements of the pathwise FIM in the transient regime while keeping the ordering of the parameters unchanged from the upper, steady state plot. The parameter sensitivity ordering is completely different meaning that the sensitivities are timedependent in the transient regime. For instance, the most sensitive parameters in the stationary regime correspond to the final products of the reaction network, however, in the time interval [0,10] these species have not been produced yet resulting to insensitive reaction constants. In terms of parameter identification and estimation, the timedependent sensitivities imply that in order to extract the maximum information content from the experimental data, we have to estimate the parameters drawing samples from different time intervals. These time intervals should be defined based on the respective sensitivity indices and selected in order to maximize the identifiability for each set of parameters.
Finally, the Pinsker inequality (9) suggests that insensitive parameters can be perturbed, even significantly, without affecting species concentrations or other observable. As an illustration of this fact, we present in Figure 8 the concentrations of various critical species of the EGFR model when the 140th (k_{65}) most sensitive parameter is perturbed (see Table S1 in Additional file 4). The rate constant k_{65} corresponds to a reaction of the Shcdependent pathway module. Solid blue lines correspond to the unperturbed parameter case while the dashed red lines correspond to the perturbed case where the perturbation is a multiplication by a factor of ten of k_{65}. We present the total number of (EGFEGFR*)2 binding species without Sch* (top, left panel) and with Sch* (top, middle panel) as well as RasGTP (top, right panel), total activated Raf or total Raf* (low, left panel), doubly phosphorated MEK or MEKPP (low, middle panel) and doubly phosphorated ERK or ERKPP. These species are important for the understanding of the system since the different modules of the EGFR reaction network communicate through them (see Figure 6). It is evident from Figure 8 that the various species concentrations remain unchanged to perturbations of the insensitive parameter k_{65} as it was expected from (9). Moreover, we notice that although the average populations become large, which implies that the maximum norm in (9) is also large, we still obtained robust results regarding k_{65}’s parameter sensitivity.
Figure 8. Timedependent concentration of various species of the EGFR network either for the unperturbed parameter vector (solid blue lines) or for the perturbed one (dashed red lines). The 140th most sensitive parameter (k_{65}) is tenfold increased and the species concentrations are not affected. For the least sensitive parameters such as k_{65}, we rigorously know from the Pinsker inequality (9) that they should not alter the concentration values or any other observable even when they are heavily perturbed.
Conclusions
In this paper, we applied and extended a recently proposed parametric sensitivity analysis methodology to complex stochastic reaction networks. This sensitivity analysis approach is based on the quantification of information loss along different parameter perturbations between timeseries distributions. This is achieved by employing the Relative Entropy Rate, which is directly computable from the propensity functions. A key aspect of the method is that we can derive rigorously an associated Fisher Information Matrix on pathspace, which in turn constitutes a gradientfree approach for parametric sensitivity analysis; as such it provides a significant advantage in stochastic systems with a large number of parameters. We demonstrated that the structure of the pathwise FIM revealed hidden parameter interdependencies between the reactions. The blockdiagonal structure of the FIM highlighted the sparsity of the matrix which resulted in further improvements in the computational efficiency of the proposed method. Therefore, parametric sensitivity analysis for highdimensional stochastic reaction systems becomes tractable since it is wellknown that in high dimensional stochastic systems sensitivity analysis techniques can involve estimators of very high variance, e.g. in finite difference methods and their recently proposed variants, which can present an overwhelming computational cost. Additionally, we proposed the use of multiscale numerical approximations of stochastic reaction networks in order to derive efficient statistical estimators for the FIM and implemented one such approximation (meanfield) in a highdimensional system.
The proposed pathwise sensitivity analysis method is tested and validated on three biological systems: (a) a simple protein production/degradation model where explicit solutions are available, (b) the p53 reaction network where quasisteady stochastic oscillations of the concentrations are observed and where multiscale stochastic approximations break down due to the persistent oscillations between low and high populations, and (c) a stochastic EGFR model which is an example of a highdimensional reaction network with more than 200 reactions and a corresponding number of parameters. In the EGFR reaction network, we combined the proposed pathwise FIM which has been derived from the stochastic network and the meanfield approximation which is used for the efficient estimation of the pathwise FIM. Moreover, our earlier rigorous analysis for the steady state regime [19] suggests suitable extensions in the transient regime which were tested and validated for the EGFR model. We will present the full rigorous theory in an upcoming publication.
Finally, we note that the relation between RER and various observables is not straightforward. However, we note that the path distribution contains all information regarding the process including the steady state and all timedependent observables: practically, our proposed sensitivity analysis represents a “conservative” sensitivity estimate in the sense that insensitive directions for the relative entropy on pathspace will yield insensitive directions for every observable. This latter statement can be justified mathematically through the Pinsker inequality (9) which was tested in the examples considered here. Based on these observations, the proposed sensitivity analysis methods can be deployed in complementary fashion with existing sensitivity analysis tools, as it can be used to narrow down the most sensitive directions in a system.
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
MAK and YP conceived the proposed sensitivity analysis methodology. YP conducted the numerical experiments. MAK and DGV designed and supervised the conducted research. All authors contributed to the preparation of the manuscript. All authors read and approved the final manuscript.
Acknowledgements
The work of MAK and YP was supported in part by the Office of Advanced Scientific Computing Research, U.S. Department of Energy under Contract No. DESC0002339 and DESC0010723. The work of DGV was supported in part by the Office of Advanced Scientific Computing Research, U.S. Department of Energy under Contract No. DEFG0205ER25702 and DESC0010549. The work of MAK was also supported in part by the European Union (European Social Fund) and Greece (National Strategic Reference Framework), under the THALES Program, grant AMOSICSS.
References

Barkai N, Leibler S: Robustness in simple biochemical networks.
Nature 1997, 387(6636):913917. PubMed Abstract  Publisher Full Text

Csete M, Doyle J: Reverse engineering of biological complexity.
Science 2002, 295(5560):16641669. PubMed Abstract  Publisher Full Text

Kitano H: Opinion  Cancer as a robust system: implications for anticancer therapy.
Nat Rev Cancer 2004, 4(3):227235. PubMed Abstract  Publisher Full Text

Donz A, Fanchon E, Gattepaille L, Maler O, Tracqui P: Robustness analysis and behavior discrimination in enzymatic reaction networks.

Hart Y, Antebi Y, Mayo A, Friedman N, Alon U: Design principles of cell circuits with paradoxical components.
Proc Nat Acad Sci USA (PNAS) 2012, 109(21):83468351. Publisher Full Text

Glynn P: Likelihood ratio gradient estimation for stochastic systems.
Commun ACM 1990, 33(10):7584. Publisher Full Text

Nakayama M, Goyal A, Glynn PW: Likelihood ratio sensitivity analysis for Markovian models of highly dependable systems.
Stochastic Models 1994, 10:701717. Publisher Full Text

Plyasunov S, Arkin AP: Efficient stochastic sensitivity analysis of discrete event systems.
J Comp Phys 2007, 221:724738. Publisher Full Text

Kim D, Debusschere B, Najm H: Spectral methods for parametric sensitivity in stochastic dynamical systems.
Biophys J 2007, 92:379393. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Rathinam M, Sheppard PW, Khammash M: Efficient computation of parameter sensitivities of discrete stochastic chemical reaction networks.

Anderson DF: An efficient finite difference method for parameter sensitivities of continuoustime Markov chains.
SIAM J Numerical Anal 2012, 50(5):22372258. Publisher Full Text

Sheppard P, Rathinam M, Khammash M: A pathwise derivative approach to the computation of parameter sensitivities in discrete stochastic chemical systems.

McGill JA, Ogunnaike BA, Vlachos DG: Efficient gradient estimation using finite differencing and likelihood ratios for kinetic Monte Carlo simulations.
J Comp Phys 2012, 231(21):71707186. Publisher Full Text

Liu H, Chen W, Sudjianto A: Relative entropy based method for probabilistic sensitivity analysis in engineering design.
J Mech Design 2006, 128:326336. Publisher Full Text

Lüdtke N, Panzeri S, Brown M, Broomhead DS, Knowles J, Montemurro MA, Kell DB: Informationtheoretic sensitivity analysis: a general method for credit assignment in complex networks.
J R Soc Interface 2008, 5:223235. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Majda AJ, Gershgorin B: Quantifying uncertainty in climate change science through empirical information theory.
Proc Natl Acad Sci 2010, 107(34):1495814963. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Majda AJ, Gershgorin B: Improving model fidelity and sensitivity for complex systems through empirical information theory.
Proc Natl Acad Sci 2011, 108(25):1004410049. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Komorowski M, Costa MJ, Rand DA, Stumpf MPH: Sensitivity, robustness, and identifiability in stochastic chemical kinetics models.
Proc Natl Acad Sci USA 2011, 108:86458650. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Pantazis Y, Katsoulakis M: A relative entropy rate method for path space sensitivity analysis of stationary complex stochastic dynamics.
J Chem Phys 2013, 138(5):054115. PubMed Abstract  Publisher Full Text

Gillespie DT: Exact stochastic simulation of coupled chemical reactions.
J Chem Phys 1977, 81:23402361. Publisher Full Text

Chatterjee A, Vlachos DG: An overview of spatial microscopic and accelerated kinetic Monte Carlo methods for materials’ simulation.
J Comput Aided Mater Design 2007, 14(2):253308. Publisher Full Text

Slepoy A, Thompson A, Plimpton S: A constanttime kinetic Monte Carlo algorithm for simulation of large biochemical reaction networks.
J Chem Phys 2008, 128(20):205101. PubMed Abstract  Publisher Full Text

Gibson MA, Bruck J: Efficient exact stochastic simulation of chemical systems with many species and many channels.
J Chem Phys 2000, 104:18761889. Publisher Full Text

Gillespie DT: Approximated accelerated stochastic simulation of chemically reacting systems.
J Chem Phys 2001, 115(4):17161733. Publisher Full Text

Gillespie DT: The chemical Langevin equation.
J Chem Phys 2000, 113:297306. Publisher Full Text

Rathinam M, Petzold LR, Cao Y, Gillespie DT: Stiffness in stochastic chemically reacting systems: The implicit tauleaping method.
J Chem Phys 2003, 119:1278412794. Publisher Full Text

Chatterjee A, Vlachos DG, Katsoulakis MA: Binomial distribution based tauleap accelerated stochastic simulation.
J Chem Phys 2005, 122:024112. PubMed Abstract  Publisher Full Text

GevaZatorsky N, Rosenfeld N, Itzkovitz S, Milo R, Sigal A, Dekel E, Yarnitzky T, Liron Y, Polak P, Lahav G, Alon U: Oscillations and variability in the p53 system.
Mol Syst Biol 2006, 2:0033. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Moghal N, Sternberg P: Multiple positive and negative regulators of signaling by the EGF receptor.
Curr Opin Cell Biol 1999, 11:190196. PubMed Abstract  Publisher Full Text

Hackel P, Zwick E, Prenzel N, Ullrich A: Epidermal growth factor receptors: critical mediators of multiple receptor pathways.
Curr Opin Cell Biol 1999, 11:184189. PubMed Abstract  Publisher Full Text

Schoeberl B, EichlerJonsson C, Gilles E, Muller G: Computational modeling of the dynamics of the MAP kinase cascade activated by surface and internalized EGF receptors.
Nat Biotechnol 2002, 20:370375. PubMed Abstract  Publisher Full Text

Kullback S: Information Theory and Statistics. New York: John Wiley and Sons; 1959.

Cover TM, Thomas JA: Elements of Information Theory.. Wiley Series in Telecommunications; 1991.

Kipnis C, Landim C: Scaling Limits of Interacting Particle Systems. Berlin, Heidelberg and New York: SpringerVerlag; 1999.

Avellaneda M, Friedman C, Holmes R, Samperi D: Calibrating volatility surfaces via relative entropy minimization.
Soc Sci Res Netw 1997.
Available at SSRN: [http://ssrn.com/abstract=648 webcite]

Dumitrescu ME: Some informational properties of Markov purejump processes.

Emery AF, Nenarokomov AV: Optimal experiment design.
Meas Sci Technol 1998, 9:864876. Publisher Full Text

Gutenkunst RN, Waterfall JJ, Casey FP, Brown KS, Myers CR, Sethna JP: Universally sloppy parameter sensitivities in systems biology models.
PLOS Comput Biol 2007, 3:18711878. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Chowdhary K, Dupuis P: Distinguishing and integrating aleatoric and epistemic variation in uncertainty quantification.
ESAIM Math Model Numerical Anal 2013, 47:635662. Publisher Full Text

Wilkinson DJ: Stochastic Modelling for Systems Biology.. Chapman & Hall; 2012.

Kay SM: Funtamentals of Statistical Signal Processing: Estimation Theory. Englewood Cliffs: PrenticeHall; 1993.

Wasserman L: All of Statistics: A Concise Course in Statistical Inference.. Springer; 2004.

Rothenberg T: Identification in parametric models.
CONOMETRICA 1971, 39:5770591. Publisher Full Text

Marangoni AG: Enzyme Kinetics: A Modern Approach.. Wiley Interscience; 2003.

Gillespie DT: A general method for numerically simulating the stochastic time evolution of coupled chemical reactions.
J Comp Phys 1976, 22:403434. Publisher Full Text

Gardiner C: Handbook of Stochastic Methods: for Physics, Chemistry and the Natural Sciences.. Springer; 1985.

van Kampen NG: Stochastic Processes in Physics and Chemistry.. North Holland; 2006.

Kurtz TG: The relationship between stochastic and deterministic models for chemical reactions.
J Chem Phys 1972, 57:2976. Publisher Full Text

Cao Y, Gillespie DT, Petzold LR: Efficient step size selection for the tauleaping simulation method.
J Chem Phys 2006, 124:044109. PubMed Abstract  Publisher Full Text

Tian T, Burrage K: Binomial leap methods for simulating stochastic chemical kinetics.
J Chem Phys 2004, 121:10356. PubMed Abstract  Publisher Full Text

Approximation of Population Processes, Society for Industrial and Applied Mathematics (SIAM). 1981.

Doering CR, Sargsyan KV, Sander LM, VandenEijnden E: Asymptotics of rare events in birth–death processes bypassing the exact solutions.

Hanggi P, Grabert H, Talkner P, Thomas H: Bistable systems: Master equation versus FokkerPlanck modeling.
Phys Rev A 1984, 29:371378. Publisher Full Text

Wilkinson DJ: Stochastic modelling for quantitative description of heterogeneous biological systems.
Nat Rev Genet 2009, 10:122133. PubMed Abstract  Publisher Full Text

Thattai M, van Oudenaarden A: Intrinsic noise in gene regulatory networks.
Proc Natl Acad Sci 2001, 98(15):86148619.
[http://www.pnas.org/content/98/15/8614.abstract webcite]
PubMed Abstract  Publisher Full Text  PubMed Central Full Text 
Prives C: Signaling to p53: breaking the MDM2p53 circuit.
Cell 1998, 95:58. PubMed Abstract  Publisher Full Text

Harris S, Levine A: The p53 pathway: positive and negative feedback loops.

Katsoulakis MA, Majda AJ, Sopasakis A: Intermittency, metastability and coarse graining for coupled deterministicstochastic lattice systems.
Nonlinearity 2006, 19(5):10211047. Publisher Full Text

Komorowski M, Zurauskiene J, Stumpf M: StochSens–Matlab package for sensitivity analysis of stochastic chemical systems.
Bioinformatics 2012, 28:731733. PubMed Abstract  Publisher Full Text

Sibilia M, Steinbach J, Stingl L, Aguzzi A, Wagner E: A strainindependent postnatal neurodegeneration in mice lacking the EGF receptor.
EMBO J 1998, 17:719731. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Kim H, Muller W: The role of the EGF receptor family in tumorigenesis and metastasis.
Exp Cell Res 1999, 253:7887. PubMed Abstract  Publisher Full Text