WIXLIB

Reduced Linear Noise Approximation for Biochemical Reaction Networks withTime-scale Separation : The stochastic tQSSA Narmada Heratha)Department of Electrical Engineering and Computer Science,Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge,MA 02139, USADomitilla Del Vecchiob)Department of Mechanical Engineering, Massachusetts Institute of Technology,77 Massachusetts Avenue, Cambridge, MA 02139, USA(Dated: 14 February 2018)Biochemical reaction networks often involve reactions that take place on differenttime-scales, giving rise to ‘slow’ and ‘fast’ system variables. This property is widelyused in the analysis of systems to obtain dynamical models with reduced dimensions.In this paper, we consider stochastic dynamics of biochemical reaction networks modeled using the Linear Noise Approximation (LNA). Under time-scale separation conditions, we obtain a reduced-order LNA that approximates both the slow and fastvariables in the system. We mathematically prove that the first and second momentsof this reduced-order model converge to those of the full system as the time-scale separation becomes large. These mathematical results, in particular, provide a rigorousjustification to the accuracy of LNA models derived using the stochastic total quasisteady state approximation (tQSSA). Since, in contrast to the stochastic tQSSA,our reduced-order model also provides approximations for the fast variable stochasticproperties, we term our method the ‘stochastic tQSSA ’. Finally, we demonstratethe application of our approach on two biochemical network motifs found in generegulatory and signal transduction networks.a)nherath@mit.edub)ddv@mit.edu1

I.INTRODUCTIONMany biochemical processes involve reactions that occur on different time-scales. Forexample, in bacterial cells, the binding of transcription factor to DNA takes place on thetime-scale of seconds, while protein production and dilution are on the order of hours1 .Such a separation in time-scales allows the system variables to be separated into slow andfast groups, and this property can be exploited to reduce the complexity of dynamicalmodels. In particular, for deterministic systems, the quasi-steady state approximation givesa reduced-order model for the slow variables, assuming that the fast variables rapidly reacha steady state2,3 . In the mathematical literature, a system of ordinary differential equations(ODEs) with multiple time-scales is represented as slow and fast subsystems by using asmall parameter to capture the separation in time-scales. The mathematical treatment ofsuch systems is given by two main methods: singular perturbation and averaging4,5 . Thesingular perturbation approach, formalized by Tikhonov’s theorem, involves setting 0 inthe system dynamics to obtain an algebraic equation that approximates the fast variables,which is in turn used to derive an approximation for the slow variable dynamics4,6,7 . Inthe averaging method, a reduced-order model for the slow variables is obtained by theelimination of the fast dynamics via integration of the system functions5 .As opposed to deterministic models, employing time-scale separation for model order reduction remains an ongoing area of research for stochastic models of biological systems. Yet,obtaining reduced descriptions of stochastic dynamics is even more critical than for deterministic dynamics in order to increase the speed of simulation and aid analytical studies. Furthermore, accurate reduced-order models are important for precise parameter estimations8 .The most prominent model used to capture the stochasticity in biological systems is theChemical Master Equation (CME), which considers the species counts as a set of discreterandom variables and describes the time evolution of their probability distributions using aset of ordinary differential equations9,10 . There have been several works that obtain reducedorder representations of the CME under time-scale separation conditions11–25 . Among these,a common approach used to reduce the complexity in the simulations of the CME is to approximate the fast variables by their deterministic quasi-steady state expressions19 . However,validity of this method still remains under investigation20,22,25–27 .The theoretical analysis of the CME is challenging due to the large system size and the2

lack of appropriate analytical tools. Therefore, several approximations to the ChemicalMaster Equation have been developed under the assumptions that the system’s volume andthe molecular counts are sufficiently large. One such approximation is the Fokker-Planckequation (FPE) where partial differential equations are used to describe the evolution of theprobability distribution of the species counts9 . Another approximation, equivalent to theFPE, is the chemical Langevin equation (CLE) where the dynamics of species counts aredescribed by stochastic differential equations28 .There have been several recent works that consider the problem of model order reductionfor biochemical reaction networks modeled by the CLE. The work by Contou-Carrere et al.obtains a reduced-order system for the slow variables of the CLE by adiabatic eliminationof the fast variable dynamics29 . They provide a numerical analysis on the error between thefull and reduced-order systems; however, an analytical error quantification is not provided.Furthermore, their work does not provide an approximation for the fast variables of thefull system. In our previous work30–32 , we considered CLE models of biochemical reactionnetworks with linear reaction rates and obtained a reduced-order system for both slowand fast variables, following a similar approach to the deterministic singular perturbationtechnique. It was mathematically demonstrated that the moments of the reduced systemare within an O( )-neighborhood of the moments of the full system.In addition to the above methods, the mathematical literature also offers model order reduction techniques for multi-scale stochastic differential equations via averaging methods5,33 .Recently, these methods have been applied in the analysis of systems modeled by the chemical Langevin equation34 . However, averaging methods require integration of system functions, which may be challenging for systems that are of high dimension or are nonlinear.Moreover, the averaging methods also do not provide an approximation to the fast variable dynamics. In the case of biochemical reaction networks, it is typically important toapproximate the fast variable dynamics, because many species are mixed - that is, their concentrations are given by the combination of slow and fast variable concentrations. Therefore,we often require both slow and fast components in the reduced-order model to correctly approximate such species dynamics.The Linear Noise Approximation (LNA) is another approximate model for the CME,where stochasticity is represented as random fluctuations around a deterministic trajectoryusing stochastic differential equations or partial differential equations10,35 . Recently, model3

order reduction methods for the LNA have been developed using projection operators27,36or singular perturbation analysis37 . However, in these works, the error between the full andreduced-order models are not analytically quantified. The work by Sootla and Anderson38provides an error quantification for model order reduction of LNA developed by Thomal etal.27,36 , but to do so they assume Lipschitz continuity of system functions, which are notLipschitz-continuous for general systems. Furthermore, the above works only provide anapproximation for the slow variable dynamics and do not approximate the fast variables.In this work, we consider biochemical reaction networks modeled using the LNA. We consider systems with separation of time-scales where the dynamics can be represented in thesingular perturbation form with as the singular perturbation parameter. We then obtaina reduced-order model that approximates both slow and fast variables. We mathematicallydemonstrate that first and second moments of the reduced system variables are within anO( )-neighborhood of the first and second moments of the full system variables. Theseresults, in turn, provide a rigorous justification for the accuracy of LNA models derivedusing the stochastic tQSSA in comparison to the standard quasi-steady state approximation (QSSA). Furthermore, in contrast to the stochastic tQSSA, our reduced system alsoprovides approximations for the fast variables stochastic properties. Hence, we term ourmethod the stochastic tQSSA . The application of our approach is then demonstrated ontwo biochemical network motifs found in gene-regulatory networks and signal transductioncascades. Through these examples, we illustrate the practical applications of the reducedorder models and the necessity of both slow and fast variable approximations for analysis.In particular, using the reduced-order model for the gene-regulatory network motif, we further investigate the parameter conditions under which the standard QSSA provides accurateresults in the stochastic setting.This paper is organized as follows. In Section II, we present the LNA model considered inthis paper. In Section III, we introduce the reduced-order system and present our results onthe error quantification between the full and reduced dynamics. In Section V, we illustratethe application of the results with two examples of biochemical network motifs.Notation: E[·] denotes the expected value of a random variable. k · k denotes theEuclidean norm for vectors and k · kF denotes the Frobenius norm for matrices.4

II.THE LINEAR NOISE APPROXIMATION WITH TIME-SCALESEPARATIONA.Linear Noise ApproximationConsider a biochemical reaction network with n species Y1 , . . . , Yn , in a volume Ω, inter-acting through m reactions of the form:kipi1 Y1 . . . pin Yn ri1 Y1 . . . rin Yn ,i 1, . . . , m,where ki denotes the rate constant of reaction i and pil ril is the change in the number ofmolecules of Yl due to the reaction i. Let y(t) [y1 , . . . , yn ]T be the state of the system ata given time t where each component yi represents the molecular count for each species asa discrete random variable. Then, the Chemical Master Equation for this system is of theformm P (y, t) X [ai (y qi , t)P (y qi , t) ai (y, t)P (y, t)], ti 1(1)where ai (y, t) is the microscopic reaction rate proportional to ki with ai (y, t)dt being theprobability that the reaction i will take place in an infinitesimal time step dt. The variableqi ri pi is the stoichiometry vector where pi [pi1 , . . . , pin ]T and ri [ri1 , . . . , ri1n ]T fori 1, . . . , m39.The Linear Noise Approximation (LNA) is an approximation to the CME, where themolecular counts are represented by continuous variables under the assumption that thesystem volume and the molecular counts are sufficiently large. As shown in the work of van Kampen10 , the LNA is derived by taking y Ωv Ωξ in the CME, where Ω is the systemvolume, v is a vector of deterministic variables and ξ is a vector of random variables thatrepresents the stochastic fluctuations. Then, performing a Taylor series expansion about thedeterministic variable Ωv and equating the terms of order Ω1/2 and Ω0 , it is shown that vgives the macroscopic concentrations and the elements of ξ are Gaussian random variableswith the dynamicsv̇ f (v, t),(2)ξ A(v, t)ξ σ(v, t)Γ,(3)5