ISAP - MATLAB Package for Sensitivity Analysis of High-Dimensional Stochastic Chemical Networks

Stochastic simulation and modeling play an important role to elucidate the fundamen-tal mechanisms in complex biochemical networks. The parametric sensitivity analysis of reaction networks becomes a powerful mathematical and computational tool, yielding information regarding the robustness and the identiﬁability of model parameters. However, due to overwhelming computational cost, parametric sensitivity analysis is a extremely challenging problem for stochastic models with a high-dimensional parameter space and for which existing approaches are very slow. Here we present an information-theoretic sensitivity analysis in path-space ( ISAP ) MATLAB package that simulates stochastic processes with various algorithms and most importantly implements a gradient-free approach to quantify the parameter sensitivities of stochastic chemical reaction network dynamics using the pathwise Fisher information matrix (PFIM; Pantazis, Katsoulakis, and Vlachos 2013). The sparse, block-diagonal structure of the PFIM makes its computational complexity scale linearly with the number of model parameters. As a result of the gradient-free and the sparse nature of the PFIM, it is highly suitable for the sensitivity analysis of stochastic reaction networks with a very large number of model parameters, which are typical in the modeling and simulation of complex biochemical phenomena. Finally, the PFIM provides a fast sensitivity screening method (Arampatzis, Katsoulakis, and Pantazis 2015) which allows it to be combined with any existing sensitivity analysis software.


Introduction
Complex reaction networks modeled by stochastic population dynamics defined on a graph representing interactions (reactions) between species, are ubiquitous across many scientific disciplines, e.g., chemistry and chemical engineering, systems biology, bioinformatics, epi-

Event Reaction
Propensity function 1  (Goutsias and Jenkinson 2013). Stochastic simulation and modeling provide an intrinsic understanding of the interplay between complexity and robustness of biochemical reaction networks and have attracted an increasing amount of attention from the modeling, computational and statistical communities. We provide a mathematical description of the exact stochastic simulation algorithm and its connections with approximate models of reaction networks in Section 2. Next we use a very simple example to introduce stochastic simulation and parametric sensitivity analysis. The robustness of a system is a property that allows a system to maintain its functions against internal and external perturbations. Due to its bimodality behavior, the Schlögl model (Vellela and Qian 2009) is often considered as a good example that motivates the need for stochastic simulations. The Schlögl model (Vellela and Qian 2009; Katsoulakis 2013) describes a well-mixed chemical reaction network among three chemical species A, B and X interacting through 4 reactions governed by the corresponding propensity functions (Table 1). The concentrations A, B are kept constant and the reaction rates k 1 , . . . , k 4 are the parameters of the model.
Sensitivity analysis provides a mathematical tool to describe the robustness of a system to perturbations by determining the (usually unknown) most/least sensitive directions in the potentially high-dimensional parameter space. The eigenvalues of the proposed pathwise Fisher information matrix (PFIM) are an approximation of the relative entropy rates (RER), which serve as a measure of the system parameter sensitivities. As shown in Figure 1, the most sensitive direction corresponds to the eigenvector with the eigenvalue 396 while the least sensitive direction corresponds to the eigenvector with eigenvalue 30. To visualize the impacts of the perturbation on different sensitive-leveled parameters towards the reaction system, we compared the stationary distributions of the population of chemical species X for the following cases: the unperturbed process, 5% of perturbation to the most sensitive parameter k 1 and 5% of perturbation to least sensitive parameter k 3 . Figure 2 indicates that the reaction network is more robust to perturbation on the less sensitive parameters as the stationary distribution departs much more from the unperturbed distribution when perturbing the most sensitive parameter comparing to the least sensitive parameter. As a result, the parameter sensitivity analysis of complex networks can be considered as an essential mathematical and computational tool to measure the robustness and the identifiability of model parameters ).
In recent years there has been a significant development in sensitivity analysis tools and software packages for low-dimensional stochastic processes. Existing modeling software packages like COPASI (Hoops et al. 2006), PottersWheel (Maiwald and Timmer 2008), SensSB (Rodriguez-Fernandez and Banga 2010) and StochSS (StochSS Org. 2013) analyze parameter sensitivity with deterministic modeling of dynamical systems. Packages COPASI and SensSB perform sensitivity analysis using finite-difference methods. Package PottersWheel introduced mean optimal transformations (Hengl, Kreutz, Timmer, and Maiwald 2007), a non- parametric bootstrap-based algorithm, and the profile likelihood (Raue et al. 2009) method for sensitivity analysis. Package StochSS implements forward sensitivity analysis (Maly and Petzold 1996;Feehery, Tolsma, and Barton 1997) for the deterministic (ODE-based) solver based on the SUNDIAL's CVODES solver (Hindmarsh et al. 2005). Later several stochastic sensitivity analysis packages have been designed for stochastic systems. For example, package SPSens (Sheppard, Rathinam, and Khammash 2013) carries out stochastic parameter sensitivity analysis with several finite-difference estimators including the independent random number method, the common random numbers method (Sheppard, Rathinam, and Khammash 2012), the common reaction path (Rathinam, Sheppard, and Khammash 2010) method and the coupled finite difference (Anderson 2012) method along with the Girsanov likelihood ratio (Plyasunov and Arkin 2007) and the regularized pathwise derivative (Sheppard et al. 2012) method. Another software package, StochSens (Komorowski, Žurauskiene, and Stumpf 2012), is developed for stochastic sensitivity analysis using the linear noise approximation (Komorowski, Costa, Rand, and Stumpf 2011) and Fisher information matrix.
However, for stochastic models with a high-dimensional parameter space as it is typically the case in biochemical reaction networks, the existing sensitivity analysis approaches are often inapplicable due to the overwhelming computational cost since they either rely primarily on finite-difference or employ various approximation in the dynamics that are not suitable for small population sizes. Traditional finite-difference methods break down for high-dimensional systems due to the complexity of the systems, i.e., too many directional derivatives combined with the high variance of each finite difference estimator (McGill, Ogunnaike, and Vlachos 2012 Figure 2: The stationary distributions of the population of chemical species X for the unperturbed process (black solid line), process after perturbed on the most sensitive parameter k 1 (red dotted line) and process after perturbed on the least sensitive parameter k 3 (green dashed line). sensitivity analysis method in  quantifies the loss of information caused by parameter perturbations between time series distributions. The proposed method is realized by employing the relative entropy rate (RER) which provides a measure of the information loss per unit time in path space after introducing an arbitrary perturbation in the parameter vector. An associated pathwise Fisher information matrix (PFIM) which is derived as an approximation of the RER provides a gradient-free approach (i.e., a sensitivity analysis methodology which does not depend on the perturbation) to quantify parameter sensitivities. The block-diagonal structure of the PFIM in reaction networks  indicates its sparsity which in turn reduces the computational complexity of the stochastic analysis using the PFIM for a K-parameter model from O(K 2 ) to O(K). Being a gradient-free sensitivity analysis method, this pathwise approach shows great advantage in computation speed especially when dealing with complex stochastic systems with high-dimensional parameter space. Due to the property that both RER and PFIM rely only on information for local dynamics rather than the equilibrium PDFs, the pathwise sensitivity method is suitable for non-equilibrium steady state systems (i.e., systems for which we do not know explicitly the equilibrium distribution), which is the case in most reaction networks, driven, or externally coupled systems. Furthermore, the pathwise sensitivity method can be also used as a computational efficient way to screen out insensitive parameters. Indeed, based on a new sensitivity bound (11) which incorporates the variance of the quantity of interest and the PFIM an accelerated sensitivity analysis strategy has been proposed (Arampatzis et al. 2015).
Here we present an information-theoretic sensitivity analysis in path-space (ISAP) MATLAB (The MathWorks Inc. 2017) package which simulates deterministic and stochastic reaction networks with three algorithms and primarily focuses on implementing symbolic and numerical procedures to access stochastic sensitivity analysis by calculating both RER and PFIM for chemical kinetics models as proposed in . In particular, this package delivers numerical realizations of time series using an array of options (e.g., mean field approximation/ordinary differential equation, chemical Langevin approximation and exact stochastic simulation algorithm; Gillespie 2007) depending on the desired accuracy, available resources and simulation regimes.
Compared to other existing software for sensitivity analysis, our package can provide efficient sensitivity analysis on complex stochastic systems especially for those with high-dimensional parameter space due to the sparsity and block-diagonal structure of PFIM which can improve scalability by reducing the computation complexity from quadratic, O(K 2 ), down to linear, O(K), with K being the number of model parameters (Arampatzis et al. 2015). Furthermore, our package is compatible to user-input time series data from other sources and can work with any other sensitivity packages in a complementary manner as a sensitivity screening tool that is cheap and efficient to implement as a first step in the sensitivity analysis using any of the existing methods and software (Arampatzis et al. 2015). Several reaction network models including p53, EGFR, Lotka-Volterra and yeast growth model are provided in the software to demonstrate the applicability of the package and more examples can be found in online databases such as the BioModels Database (Le Novère et al. 2006, http://www.ebi.ac.uk/ biomodels-main/publmodels).

Mathematical background
In this section we briefly introduce the mathematical background of stochastic chemical kinetics simulation and stochastic sensitivity analysis.

Stochastic chemical kinetics simulation
In the same mathematical context of continuous time Markov chains we can also consider well-mixed chemical reaction networks. With these models, however, one ignores the spatial locations of individual molecules and it is only the molecular populations for individual species that are needed to define the dynamical model.
We consider a well-mixed chemical reaction network with N chemical species, C = {C 1 , . . . , C N }, which can interact through M specified chemical reaction channels, R = {R 1 , . . . , R M }, governed by K model parameters, θ = {θ 1 , . . . , θ K }. The state of the system at any time t ≥ 0 is denoted by a N -dimensional vector X t = [X t,1 , . . . , X t,N ] where X t,i is the population of chemical species C i at time t. Let ν be the N ×M -dimensional stoichiometry matrix consisting of ν i,j , the stoichiometric coefficient of species C i in reaction R j . For each reaction j, we have a corresponding N -dimensional stoichiometry vector ν j that identifies the numbers of molecules of every species that are produced and consumed in each instance of a reaction of type j. Provided that the reaction network at time t is in state X t = x, a propensity function, a θ j (x), is defined to identify the approximate probability a θ j (x)dt that the jth reaction occurs in the time interval [t, t + dt]. For example in mass action kinetics the propensities have the form, where g is a polynomial. Other types of propensities include the Michaelis-Menten kinetics with propensity function Here we consider its simplest form, where it gives the rate for what is called a single-substrate reaction without intermediate. It identifies the rate at which molecules of species A transform to species B (i.e., the reaction A → B).
Mathematically, {X t } t∈R + is a countable state space continuous time Markov chain (CTMC) with transition rates, a θ j (·), j = 1, . . . , M . The clock of the jumps from a current state x to a new state x = x + ν j is determined by the transition rates with transition probabilities is the total rate. To fully understand the reaction network, the initial population of all chemical species should also be given.
ISAP considers problems modeled in the framework of stochastic chemical kinetics where the reaction networks can be decomposed to coupled systems of single directional chemical reactions of the form: with the corresponding propensity function, a(x), where x's are the population of chemical species.
The following example is a system of 2 chemical species (X 1 , X 2 ) with 3 reactions: with the 3 propensity functions a j (X) where j = 1, 2, 3: which in this example (not necessary in general) follow mass-action kinetics law. According to the coupled biochemical reaction system, the package extracts a stoichiometry matrix which consists of the coefficients of the chemical species in the reactions and propensity functions which determine the reaction rates. Given the stoichiometry matrix and the propensity functions, realization of all reactions taken place in the given time period can be generated according to any simulation algorithms chosen by the user. In our package, there are three different available options: • Stochastic simulation algorithm (SSA; Gillespie 2007): generates the exact, discreteevent stochastic process.
• Chemical Langevin equation (CLE; Gillespie 2007): generates an approximated, continuous time stochastic process governed by stochastic differential equations.
Although the SSA method provides exact numerical realizations on the time evolution of a well-stirred chemically reacting system, it requires a great amount of simulation time especially when multiscale reaction networks are studied. Many multiscale approximations of the original SSA have been developed to accelerate the simulation algorithms. The CLE method approximates the stochastic process by assuming that all reactions fire simultaneously on a fixed time interval while the MFA method can be regarded as the CLE method without the noise term ).
Mean-field models for reaction networks are commonplace due to their computational efficiency. Here the stochastic process is written as where x(t) is the mean of the process, ξ(t) is the stochastic zero-mean part while η is the amplitude of the stochastic term which is proportional to the inverse square root of the reactant populations. As a consequence, the fluctuations of the time-evolving species populations become vanishingly small compared to the deterministic contributions for large populations and the deterministic term whose dynamics are governed by the ordinary differential equation (ODE) system (Equation 4) dominates the process.
The mean field approximation (Equation 4), as well as a higher order stochastic differential equation (SDE) approximation called chemical Langevin can be readily obtained from the stochastic dynamics.

Pathwise sensitivity analysis
The parameter sensitivities are estimated via calculating the relative entropy between time series (or path) distributions  which quantifies the loss of information when assuming a perturbed distribution instead of employing the true distribution. In the stationary regime, pathwise relative entropy increases linearly in time. We refer to  and  for definitions, derivations and properties of the information quantities discussed here. Let Q θ denote the path distribution of the process with parameter vector θ ∈ R K . The relative entropy rate (RER; Pantazis and Katsoulakis 2013), H(Q θ |Q θ+ ), which is the time average of the relative entropy between path distributions, is a time independent measure of sensitivity since it measures the loss of information per unit time in path space due to an arbitrary -perturbation of parameter combinations in any direction in the parameter space. A Taylor series expansion of RER in terms of reveals that RER is locally a quadratic function of the parameter perturbation vector ∈ R K : Reactions Parameters Pathwise Fisher Information Matrix In chemical reaction networks, reactions typically depend only on a small subset of the parameter vector. This parametric dependence of the propensities can be directly reflected on the PFIM. By grouping the reactions into subsets in such a way that each subset contains the minimum number of reactions having common parameters, the PFIM becomes a block-diagonal matrix upon rearrangement of the parameter vector . In this example, the grouping of the 11 parameters is carried out by noting the connected parts of the graph and that the largest dimension of the blocks is 3.
The PFIM , F H (Q θ ), contains up to third order accuracy all the sensitivity information for the path distribution Q θ with parameter vector θ for any perturbation direction . Therefore, the computation of the PFIM is sufficient up to third order for the evaluation of all the local sensitivities of the path distribution around the parameter vector θ. Unlike the traditional FIM which is computed only based on the stationary distribution, the PFIM gathers additional information by taking into consideration the dynamical aspects of the process as well. Most importantly, an explicit formula for the PFIM is given by This explicit formula reveals that for a typical reaction network the PFIM has block-diagonal structure upon rearrangement of the parameter vector. Moreover, PFIM is a steady state average which only depends on the propensity functions a θ j (x). Figure 3 describes the blockdiagonal structure for the dynamics of the crypt compartmental model (Smallbone and Corfe 2014). Concentrating on the PFIM estimation, since the stationary distribution, µ θ , is usually unknown, an unbiased estimator for the PFIM (Pantazis and Katsoulakis 2013) is computed numerically as an ergodic average. RER and the corresponding PFIM become computationally feasible as they are directly computable from the propensity functions.
For biochemical reaction networks where model parameters may differ by orders of magnitude, it is reasonable to carry out sensitivity analysis by adding perturbations which are proportional to the parameter magnitude. By perturbing the logarithm of the model parameters and using the chain rule, we obtain the logarithmically-scaled PFIM , Since the stationary distribution, µ θ , is often unknown, an unbiased estimator for FIM is given as ergodic averages , where ∆t i is an exponential random variable with parameter given by the total rate, a θ 0 ( . Finally, we would like to address how the PFIM can be utilized for fast parameter screening. We refer to Arampatzis et al. (2015) for a detailed discussion. Proceeding, let us denote by G(·) = [G 1 (·), . . . , G L (·)] the vector with L state-dependent observable functions, G l : X → R, l = 1, . . . , L. One typical option for the observable function is the time-average of a function which is defined in a general setting as where the most common observable is the population of the lth species, i.e., the projection of the state vector to the lth direction (g l (x) = x l ). The sensitivity matrix of the observable functions, S ∈ R K×L , is defined by The element, S k,l , is the sensitivity index (SI) of the lth observable with respect to the kth parameter. For the stationary regime and time-averaged observables, a sensitivity bound for the absolute value of the SIs can be derived either from rearranging the Cramer-Rao inequality (Arampatzis et al. 2015) or from variational formulas (Dupuis, Katsoulakis, Pantazis, and Plecháč 2016). Particularly, the sensitivity bound reads where F H (Q θ ) is the PFIM while C l is a constant related only to the lth observable function, G l . Apparently, this inequality can be utilized as a screening tool to discard the most insensitive SIs since diagonal PFIM elements with small values imply low SIs (Arampatzis et al. 2015).
Classical sensitivity analysis approaches scale poorly in the presence of high dimensions and many parameters. Take the mean field approximations (Equation 4) for example, where we consider a system of N ordinary differential equations, when performing a local sensitivity analysis on the model parameters θ ∈ R K , we define sensitivity indices measuring the impact of small parameter perturbations on the model by Then, a new system of ordinary differential equations is derived and augmented with the original system (Equation 12): Therefore, in order to study the local sensitivity of the system of ordinary differential equations we need to solve K × N additional equations, where both K, N 1. As a result, to estimate parameter sensitivity, the classical sensitivity analysis approaches have to solve an adjunct system of (K + 1) × N differential equations (i.e., the ordinary differential equation system (Equation 12) that governs the dynamics of the reaction network augmented by the ordinary differential equations for the derivatives with respect to the parameters (Equation 14)). In comparison, the proposed pathwise sensitivity analysis approach can estimate parameter sensitivity in a more efficient manner as it only needs to solve N differential equations (Equation 12) and calculate the unbiased estimator of PFIM (Equation 8).

Structure of the package
The ISAP package is implemented as a set of MATLAB functions. By taking the user-defined biochemical reaction system, the package extracts a stoichiometry matrix which consists of the coefficients of the chemical species in the reactions and the propensity functions which determine the reaction rates. With the stoichiometry matrix and the propensity functions, realizations of time series in the given time period can be generated according to different simulation algorithms (SSA, CLE and MFA) by users' choice.
This package provides both, simulators for various processes as well as tools for sensitivity, robustness and identifiability of model parameters (see the schematic in Figure 4). Our package using time-series either locally generated or externally provided by another simulator is able to efficiently compute the unbiased estimator of both RER and PFIM . Furthermore, this package offers functions to extract sensitivity coefficients from the spectral analysis of the PFIM and parameter sensitivities can be ranked accordingly. In our package, simultaneous perturbations can be added to any subset of model parameters to test the robustness of the system and the analysis of eigenvalues of the PFIM provides information on the identifiability of model parameters

Create examples
In order to build an example, users need to create three .txt files and store them in the \examples\example_name\input_files folder. The three input files include the following information of the reaction network: • Coupled one direction reaction functions.
• Initial populations of biochemical species.
• Reaction rate constants for all reactions.
The following three ways are possible to create input files: • Several examples (i.e., crypt compartmental model, NSCLC model, EGFR model, Lotka-Volterra equation, Michaelis-Menten kinetics, p53 gene network, Schlögl's model, and yeast model) are provided in this software.
• Users can also download examples from online databases, such as the BioModels (Le Novère et al. 2006). Our software provides a parser function, sbml2txt.m, to convert the downloaded .xml file to our input files using the following steps: 1. Go to BioModels Database (http://www.ebi.ac.uk/biomodels-main/publmodels) and select the desired BioModel ID (for example BIOMD0000000520).
3. Use the sbml2txt.m function to create the three input files.
Note: The current version of the ISAP package is able to handle the SBML file examples with the following features, otherwise users have to create examples manually: -No 'Rules': assignment rule and rate rule.
-Names of the parameter are distinct to each other.
-Only single compartment size is allowed.
-Parameter value must be constant.
• Users can also follow the formats in Section 4.3 to create examples (three input .txt files: reaction functions, initial populations and reaction rate constants) manually.

Relevant functions
• sbml2txt(file, outputfile) Description: Remove the 'listOfModifiers' tag from the .xml file and rewrite if necessary. Load the .xml file into MATLAB and create three input files: reactions mechanism, reaction constants and initial populations.

Inputs:
-'file': Complete path of the location of .xml file.
-'reaction constants': Reaction rate constants for all reactions.

Create basic data structures
In order to simplify the computation, the package takes the three input files and further generates three .txt files as below.
• 'connectivity matrix': associates each species with an index and represents the reaction network as a matrix.
• 'propensity functions': extracts propensity functions from the reaction mechanism.
• 'stoichiometry matrix': consists of the coefficients of the chemical species in the reactions.

Relevant functions
• rxn2csm(fname) Description: Read three files with the reactions mechanism, reaction constants and initial populations to creates the connectivity matrix, stoichiometric matrix and the propensity functions.
• initialization_lin(rname) Description: Initialization for simulation and sensitivity analysis. Reactions, species, reaction rates, connectivity and initial configuration are specified. Linear-scale sensitivity analysis is considered.
-'get_rate.m' file: Compute reaction rate for the given reaction.
-'get_PFIM.m' file: Compute element of PFIM value for the given reaction.
• initialization_log(rname) Description: Initialization for simulation and sensitivity analysis. Reactions, species, reaction rates, connectivity and initial configuration are specified. Logarithmicscale sensitivity analysis is considered.

Outputs:
-'data': Data structure array (same as the linear case in this section).

Stochastic simulation
Our package provides three different ways to conduct stochastic simulation.

Stochastic simulation algorithm
A stochastic simulation algorithm generates an exact, discrete, single event, stochastic process.

Relevant functions
• reaction_ssa(data) Description: Simulation for well-mixed reaction systems with a stochastic simulation algorithm.
-'output.species_ts': The concentrations of each species at each time step.

Chemical Langevin equation
A chemical Langevin equation generates an approximated, continuous, stochastic differential trajectories.

Relevant functions
• reaction_cle(data) Description: Simulation for well-mixed reaction systems with a chemical Langevin equation.
-'output.species_ts': The concentrations of each species at each time step.

Mean field approximation
A mean field approximation/ordinary differential equation generates approximated, continuous, deterministic differential trajectories.

Relevant functions
• reaction_mfa(data) Description: Simulation for well-mixed reaction systems with an ordinary differential equation.
-'output.species_ts': The concentrations of each species at each time step.

Interpolating simulated time series
Different simulation algorithms provide sample points from different sampling time instants. The package can homogenize the time interval into equal partitions by interpolation and re-sampling.

Relevant functions
• sim_rxn(data, dname, method) Description: Simulate the reaction network. Each time series is interpolated at given time instants and stored in .mat files.

Pathwise sensitivity analysis
Package ISAP makes stochastic sensitivity analysis flexible by providing functions to compute the RER only, PFIM only or both.

Compute both RER and PFIM
With an input time series, package ISAP computes both, RER and PFIM, for all parameters assuming some perturbations are added to the original parameter values.

Relevant functions
• compute_PFIM_RER(data, dname) Description: Compute both PFIM and RER from the stored time series.

Compute RER
With an input time series, package ISAP computes the relative entropy rates for all parameters assuming some perturbations are added to the original parameter values.

Relevant functions
• comp_RER(X, data) Description: Compute a vector of RERs with input time series (species population only, no time) and perturbation matrix (each row represents a perturbation vector to all parameters).

Inputs:
-'X': Time series of the populations of chemical species.

Compute PFIM
With an input time series, package ISAP computes the Fisher information matrix for all parameters assuming some perturbations are added to the original parameter values.

Relevant functions
• comp_PFIM(X, data) Description: Compute PFIM with input time series (species population only, no time).

Inputs:
-'X': Time series of the populations of chemical species.

Example problems
Several examples are included in the examples directory: BIO520: Crypt compartmental model (3 species, 7 reactions, 11 parameters).

BIO427:
Model of epidermal growth factor receptor and type 1 insulin-like growth factor pathways in non-small cell lung cancer (NSCLC) (21 species, 22 reactions, 54 parameters).

Crypt compartmental model
It is widely believed that colonic crypts, millions of invaginations in the lining of the colon, are the places where colorectal cancer first appears. The crypt compartmental approach which accounts for populations of stem cells, differentiated cells, and transit cells helps us to model the development and early progression of colorectal cancer. Here we study a reaction network model for the dynamics of the crypt compartmental model (Smallbone and Corfe 2014) which consists of 3 cell species with interactions through 7 reactions with 11 parameters. First we will take the crypt compartmental model, the model of the colon crypt, as an example to explain how to download and transform examples from the BioModels Database and then introduce the test functions in package ISAP. In general, it is easy to modify the test functions to accommodate a user defined reaction network.
7. Run example_methods_BIO520 and example_isap_BIO520 to do simulation and perform sensitivity analysis respectively.
The main components of the test functions are listed as below: Step 1: Set the directory of the reaction network.
>> dname = examples/BIO520/ ; >> fname = [dname input_files/BIO520 ]; Step 2   Step 5: Plot the species of the reaction network. The time series of the reaction network, see Figure 5, can be realized with SSA, MFA or CLE. A comparison of the simulation time with different algorithms is listed in Table 2. The SSA method requires a great amount of simulation time (71.24 seconds) to provide exact numerical realizations of the system. Assuming that all reactions fire simultaneously on a fixed time interval, the CLE method approximates the SSA method and accelerates the simulation time to 3.57 seconds. Being regarded as the CLE method without the noise terms, the MFA method provides the fastest but least accurate results in 0.22 seconds.
For the purposes of pathwise RER and PFIM calculations, we use the simulation data from the CLE method.
>> PFIM = compute_PFIM_RER(data, dname); >> figure; >> plot_PFIM(PFIM, data); The 3D barplot shows a block-diagonal structure of elements of the PFIM with the blocks denoting the reaction subsets with the same parametric dependence . The height of the bars corresponds to the values for the elements of PFIM. By visually inspecting Figure 4, the sensitive parameters are d 1 , b 1 , c 1 , a 1 and d 2 . In particular, d 2 is the most sensitive parameter. The plot demonstrates that perturbations may affect the dynamics to those sensitive model parameters and the rest of the parameters are robust to variations.

Relevant functions
• example_methods_BIO520 Description: Simulate the crypt compartmental model reaction network with SSA, CLE and MFA. Interpolate the simulated time series.
Outputs: -Plot the simulation time series of all species from different algorithms side by side.
• example_isap_BIO520 Description: Simulate the crypt compartmental model reaction network with SSA, CLE or MFA and interpolate the simulated time series. Add 10 percent perturbation on the parameter values and conduct either linear-scale or log-scale parameter sensitivity analysis by computing RER and PFIM.
Outputs: -Plot 3D barplot for the elements of PFIM.
-Plot the RER computed directly and utilizing PFIM side by side with confidence intervals.

NSCLC model
The second example comes from a Systems Biology approach which aims to understand the molecular biology of the epidermal growth factor receptor (EGFR, also known as ErbB1/HER1) and type 1 insulin-like growth factor (IGF1R) pathways in non-small cell lung cancer (NSCLC) (Bianconi, Baldelli, Ludovini, Crinò, Flacco, and Valigi 2012). This reaction network consists of 21 species which can interact through 22 reactions with 54 parameters.
Step 1: Set the directory of the reaction network.
Step 4: Simulate the reaction network and save the species time series.
>> sim_rxn(data, [dname data/ ], MFA ); Solving the derived system of ODEs with MATLAB routine ode15s, we compute the PFIM at the steady state regime. The associated PFIM has dimension 54 × 54. The block diagonal structure of the PFIM (Figures 7 and 8) indicates sparsity which can be used to reduce the computational costs.

Summary of functions
We summarize the implemented functions showing their hierarchical dependencies.

Conclusion
Package ISAP is a novel and extendable MATLAB package for simulation and computationallyefficient sensitivity analysis of complex stochastic reaction networks. The wide range of capabilities and the detailed documentation make package ISAP ideal for helping the modeler to create user-defined models and generate simulation data. By employing the RER, the implemented sensitivity analysis approach quantifies the information loss per unit time along different parameter perturbations between time series distributions. An associated PFIM addresses a gradient-free approach for parametric sensitivity analysis. Compared to other sensitivity software packages, our package provides a significant advantage in reducing the computational cost especially for complex stochastic networks with a high-dimensional parameter space and can be used as a screening tool for insensitive parameters.