Parametric sensitivity analysis for biochemical reaction networks based on pathwise information theory

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 high-dimensional 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 time-series distributions. For this reason, we need to work on path-space, 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 rigorously-derived 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 gradient-free approach to quantifying parameter sensitivities. The structure of the FIM turns out to be block-diagonal, revealing hidden parameter dependencies and sensitivities in reaction networks. Conclusions As a gradient-free 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 quasi-steady 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 high-dimensional stochastic reaction network with more than 200 reactions and a corresponding number of parameters.


Background
The need of an intrinsic understanding of the interplay between complexity and robustness of biological processes and their corresponding design principles is well-documented, see for instance [1][2][3][4][5].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 low-dimensional stochastic processes, modeling well-mixed chemical reactions and biological networks.Some of the mathematical tools included log-likelihood methods and Girsanov transformations [6][7][8], 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 high-dimensional 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 non-Gaussian 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 [14][15][16][17][18] 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 time-series distribution.Typically, the space of all such time-series 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 long-time, stationary dynamics in path space, including for example bistable, periodic and pulse-like dynamics.Second, it is a gradient-free sensitivity analysis method suitable for high-dimensional 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 non-equilibrium steady state systems.In [19], we demonstrated these features by focusing on two classes of problems: Langevin particle systems with either reversible (gradient) or non-reversible (non-gradient) forcing, highlighting the ability of the method to carry out sensitivity analysis in non-equilibrium systems; and spatially extended Kinetic Monte Carlo models, showing that the method can handle high-dimensional 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), [20][21][22] and the next-reaction method [23], or by employing approximations such as mean field ODEs, tau-leap [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 block-diagonal 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 high-dimensional, 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,[24][25][26][27].
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 single-species birth/death model, with explicitly known formulas for the stationary and the time-dependent 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 [29][30][31].This is a high-dimensional system both in the number of variables and parameters, including 94 species and 207 reactions.Having a gradient-free method such as FIM 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 continuous-time jump Markov processes as well as we reveal the block-diagonal 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 well-mixed 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 N -dimensional 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 N -dimensional vector ν j correspond to the stoichiometry vector of j-th 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 j-th 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, {X(t)} t∈R+ is a continuous-time, time-homogeneous, jump Markov process with countable state space E ⊂ N N .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 a 0 (x) := M j=1 a j (x) while the transition probabilities of the process are determined by the ratio aj (x) a0(x) .We refer to Algorithm 1 for the details of the stochastic simulation.

Relative Entropy
Assume that two probability distributions (or more generally probability measures) P and P have corresponding probability densities p = p(x) and p = p(x).Then, the Relative Entropy or Kullback-Leibler divergence of P with respect to P is defined as [32,33] In a more general setting, relative entropy is defined as R P | P := log dP d P dP where dP d P is a function known as Radon-Nikodym derivative while the integration is performed with respect to the probability measure P, [34].A necessary condition for the relative entropy to be well-defined is that the Radon-Nikodym derivative exists which is satisfied when P is absolutely continuous with respect to P. 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 non-negative, (ii) it equals to zero if and only if P = P P-almost everywhere, and, (iii) R P | P < ∞ if and only if P and P are absolutely continuous with respect to each other.
From an information theory perspective, relative entropy quantifies the loss of information when P is utilized instead of P, [33].In other words, relative entropy quantifies the inefficiency of assuming an incorrect or perturbed distribution P instead of employing the true distribution P. 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 θ ∈ R K (i.e., a j (x) ≡ a θ j (x)) while the continuous-time jump Markov process X(t) t∈R+ 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., time-series distributions) are high-dimensional complex objects; for instance, if we consider the simpler discrete-time Markov chain case {Z n } n∈Z+ , defined by the transition probability density p(z, z ′ ), then, utilizing repeatedly the Markov property, the stationary path distribution of the time-series (z 0 , z 1 , . . ., z T ) is given by Proceeding, we consider another continuous-time jump Markov process { X(t)} t∈R+ defined by perturbing the propensity functions by a small vector ǫ ∈ R K .The corresponding steady state and path distributions of { X(t)} t∈R+ are denoted by µ θ+ǫ (x) and Q θ+ǫ [0,T ] , respectively.Let the two path distributions Q θ [0,T ] and Q θ+ǫ [0,T ] be absolutely continuous with respect to each other which is satisfied when a θ j (x) = 0 if and only if a θ+ǫ j (x) = 0 holds for all x ∈ E and j = 1, ..., M .Then, the Relative Entropy of the path distribution Q θ [0,T ] with respect to Q θ+ǫ [0,T ] is defined similarly to (1) as where . In fact, using the Girsanov's formula, we can obtain an explicit expression for the Radon-Nikodym 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 well-defined quantity, [36].As first proposed in [19], H Q θ | Q θ+ǫ is a suitable time-independent measure of sensitivity: it measures the rate of the loss of information due to an ǫ-perturbation of the model parameters, in the long-time, stationary dynamics regime of the stochastic process.Furthermore, RER admits an explicit formula given by (see Supplementary Information for a rigorous derivation) 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 high-dimensional parameter spaces.Thus, a sensitivity analysis methodology which does not depend on ǫ's -such methods are called "gradient-free"-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 ǫ ∈ R K .Indeed, RER is non-negative 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 in the parameter vector θ-as [19]: where F H (Q θ ) 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], F H (Q θ ) 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 Supplementary Information 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 block-diagonal 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, F H (Q θ ), 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) most/least sensitive directions of the system around θ. Indeed, by ordering the eigenvalues of the FIM as λ θ 1 ≥ ... ≥ λ θ K ≥ 0 , it can be inferred that the most sensitive direction corresponds to the eigenvector with eigenvalue λ θ 1 while the least sensitive direction corresponds to the eigenvector with eigenvalue λ θ K .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 A-optimality 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 logarithmically-scaled FIM: where F H (Q θ ) 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 path-space level, i.e., for the entire stationary time-series.However, in most simulations of chemical and biological networks, the main focus of interest is observables such as mean populations, population correlations, population variance as well as path-space 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 Csiszar-Kullback-Pinsker) 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 R Q θ | Q θ+ǫ 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 path-dependent observables such as long-time averages and autocorrelations.
From a practical perspective, (9) can be used as an indicator that suggests -even in the presence of a very high-dimensional parameter space-which are the insensitive parameter directions for observables of stochastic dynamical systems, based on the proposed relative entropy methods.The suggested least-sensitive parameter-space directions can be subsequently verified computationally.We provide two examples of this practical strategy in the p53 and the EGFR models discussed in the sequel.
Remark 1: However, we note that in order to carry out such an analysis in a mathematically rigorous manner, 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.Nevertheless, it is possible to derive a tighter inequality where the maximum norm of f in ( 9) is replaced by the standard deviation of f , we refer to the forthcoming publication [38] and to recent related work on uncertainty quantification and relative entropy, [39].

Block-diagonal 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 k 1 , ..., k Lj ∈ {1, ..., K} 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 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 vectorbecomes a block-diagonal matrix.The pathwise FIM is then written as where A θ 1 , ..., A θ I 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 block-diagonal structure of the FIM is shown on the right panel of Figure 1.Before proceeding with the theoretical computation of the FIM for various well-known 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 block-diagonal 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 complete-data case [40], i.e., obtain Cramer-Rao 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 D-optimality test is given by the relation det(F H ) = I i=1 det(A i ), while the trace of the inverse of the FIM utilized in the A-optimality reduces to tr(F H −1 ) = I i=1 tr(A −1 i ).
(iv) Given that parametric identifiability is characterized by the magnitude of the eigenvalues of the FIM, e.g. a zero eigenvalue corresponds to a non-identifiable direction in parameter space, [18,43], then the block-diagonal structure (11) can provide additional insights in parameter identification.For instance, the identifiability of the parameters of the group that corresponds to the i-th block, A i , will be increased if the smaller eigenvalues of the i-th block can be increased.On the other hand, if the determinant of the i-th block equals to zero then at least one of its eigenvalues is zero and thus the corresponding linear combinations of parameters are non-identifiable.Similarly, the robustness of the system to perturbations of the parameters of the i-th group will be increased if there is a way to decrease the larger eigenvalues of the i-th 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 block-diagonal FIM is demonstrated.

Reactions with independent reaction constants
An important class of well-mixed reaction networks take the general form "α j A j + β j B j θj → ..." 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 j-th reaction.The propensity function for the j-th reaction is given as the product between a rate constant and a function of the current state x: Typically, g j (x) = xA j αj xB j βj 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 explicitly given by where µ θ is the stationary distribution of the stochastic process.Furthermore, it holds that ∂ θ k log a θ j (x) = 1 θ k δ k (j) 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 R K .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 time-series 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 complete-data 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 production-degradation 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, see (25) and (27).
Finally, we would like to remark that even though E µ θ a θ k (x) = θ k E µ θ g θ k (x) 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 production-degradation model (e.g., compare (24) and ( 27)).

Michaelis-Menten kinetics
Another important class of reaction networks is the Michaelis-Menten kinetics.In its simplest form, this chemical network contains two reactions between species A and B (i.e., A ↔ B) with propensity functions given by This reaction sub-network which is derived under a quasi-steady-state assumption is one of the best-known models of enzyme kinetics in biochemistry.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 θ k+1 (usually denoted by K m ) is the substrate concentration at which the reaction rate is half the maximum value.The propensities of this Michaelis-Menten sub-network depend on two parameters (L k = L k+1 = 2 in ( 10)) thus the corresponding FIM block is a 2 × 2 matrix.The elements of the FIM matrix are given by for the k-th row while the k + 1-th row is given by In general, biochemical reaction networks may have significantly more complex propensities, nevertheless, the computation of the FIM elements 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 long-time 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, a θ 0 (x i ), while T = n i=1 ∆t i is the total simulation time.The sequence {x i } n i=0 is the embedded Markov chain with transition probabilities from state x i to state x i+1 is given by the ratio a θ j (xi) a θ 0 (xi) .The weight ∆t i , which is the waiting time at state x i , is necessary for the unbiased estimation of the observable, [44].Similarly, the unbiased estimator for the FIM is Noticing that the computation of the local propensity functions a θ j (x i ) for all j = 1, ..., M is needed for the simulation of the jump Markov process when Monte Carlo methods such as SSA [44] 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.(e) Update FIM:

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 multi-scale approximations of the original SSA have been developed in order to handle such issues resulting to accelerated simulation algorithms.For example, mean-field 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 [45,46].Stochastic corrections to the mean-field model such as stochastic Langevin [25] and linear noise approximation [47] can be applied in order to improve the accuracy of the simulation.An alternative approximation of the jump Markov process is the tau-leap method proposed by Gillespie [24] where a batch of events occurs at each time-increment, τ .Several improvements of the basic tau-leap algorithm on how to select adaptively the τ [48] or avoiding negative populations [27,49] have been proposed, however, their performance is heavily model-dependent.
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 well-known mean-field approximation.The popularity of the mean-field 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 zero-mean 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,47,50].Thus, for large populations, the fluctuations of the time-evolving 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 time-step numerical integrator up to final time T = n i=0 ∆t i .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: Multi-scale 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 [51] or even explicitly available formulas for escape times [52] 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 time-series, 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 time-normalized 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 time dependent relative entropy rate will follow.

Results
A simple protein production/degradation model We first consider an elementary stochastic model for protein production and degradation, [53], which is also a component of more complex models for gene regulatory networks, [54].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, time-dependent moments and autocorrelations, [45,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 k1 k2 .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 supplementary material.Some of the implications of the differences between these two FIMs are discussed next.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, k1 k2 , 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 [45, 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 k1 k2 .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 two-dimensional parameter vector of the system is non-identifiable.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) time-series 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,55,56].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, Mdm2-precursor 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 Mdm2-precursor while y corresponds to Mdm2.The state of the reaction model is defined as x = [x, y 0 , y] T while the parameter vector is defined as Figure 2 shows the time-series 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.
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.Comparison to the LNA-based sensitivity approach In [18], the authors suggested a linear noise approximation (LNA) for the stochastic evolution around the nonlinear mean-field 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 long-time dynamics, e.g.metastable, it is not Figure 2: Molecule concentration of p53, Mdm2-precursor 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.correct as large deviation arguments [51] and explicit formulas for escape times [52] show.Similar issues with non-gaussianity in the long-time dynamics arise in stochastic systems with strongly intermittent (pulse-like) or random oscillatory behavior [57].In the p53 model considered in [18] which had the same parameter values as here, Figure 2 reveals that the time-series 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 3 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 specific striking difference between the two sensitivity approaches is that in our proposed method the sensitivity of parameter b x is relatively high while in the LNA-based method it is lower, see Figure 4 (dark blue) and also compare Figure 3 of this publication and Figure 3 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' time-series.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 time-series, [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 pathwise 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%.It is evident that averaged PSD is more sensitive to perturbations of b x rather than to perturbations of b y as our sensitivity analysis method predicted while the LNA-based method suggested the reverse order of sensitivity.Moreover, the L 1 norm between the unperturbed PSD and 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.An explanation of the performance of the LNA-based 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.

Epidermal Growth Factor Receptor model
The EGFR model is a well-studied system describing signaling phenomena of (mammalian) cells [29][30][31].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 [59,60].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 at an abstract level.Initially, the extracellular binding of EGF with the EGF receptors induce receptor dimerization.Then, two principal pathways, Shc-dependent and Sch-independend, are initiated leading to activation of Ras-GTP.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 1 & 2 of supplementary information in [31].For completeness, all the reactions along with their rates are provided in the supplementary materials of this publication.The propensity functions for the reactions x Bj β j , j = 1, ..., 97, 100, ..., 207 with the exception of reaction pair R 98 , R 99 where their propensity functions are governed by the Michaelis-Menten kinetics a j (x) = V max x Aj /(K m + x Aj ), j = 98, 99 (30) 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 supplementary materials.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].As Figure 8 suggests, 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 gradient-free manner.
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 k-th diagonal element of the FIM corresponds to RER where the perturbation takes place only to the k-th parameter (see (5)). Figure 7 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 multi-parameter models in systems biology and which is known as "sloppiness", [61].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.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 time-dependent 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 time-dependent 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 140-th (k 65 ) most sensitive parameter is perturbed (see Table 1 in supplementary materials).The rate constant k 65 corresponds to a reaction of the Shc-dependent 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 (EGF-EGFR*)2 binding species without Sch* (top, left panel) and with Sch* (top, middle panel) as well as Ras-GTP (top, right panel), total activated Raf or total Raf* (low, left panel), doubly phosphorated MEK or MEK-PP (low, middle panel) and doubly phosphorated ERK or ERK-PP.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.

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 time-series 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 path-space, which in turn constitutes a gradient-free 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 block-diagonal 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 high-dimensional stochastic reaction systems becomes tractable since it is well-known 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 (mean-field) in a high-dimensional 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 quasi-steady 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 high-dimensional 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 mean-field 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 time-dependent observables: practically, our proposed sensitivity analysis represents a conservative sensitivity estimate in the sense that insensitive directions for the relative entropy on path-space, 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.

Figure 1 :
Figure 1: Left panel: The graph representation of the dependencies between the reactions (left column) and the model parameters (right column).The grouping of the parameters is then carried out by noting the connected parts of the graph.Right panel: The corresponding block-diagonal structure of the FIM.In this example K = 7 while the largest dimension of the blocks is L = 3.

Figure 3 :
Figure 3: Upper panel: RER in time for the parameter perturbation of b x (blue), k (green) and a y (red) by +10% (i.e., ǫ 0 = 0.1).The relaxation time of the RER as an observable is very fast.Lower panel: RER for various perturbation directions computed either directly (blue and red bars) or based on FIM (green bars).Direction e k corresponds to perturbation of parameter θ k .

Figure 4 :
Figure 4: The proposed pathwise FIM (left) based on RER as well as the (scaled) FIM based on LNA computed from the StochSens package[58].Evidently, the proposed method uncouples the parameter correlations since most of the off-diagonal elements are zero.

Figure 5 :
Figure 5: Power Spectral Densities (PSD) of the time-series 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 .

Figure 6 :
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.

Figure 7 :
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.

Figure 8 :
Figure8: Time-dependent 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 140-th most sensitive parameter (k 65 ) is ten-fold 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.

Table 2 :
Parameter values for the p53 model.
(upper plot) in conjunction with Table1of the supplementary file EGFR table.pdffullydescribe the (local) sensitivities of the reaction network.Table1in EGFR table.pdfpresents 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 Michaelis-Menten reactions-therefore the diagonal elements correspond to the eigenvalues of the FIM.The

Table 3 :
Initial population of the species for the EGFR network.