On a theory of stability for nonlinear stochastic chemical reaction networks

We present elements of a stability theory for small, stochastic, nonlinear chemical reaction networks. Steady state probability distributions are computed with zero-information (ZI) closure, a closure algorithm that solves chemical master equations of small arbitrary nonlinear reactions. Stochastic models can be linearized around the steady state with ZI-closure, and the eigenvalues of the Jacobian matrix can be readily computed. Eigenvalues govern the relaxation of fluctuation autocorrelation functions at steady state. Autocorrelation functions reveal the time scales of phenomena underlying the dynamics of nonlinear reaction networks. In accord with the fluctuation-dissipation theorem, these functions are found to be congruent to response functions to small perturbations. Significant differences are observed in the stability of nonlinear reacting systems between deterministic and stochastic modeling formalisms.


I. INTRODUCTION
Improved understanding of compelling dynamic phenomena, from molecular chaos to relativistic cosmological models, has been possible to large extent, thanks to such mathematical theories as Lyapunov's theory of stability of nonlinear dynamic systems. 1 First developed over 120 yr ago, Lyapunov's elegant theories also propelled the development of the disciplines of systems and control engineering. 2,3n the area of chemical dynamics, Lyapunov's stability theory laid the foundation for understanding nonlinear chemical reaction systems, e.g., the Belousov-Zhabotinsky reaction pattern formation, 4 or the non-isothermal, continuously stirred chemical reactor. 5tochasticity is often a significant, impactful feature in dynamic systems. 6Stochastic processes in physics and chemistry have been studied since the days of Brown, Einstein, and Langevin. 7In the recent past, stochasticity has become a significant focus of studies of biochemical reaction systems. 8Compellingly, understanding the principles that underlie cellular and biomolecular behaviors often requires accounting for the influence of fluctuations inherently present in small biological systems.
Powerful mathematical formalisms have been developed to model and analyze stochastic dynamic systems, including models based on master equations, 9 Langevin equations, 10 and Fokker-Planck equations. 11A prominent place in the theory of stochastic dynamics is occupied by the fluctuation-dissipation theorem whereby a linear response of a given system to an external perturbation is expressed in terms of fluctuation properties of the system in thermal equilibrium. 12,13or chemical reaction systems, Gillespie's stochastic simulation algorithm has sparked the development of a plethora of modeling approaches. 14Yet, despite the significance of stochasticity in nonlinear dynamics, attempts to develop a theory for quantifying the stability of stochastic reacting systems have largely remained unfruitful.In other words, there has been no universally applicable framework for describing the stability of stochastic chemical reaction networks, analogous to Lyapunov's theory for deterministic dynamical systems.
Here, we discuss a theoretical approach for stability analysis of stochastic nonlinear chemical dynamics.The theory is based on the recently developed zero-information (ZI) closure scheme of chemical master equations. 15With ZI-closure, the steady-state probability distribution of small nonlinear chemical reaction networks may be obtained accurately and quickly.We model the master probability with its moments, i.e., the mean, variance, and skewness.Although lower order moments depend on higher order ones for nonlinear reaction systems, we have found that higher order moments contribute little information needed in reconstructing the master probability.The ZI-closure scheme then finds the higher order moments in terms of known lower order moments by maximizing the information entropy of any reaction network.This is not to say that the ZI-closure scheme is the only potential approach to analyzing stochastic stability, or even moment closure.For example, the finite-projection algorithm 16 is a well known method for approximating solutions to the chemical master equation (CME), and numerous closure schemes both analytical 17 and numerical 18 are currently in use.We utilize ZI-closure because it attacks the problem from a moment-based perspective, utilizing the moment equations as a deterministic analogue to an otherwise stochastic simulation.Additionally, among the many closure schemes available, ZIclosure has demonstrated the level of accuracy required to properly pursue the presented stability analysis methods.
Herein, a process is described that resembles Lyapunov's first method for computing eigenvalues directly from the Jacobian matrix of probability moments at the steady state.
In the examples provided, we show that steady-state solutions are stable with the eigenvalues being negative real numbers.This is certainly not unexpected for stochastic models with Markovian dynamics, but to our knowledge, the connection has not been made before between the inherent stability of probability distributions of Markov models and the eigenvalues of the Jacobian matrix of probability distribution moments.We discuss stochastic oscillations, when eigenvalues appear with non-zero imaginary parts, along with the possibility of observing positive eigenvalues in stochastic reacting systems.
Intriguingly, parallels are drawn to the fluctuationdissipation theorem. 12We observe that time correlation functions of probability distribution moments are congruent to transient relaxation dynamics.Transient dynamics are computed when probability distributions of chemical reaction networks are perturbed from their steady state.
In what follows, we describe the ZI-closure method briefly.We present the theory for computing the Jacobian matrix of steady-state probability moments and for computing the Jacobian eigenvalues.We discuss the calculation of correlation functions of probability distribution moment fluctuations.We model and analyze prototypical nonlinear chemical reaction network models with stochastic dynamics.We compute steady-state solutions and analyze their stability.We conclude by stressing the inherent limitations of the proposed methods and by speculating about their significance in the field of stochastic nonlinear dynamics.

A. Zero information closure scheme
A complete model of randomly evolving chemical reactions is one based on the CME, Here, X is the state of the system, an N-dimensional vector with the concentrations of the N reactants and products in the chemical reaction network.P X; t is the probability of being in the state X at time t.T X |X ′ is the transition probability per unit time of reaction events changing the state from X ′ to X.
If the reaction rate laws for the reaction network are elementary, the CME can be used to generate a set of equations of the form 19,20 Here, µ is a vector of probability moments up to order M (length N M ) and µ ′ is a vector of the higher-order moments needed to close the system (moment order M ′ and vector length N ′ M ).A is thus an N M × N M matrix and A ′ is an N M × N ′ M matrix.Solving Eq. ( 2) is equivalent to solving the CME.The presence of higher-order moments (µ ′ ) in the set of equations results in the moment closure challenge.These terms are present if there are second-order reactions or higher; A ′ is then a non-empty matrix.
We recently proposed a moment closure scheme.ZIclosure rests on the simple notion that higher order moments, although certainly non-negligible in arithmetic value, offer little in the way of information necessary to reconstruct the underlying probability distribution. 15We have shown that ZI-closure, based on the maximization of the system's information entropy, offers a closure scheme for the CME of small, nonlinear, stochastic chemical reaction networks.
The entropy of a distribution is defined by Shannon as follows (for simplicity's sake, we continue the theoretical discussion for a single-component system, with N = 1; the derivation of equations can be extended trivially for higherdimensional, multi-component systems): 21 The entropy defined in Eq. ( 3) can be maximized subject to constraints (the assumed known lower-order moments) to define the maximum entropy distribution, 22 Note again that this is for a single component system with M lower-order moments chosen for the closure scheme; the choice of M depends on the complexity of the underlying probability distribution (in examples, we discuss how up to twelve moments may be necessary to capture some bimodal distributions).The parameters α i are Lagrange multipliers.With the now known maximum entropy probability distribution, higher-order moments can be computed as where f µ ′ i (x) is the function for the i-th higher-order moment (for the above example, µ ′ might be x M +1 and f µ ′ (x) = x M +1 ).

B. Steady-state probability distributions, Jacobian matrices, and eigenvalues
The steady-state distribution of a stochastic chemical reaction network can be determined by setting the left side of Eq. (2) to zero.This calculation can be fast for small reaction networks, certainly faster than kinetic Monte Carlo sampling of the probability distribution. 15 linearized Jacobian (J SS ) may be computed when the steady-state distribution is calculated with ZI-closure, and Eq. ( 2) can be approximated around the steady state as Nonlinear stability analysis can then be performed by calculating the eigenvalues and eigenvectors of the Jacobian.Before we discuss how these eigenvalues and eigenvectors dictate the dynamic behavior and the stability of the reaction network around steady states, we briefly present the calculation of J SS .
The Jacobian matrix can be computed with the following equation (from Eq. ( 2)): The higher-order moments are related to the lower-order moments through the Lagrange parameters α in P H (x). Thus, for µ ′ j with j ≥ M + 1, we write and taking derivatives at the steady state, Note that this is an N ′ M × N M matrix.From a computational standpoint, it is easier to obtain the following two matrices at steady state: Consequently, the Jacobian in Eq. ( 7) is simply Herein, we compute the steady-state probability distribution of small chemical reaction networks.We compute the Jacobian matrix, and the eigenvalues/eigenvectors of J SS , which are now trivial to calculate.See the Appendix for a detailed derivation of J SS and the supplementary material 27 for an example.
The significance of eigenvalues and eigenvectors becomes evident with the introduction of time correlation functions and response functions, which we introduce briefly next.

C. Time correlation functions and response functions
Time correlation functions of stochastic processes reveal insightful information about the characteristic time scales of dynamic phenomena.
Considering a chemical reaction network comprised of N components, which has reached a steady state distribution, the time autocorrelation function of a chemical component x i (t) is defined as Cross-correlations C i j between different components x i and x j can be computed similarly.
For linear stochastic systems, it is trivial to show that the entire correlation matrix, C, of all pairwise correlation functions satisfies 23 ∂C ∂t = AC.
Since this is a set of linear equations, eigenvalue analysis provides an exact solution for the autocorrelation function The prefactors A k depend on the eigenvectors of A and the initial conditions, and λ k are the eigenvalues of A.
With ZI-closure, this theoretical treatment is extended for nonlinear stochastic chemical dynamics.We demonstrate that the correlation functions of probability distribution moments can be computed near the steady state as Analogously to linear system cases, using J SS , one can obtain eigenvalues (λ) for nonlinear systems at steady state and determine time correlation functions using the eigenvalues and eigenvectors of J SS , as in Eq. ( 16).A detailed derivation of the correlation functions is provided in the Appendix and an example is provided in the supplementary material. 27urthermore, we demonstrate that the fluctuationdissipation theorem applies to nonlinear stochastic chemical reactions.Chemical reaction networks can be perturbed to new steady states by changing the value of system parameters.The system's probability distribution then relaxes from one steady state to another.We find that these response functions are congruent to steady state correlation functions of probability moment fluctuations, so long as the perturbation results in a linear system response.
The normalized response function for any moment µ i is calculated as In this case, the system starts in steady state 1.The system parameters are then changed and relaxation occurs towards a new steady-state 2. Typically, we change a single parameter value by 5%.In Secs.III-V, the eigenvalues for several simple models across a range of kinetic values are obtained.We evaluate the accuracy of the eigenvalue analysis by comparing the correlation functions calculated via ZI-closure to the correlation functions calculated by kinetic Monte Carlo sampling (Eq.( 16)).We then confirm fluctuation-dissipation congruency relations between steady state correlations and non-steady state linear responses.
To assess the accuracy of the methods, results from ZIclosure are compared to results from kinetic Monte Carlo simulations.The kinetic Monte Carlo results are obtained from 10 7 stochastic simulation trajectories.This ensemble size ensures sampling convergence for probability distributions of small chemical reaction networks. 24The stochastic simulations were conducted using the Hy3S algorithm. 2484101-4 P. Smadbeck

III. RESULTS
We examine four small reaction network models: a simple reversible dimerization system (model 1), a two component Michaelis-Menten system (model 2), 25 the Schlögl model (model 3), 26 and an N-stage linear cycle.Each model serves a particular purpose in showing that the proposed analysis of stochastic chemical reacting systems is accurate and generally applicable.The description, initial conditions, and parameters for each model are provided in Table I.
Proof of concept that ZI-closure and eigenvalue analysis capture nonlinear stochastic chemical dynamics can be demonstrated with the reversible dimerization model.The reversible dimerization system requires no more than 4thorder moment closure in order to accurately construct the probability distribution. 15Consequently, the derivation of the Jacobian is straightforward and can be detailed analytically for demonstration purposes (see the Appendix and the example in the supplementary material). 27The eigenvalues of the Jacobian are trivial to compute using Eq. ( 13), leading to the calculation of autocorrelation functions.
In Figure 1(a), autocorrelation functions for component A, C A A (t), are plotted for varying kinetic constant values.These functions are calculated via eigenvalue analysis (solid black lines) and kinetic Monte Carlo sampling (shapes).The results clearly demonstrate the accuracy of eigenvalue analysis based on ZI-closure.
Not unexpectedly, the decay rate of probability moment correlations decreases with increasing kinetic constant values.The relaxation of autocorrelation functions for this simple system is captured by a single exponential term, whose exponent is the dominant eigenvalue of the Jacobian matrix.The inset to Figure 1 This article is copyrighted as indicated in the article.Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.84.165.204On: Fri, 08 May 2015 14:47:32 are negative real numbers in all studied ranges of kinetic constants.Because the deterministic model is always stable, it is unsurprising that the stochastic steady states of reversible dimerizations are always stable as well.
In Figure 1(a), the moment response curve subject to perturbation is also plotted (x's).In accordance with the fluctuation-dissipation theorem, the non-steady state response function is congruent to the steady state correlation function.In other words, fluctuations at steady state dictate the response of a reaction network to small external perturbations.We find this to be true for any tested parameter combination.
Figure 1(b) focuses on the response of the probability distribution to perturbation for the case of k 1 = 1 (1/molecule s) and k 2 = 1 (1/s).In this example, the system is perturbed away from the first steady state at t = 2, allowed to relax to a second steady state, and then perturbed back at t = 3.The two responses are congruent, if inverted, confirming that the system exhibits a linear response.All perturbation results presented herein exhibit similar behavior.
Analysis of the Michaelis-Menten reaction network demonstrates that the presented theory is applicable to multi-component systems.This reaction network has two independent components, the substrate S and the enzyme E.This example system exhibits symmetry, such that the eigenvalues and eigenvectors computed when changing k 1 and k 2 are identical to those computed when changing k 4 and k 3 , respectively.Because of this, presented results are limited to ranging k 1 and k 2 .
Figure 2(a) shows the autocorrelation function for the substrate S, C SS (t), calculated with eigenvalue analysis (solid black lines).The same autocorrelation function is calculated with kinetic Monte Carlo sampling (shapes).In addition, the response of the mean substrate ⟨S⟩, F S (t), is shown when the system is subjected to a small perturbation (x's).Figure 2(b) shows the corresponding results for the enzyme, E. Again, steady-state correlations, computed by either ZI-closure or kinetic Monte Carlo, and response functions display relaxation congruency.
The inset to Figure 2(b) shows the dominant eigenvalues calculated for a range of k 1 and k 2 parameters.In order to generate the inset plot, nearly 1000 steady states were calculated varying the kinetic constant values, and the eigenvalues of steady-state Jacobians were computed.This detailed description of the stochastic dynamics of a Michaelis-Menten model is produced in a matter of seconds with ZI-closure, in contrast to an equivalent calculation using kinetic Monte Carlo simulations, which requires many CPU hours.Again, the eigenvalues of the Jacobian are real negative numbers for physically relevant kinetic constant values.
The final nonlinear model is the Schlögl model.As in our prior work, 15 the Schlögl model is studied because of its underlying bimodal distribution.Figure 3(a) shows that with slight changes to a single kinetic parameter (here k 4 ), a range of interesting distributions can be sampled.At high kinetic constant values, k 4 > 3.8 (1/s), a single peak is observed at high numbers of molecules of X.There is a transition range from k 4 = 3.8 to 3.2 (1/s), where the system exhibits a bimodal distribution.At low values, k 4 < 3.2, a single peak is observed at low numbers of molecules of X.Because of its complex form, the probability distribution can be captured accurately with at least a 9-order ZI-closure. 15On the X-k 4 plane, the steady-state results for the deterministic model are also plotted as a solid black line.
The autocorrelation functions for the X component, C X X (t), are computed for numerous steady states, varying the value of k 4 and shown in Figure 3 unimodal and bimodal distributions.The system dynamics is faster when the distribution is unimodal than when it is bimodal.This may be attributed to the time necessary for a system to cross the barrier and transition from one probability mode to another.
In all previous examples, the eigenvalues of the Jacobian matrix are negative real numbers, resulting in monotonic decays for the correlation functions.In order to investigate periodic dynamic behaviors, we model a cyclical system, model 4. In this linear cycle, a large number of first-order reactions are linked to form a cycle (with size N and k i = N (1/s) such that the average cycle time is always 1 s).The moment equations for such a system are analytically solvable, and thus, autocorrelation functions can be produced without moment closure (see the Appendix, where the autocorrelation functions are calculated analytically for a linear system).In Figure 5, the power spectral densities 23 are shown for N = 10 (green) and N = 50 (purple).The spectral densities are computed as a Fourier transform of autocorrelation functions.The moment equation system (solid black lines) and kinetic Monte Carlo (circles) results show a good match.Note that the imaginary part of the eigenvalues (marked out as dashed lines) matches up exactly with the frequency of oscillation for the system.Again, this is not now unexpected, but we have not found in the literature a similar analysis of stochastic oscillations based on eigenvalues of probability moment Jacobians.
In particular, for N = 10, the peak corresponds well with the calculated eigenvalue.For N = 50, two eigenvalues are marked which correspond to the first two power-spectral Results are obtained for both kinetic Monte Carlo results (black solid line) and calculations via moment equations (dots).The imaginary part of the dominant eigenvalues (dashed colored lines) is also shown.With N = 10, the peak corresponds well with the calculated eigenvalue.For N = 50, two eigenvalues are marked which correspond to the first two power-spectral density peaks.For N = 50, there are, in fact, 12 unique imaginary eigenvalue parts (frequency modes) calculated with the moment equations.This suggests that higher order modes may be discernible with moment equation analysis, even when they are impossible to determine numerically.The inset shows the imaginary part to the dominant eigenvalue for a range of N values asymptotically approaching 2π as N goes to infinity.
This article is copyrighted as indicated in the article.Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.84.165.204On: Fri, 08 May 2015 14:47:32 density peaks.For N = 50, there are, in fact, 12 unique imaginary eigenvalue parts (frequency modes) calculated with moment equations.This suggests that higher order modes may be discernible with moment equation analysis, even when they are impossible to determine numerically.The inset shows the eigenvalues for the cycle as N is increased from 3 to 60, with N = 10 and 50 highlighted.As expected, as the cycle size increases, the imaginary part of the dominant eigenvalue converges onto one cycle.

IV. DISCUSSION
An important step toward understanding nonlinear stochastic systems is the calculation of steady-state probability distributions.A panoply of tools, first developed for stability analysis of deterministic dynamic systems, can then find use in stochastic dynamic systems.
The steady-state probability distribution is directly obtained with ZI-closure for small reaction networks.This calculation is challenging with kinetic Monte Carlo algorithms, such as Gillespie's stochastic simulation algorithm, not the least because there is no certainty that an ensemble of stochastic trajectories has indeed relaxed to steady-state.
Probability moments are computed up to a finite order with ZI-closure, with higher-order moments computed by maximizing the probability distribution information entropy.To our knowledge, this is the first time that a master equation closure scheme can be used for stability analysis of nonlinear stochastic chemical reactions.
Nonlinear stochastic chemical reaction dynamics can now be linearized around the steady state.With the resulting Jacobian matrix, a stability analysis is possible around steadystate probability distributions.The eigenvalues of the Jacobian provide direct insight into the dynamic behavior and stability of chemical reacting systems and into possible oscillatory behaviors.
With Jacobian eigenvalues and eigenvectors at steady state, the autocorrelation functions for the chemical components can be computed.These correlations capture the entirety of dynamic features of nonlinear reaction networks.Useful information can be obtained about the time scales of stochastic processes at or near steady states.
In the examples shown, the eigenvalues decrease smoothly with increasing kinetic constants.As expected, higher kinetic constant values result in faster reaction kinetics, captured quantitatively by eigenvalues.
For the Schlögl model, the behavior of stochastic nonlinear dynamics is less apparent.The deterministic Schlögl model has three steady states for a range of parameter values, as shown in Figure 3(a).Two of the steady states are stable, as determined using Lyapunov methods, but the third is not.In notable contrast, the stochastic Schlögl model always has a single stable, steady state probability distribution.This probability distribution at steady state may be bimodal but the Schlögl model exhibits ergodicity and visits the two distinct modes with a frequency proportional to their probability.The system time scales reside within each mode and the time scale for crossing through the probability trough is captured by time correlation functions of probability moments at steady state.
This is an important difference between deterministic and stochastic models with implications on the relevance of different modeling formalisms.Chemical reacting systems away from the thermodynamic limit exhibit fundamentally different stability characteristics as compared to their deterministic counterparts at the thermodynamic limit.
We find that all autocorrelation functions determined using eigenvalue analysis match the autocorrelation functions determined with kinetic Monte Carlo simulations.However, kinetic Monte Carlo simulations require simulating up to 10 7 stochastic trajectories in order to obtain equally accurate results.Such simulations require considerable computational resources, even for the small reaction networks studied here.For ZI-closure and the subsequent analysis, only one optimization step is necessary for the information entropy.As long as the size of the discrete state space of reaction network can be enumerated, ZI-closure may outperform kinetic Monte Carlo simulations.For larger reaction networks or for larger numbers of molecules than the ones studied here, ZI-closure will be computationally challenging to implement.ZI-closure may also be computationally challenging to implement in cases of complex probability distributions.Eigenvalue analysis of the Schlögl model using 9th-order ZI-closure, for example, required the accurate calculation of 18th-order moments.While higher-order closure would have been desirable (12thorder closure was used in Ref. 15), numerical challenges made such analysis impractical.In all of these cases when ZI-closure is difficult to implement, model and state space reduction techniques may prove useful.
Intriguingly, steady state autocorrelation functions are found to be congruent to perturbation responses.This means that steady-state fluctuations provide information not only on the stability of a system around the steady state, but also on how the system will respond dynamically to both external forcing functions and parameter variations.In linear response regimes, the fluctuation-dissipation theorem predicts such congruency, but to our knowledge, this has never before been shown for nonlinear stochastic reactions.
We believe that the analysis is, in principle, applicable to any nonlinear stochastic dynamic system for which the master probability equation can be expressed in terms of probability moments.In practice, we will be limited to small reaction networks barring a major theoretical advance in reducing the state space of multicomponent reacting systems.

V. CONCLUSIONS
ZI-closure solves the chemical master equation of small nonlinear reaction networks.The steady state of nonlinear stochastic chemical reaction networks can be calculated with ZI-closure in a single step maximization of the system's information entropy.The alternative, with a significant computational cost, has been to sample the probability with kinetic Monte Carlo simulations.
Nonlinear dynamics can be linearized around the steadystate probability distribution with ZI-closure.Stability analysis tools, such as Lyapunov's first method, can then be employed This article is copyrighted as indicated in the article.Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 134.84.165.204On: Fri, 08 May 2015 14:47:32 to reveal the stability characteristics of stochastic nonlinear reacting systems.Parameter sensitivity analysis is possible and the design of control systems for nonlinear reactions under the influence of noise becomes conceivable.With these theoretical advances, not possible before the discovery of ZI-closure, small chemical and biochemical systems may be rigorously studied.The possibility also opens up for designing stochastic reaction control systems.
Important differences are observed with this analysis between the stability of deterministic and stochastic modeling formalisms.In stochastic systems, steady states may be observed with finite probability, whereas these states may be unstable according to deterministic models and thus, practically impossible to attain.These are important considerations that impact the choice of modeling formalisms used in capturing the dynamic behavior of small reacting systems, for example, ones often observed in biological systems.

ACKNOWLEDGMENTS
This work was supported by grants from the National Institutes of Health (Grant No. GM111358) and by a grant from the National Science Foundation (Grant No. CBET-1412283).We acknowledge computational support from the Minnesota Supercomputing Institute (MSI) and from the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation Grant No. ACI-10535753.Support from the University of Minnesota Digital Technology Center and the University of Minnesota Institute for Engineering in Medicine is also

APPENDIX: THEORY DERIVATION
In this Appendix, we detail the theory used in the development of the steady-state Jacobian (J SS ) and autocorrelation calculation for ZI-closure non-linear analysis.To start, we will describe the form and derivation of the Jacobian matrix J SS that is central to the eigenvalue analysis.The second part of the Appendix will connect the moment equations to the correlation equation dynamics, used to match the data calculated exclusively using the eigenvalues calculated using ZI-closure to stochastic simulation algorithm (SSA) results.We close with a minor note concerning the non-linear correlation functions which take, in our opinion, a slightly unexpected form.

Steady-state Jacobian for ZI-closure
For a non-linear reaction network with polynomial rate laws, the moment equations take the following form: 20 where µ is the lower-order moment vector and µ ′ is the higherorder moment vector.A and A ′ are matrices corresponding to the linear and non-linear parts of the moment dynamics, respectively.What we need for this analysis is the Jacobian, denoted as J SS , near the steady-state distribution (obtained by setting the left side of Eq. (A1) to zero).Through a closure scheme, the unknown higher-order moments are related to the lower-order moments (µ ′ = F ( µ ) ).Although the analysis presented can be applied to arbitrarily complex networks, for illustrative purposes, let us assume a single component (x) system and the lower-order moments are up to order two (µ = x 0 , ⟨x⟩ , x 2 T ) and the only higher-order moment necessary for closure is the third-moment (µ ′ = x 3 ).In this case, Eq. (A1) is more clearly defined as Given three variables (the lower-order moments) and three functions, the Jacobian is defined as and thus, Three unknowns (sensitivities of the higher-order moments with regard to the lower-order moments) need to be calculated at the steady-state distribution.These relations can be found while determining the steady-state distribution through the Lagrange parameters.
In the steady-state distribution, Newton-Raphson optimization what we call an "expanded" Jacobian is determined. 15art of this matrix for the above example is as follows: Note that the Lagrange parameters are denoted as α i here (in prior publications, they were denoted as λ i ).This was done in order to more clearly differentiate between the Lagrange parameters and the eventual calculation of eigenvalues that are traditionally denoted as λ.With minor mathematical manipulation, these matrices determine these unknown factors explicitly, As can be seen with Eq. (A4), these are all of the unknowns necessary to calculate J SS .Ultimately, J SS can be written as More generally, for any reaction network, the steady-state Jacobian can be calculated by the compact matrix formula as follows: This can be applied to any size reaction network for any order moment closure (although with potentially great computationally difficulty and/or error).Non-linear analysis is performed on the matrix J SS .

Correlation equation dynamics
The next step is to show the connection between the moment equations and the correlation equation dynamics.The moment equations for this particular analysis will take on a slightly different form as to how it has been written previously, This is for a linear system (thus, A ′ and µ ′ are absent), but also in this case, the zeroth-moment terms are taken out of the matrix A into a separate vector b.For an Ncomponent system, the moment vector (for up to order-one) is µ = [⟨x 1 ⟩ , ⟨x 2 ⟩ , . . ., ⟨x N ⟩] T .The regression theorem states simply that with these moment equations, the following is also true: The correlation functions are defined as , and thus, the following is also true: Note that the zeroth-order moment terms (b) disappears.
Ultimately that vector is always multiplied by zero and can be eliminated.It is important to note that this is not done in Subsection 1 of the Appendix with the steady-state Jacobian.
The first Lagrange parameter (for the zeroth moment) is important when calculating J SS , but then is eliminated with regard to the correlation dynamics since its corresponding eigenvalue is zero and represents a steady-state which is subtracted out.Now, solely using ZI-closure, the full exact correlation function trajectory can be calculated.Assume that the system is an N-component system, and you are interested in performing analysis on the first random variable (x 1 ) in particular.In this case, you need to obtain an analytical solution to C x 1 , x 1 (t), and thus, the correlation matrix can be reduced to a single vector, where C x 1 = C x 1 , x 1 ,C x 2 , x 1 , . . .,C x N , x 1 T .
Analysis on matrix A can produce a vector of eigenvalues λ = [λ 1 , λ 2 , . . ., λ N ] T and a corresponding matrix of eigenvectors S = ν 1 ν 2 • • • ν N .Spectral resolution dictates that the analytical solution (near the steady-state) for C x 1 , x 1 (t) can be written as where I is an N × N identity matrix.Since we assume that the probability distribution is stationary, the initial condition of the correlation functions will relate to the steady-state covariances, The key relationship between the correlation functions and the eigenvalues is that the correlation functions can be reduced to This function is calculated exclusively with ZI-closure and does not require the use of the SSA at all.In the manuscript main text, we use the autocorrelation functions to compare to SSA results in order to make a quantitative statement about whether the eigenvalues calculated using moment closure are truly representative of the CME behavior.

Non-linear correlation function notes
It has been suggested 23 (and we are attempting to show) that the following equation will hold near the steady-state distribution for non-linear systems: Here, we are concerned with component x 1 in particular.In the non-linear systems, a closure order is chosen, denoted as M, and there are N M lower-order moments.For illustrative purposes, assume that we are dealing with a two component system (like the Michaelis-Menten system) and choose the closure order as two.Thus, the moment vector will be µ = ⟨x 1 ⟩ , ⟨x 2 ⟩ , x 2 1 , ⟨x 1 x 2 ⟩ , x 2 2 T .If we are concerned with component one in particular then the corresponding correlation vector will be The reason this introduces a subtle complication to the non-linear method has to do with the spectral resolution which relies on the initial condition of the correlation vector (C x 1 (0)).For this particular example, the initial condition is as follows: As can be seen, the covariances cover multiple moment orders and result in somewhat unexpected form for the initial conditions.We felt like this subtlety was important to note for reproducibility.

FIG. 4 .
FIG. 4. Schlögl model correlation and response functions.Results are shown for the three extreme points of the chosen kinetic range (k 4 = 3.0 (1/s) as blue circles, k 4 = 3.5 (1/s) as purple circles, and k 4 = 4.0 (1/s) as red circles).The 9th-order ZI-closure results are shown as solid black lines, the kinetic Monte Carlo results are shown as circles and the response functions are shown with x's.The inset shows the dominant eigenvalue determined with 9th-order ZI-closure.Slower dynamic behavior is observed in the center of the kinetic range (the bimodal range for the probability distribution) than at either end of the range (the unimodal range).

FIG. 5 .
FIG. 5. Linear cycle reaction network.Results are for cycle sizes of N = 10 (green) and N = 50 (purple) (all steps having a rate constant of k = N (1/s)).Results are obtained for both kinetic Monte Carlo results (black solid line) and calculations via moment equations (dots).The imaginary part of the dominant eigenvalues (dashed colored lines) is also shown.With N = 10, the peak corresponds well with the calculated eigenvalue.For N = 50, two eigenvalues are marked which correspond to the first two power-spectral density peaks.For N = 50, there are, in fact, 12 unique imaginary eigenvalue parts (frequency modes) calculated with the moment equations.This suggests that higher order modes may be discernible with moment equation analysis, even when they are impossible to determine numerically.The inset shows the imaginary part to the dominant eigenvalue for a range of N values asymptotically approaching 2π as N goes to infinity.
a 11 a 12 a 13 a 21 a 22 a 23 a 31 a 32 a 33