Graph-facilitated resonant mode counting in stochastic interaction networks

Oscillations in dynamical systems are widely reported in multiple branches of applied mathematics. Critically, even a non-oscillatory deterministic system can produce cyclic trajectories when it is in a low copy number, stochastic regime. Common methods of finding parameter ranges for stochastically driven resonances, such as direct calculation, are cumbersome for any but the smallest networks. In this paper, we provide a systematic framework to efficiently determine the number of resonant modes and parameter ranges for stochastic oscillations relying on real root counting algorithms and graph theoretic methods. We argue that stochastic resonance is a network property by showing that resonant modes only depend on the squared Jacobian matrix J2, unlike deterministic oscillations which are determined by J. By using graph theoretic tools, analysis of stochastic behaviour for larger interaction networks is simplified and stochastic dynamical systems with multiple resonant modes can be identified easily.

The main tools for investigating stochastic cycles are based on the calculation of the exact power spectra for the constituents of the network from a Langevin equation [3,13,14], which demands knowledge of noise covariances.The determination of noise covariance requires extensive coarse graining, starting from a master equation formulation of the interaction system, and via weak noise expansions the deterministic equations, and a Fokker-Planck equation can be calculated.Eventually coarse graining allows the use of the simpler chemical Langevin equation [3].We seek to streamline the coarse graining process by showing how the desired information, namely the number of resonant frequencies of a network, can be extracted from the deterministic equations only.We also find the parameter ranges associated with a number of resonant modes using graph theoretical approaches developed for chemical reaction networks.
There is a large body of algebraic and graph theoretic techniques for studying deterministic mathematical models.Usually these interaction networks have a large number of parameters, typically one rate constant per interaction and the model parameters are responsible for the dynamics of the system [9,15].Past research focussed successfully on exploiting the network structure of an interaction system for determining its dynamical behaviour, as network structure is a feature of a model and unaffected by the choice of rate constants [9][10][11].
In [9] it was shown how network structure can be used to determine whether a given chemical reaction network has stable steady states, a useful tool to rule out mul-tistationarity in a network.More recently graph theoretical methods have been employed to show how network features such as feedback cycles can lead to oscillations and multistationarity in chemical reaction networks [11].Graph theoretical methods provide the additional advantage over the approach in [9] that they allow one to explore the bifurcation structure of the network.Despite the apparent advantage of using graph theoretical methods for the investigation of dynamical capabilities of interaction networks the graph based investigation of stochastic models is still in its infancy [16].
In this letter we provide an alternative route for calculating the resonant frequencies (and its parameter range) of stochastically-driven oscillating systems.Instead of solving the roots of a rational function of the power spectrum from the weak noise approximation, we investigate the maxima of this function.To do this, we adapt algebraic techniques (e.g.Sturm chains) and a graph theoretic formulation for finding the coefficients of the characteristic polynomial and thereby offering a methodology for studying stochastically-driven oscillations without requiring excessive expansions.
Weak noise and power spectrum.-TheBrusselator is a model for an autocatalytic reaction such as the Belousov-Zhabotinsky reaction [17] and follows the reaction scheme In the Brusselator model the chemical species A and B are assumed to be constant and hence represent the model parameters.Using the stochastic version of the law of mass action [18] we can determine the reaction arXiv:1702.08747v1[q-bio.MN] 28 Feb 2017 rates where Ω represents the total volume of the system.To formulate the chemical master equation we define the operators E ± X and E ± Y on functions of X, Y and time t as The definition in (2) allows us to write the master equation for the Brusselator in the compact form For large Ω a Van Kampen expansion [3] of equation ( 3) is of the form with x 1 as a new stochastic variable.A similar expansion for Y yields the deterministic equations for the Brusselator It is well known that the Brusselator exhibits a supercritical Hopf bifurcation when A 2 + 1 = B [17] and that the system has a stable steady state at (A, B/A) when A 2 + 1 < B, which is the focus of our analysis.The deterministic equations represent the leading order of the expansion in the limit where Ω is large and at the next order we obtain a Fokker-Planck equation [3].At steady state it is, however, simpler to use the equivalent representation of a chemical Langevin equation [3,4] where bold quantities represent vectors, J is the Jacobian of (5) evaluated at the fixed point, and λ is a vector of Gaussian Markov processes.Equation (6) determines the stochastic behaviour of the Brusselator at large, but finite Ω.
A useful tool to find oscillations in stochastic trajectories is the power spectrum P k (ω 2 ) = |x k | 2 where xk is the Fourier transform of the k th element of ( 6) and • denotes the average over a number of realisations [13].The general form of the power spectrum of the k th species of any interaction network whose stochastic behaviour can be described by equation ( 6) is with where I is the identity matrix, adj(•) is the adjugate matrix, det(•) is the determinant and • denotes the average.R(ω 2 ) and Q k (ω 2 ) are polynomials of degree n and n − 1 with n being the number of species in the network, in the case of the Brusselator n = 2.Note that R(ω 2 ) reduces to the characteristic polynomial of J 2 if we let ω 2 = −λ.Previous approaches proceeded by analysing all n rational functions (8) to determine the exact shape of the power spectra, and hence prove the existence of maxima.We will show how to determine the number of peaks and their parameter ranges by considering a single polynomial equation.Stochastic oscillations manifest themselves as peaks in the power spectra.From equation ( 8) it becomes apparent that peaks may either arise either from maxima of Q k (ω 2 ) or minima of R(ω 2 ) or both.
In analogy with the damped harmonic oscillator we define ω R as a resonance frequency or resonant mode such that R(ω 2 R ) is a minimum.Our definition implies further that the resonance frequencies are properties of the underlying network structure, represented by J 2 , rather than the individual network constituents.Surprisingly, the number of resonant modes is independent of the noise covariances λ i λ j , even though resonance in interaction networks is a stochastic effect, giving further indication that resonance is a network property.Sturm chains for counting the maxima of power spectrum.-Wenow turn to determine the number of resonant modes in a given network and show how parameter ranges for stochastic oscillations can be calculated in the Brusselator.At resonance the polynomial R(ω 2 ) has a minimum which translates into the condition and, since the angular frequency ω is a real number, we are interested in finding all distinct, real, positive solutions to equation (11).A method to determine an upper bound of such solutions is given by 'Descartes' rule of signs' [19] which states that the maximum number of real, positive roots of a polynomial is given by the number of sign changes of consecutive non-zero coefficients, if the terms of the polynomial are ordered with descending variable exponent.Descartes' rule, however, only gives an upper bound and counts multiple roots as distinct roots.
An exact root counting algorithm is given through the computation of Sturm sequences and the use of Sturm's theorem [20].For a univariate polynomial p(x) Sturm's theorem gives the number of distinct real roots in an interval (a, b] with a < b.To apply Sturm's theorem we compute a Sturm chain for p(x) p 0 = p(x), p 1 = dp(x) dx = p (x), p 2 = −rem(p 0 , p 1 ), . . .
for which we can build the Sturm chain hence, and FIG. 1.The phase diagram for the stochastic Brusselator.In the stochastic weak noise regime there exists a band (blue, colour online) between the stable oscillations and the stable steady state where stochastic oscillations can be seen.
Hence, necessary and sufficient conditions for the Brusselator to show stochastic oscillations are This system of inequalities can be solved in a computational mathematics software such as Mathematica to give the region of stochastic oscillations shown in Figure 1.
From equations ( 14) and ( 15) it becomes apparent that often we only need to evaluate specific coefficients of R(ω 2 ) rather than find the polynomial itself.Often, unless exact parameter ranges are needed, even fewer polynomial coefficients need to be considered due to some coefficients' inability to change sign, a feature easily identified from network motives in the graph of J 2 .In the remainder of this letter we will outline a graph-based method to facilitate the finding of coefficients of R(ω 2 ) based on [11].Graph theoretic formula for the coefficients of a characteristic polynomial.-Paper[11] gives a graph theoretic formula for the coefficients of characteristic polynomial of the Jacobian matrix of a chemical reaction network, an application of the earlier work of Maybee et al. [21] who consider a general square matrix A. Following their approach we use the squared Jacobian J 2 as an adjacency matrix for a directed graph G.We use a vertex set V (G) = {1, • • • , n} for an n species interaction network.There is an edge from vertex i to vertex j if J 2 ji = 0.The convention used in [11,21] is to only draw self loops if A ii > 0, however, for convenience, we will always draw a self loop if J 2 ii = 0. Using these conventions we can draw the directed graph for the Brusselator as shown in Figure 2. We define a cycle c of length k in G as a series of distinct vertices {v i1 , Consider the characteristic polynomial p(x) = n i=0 a i x i of a matrix A. We can now apply a graph theoretic formula for the coefficients a i , derived in [21] and applied to interaction networks in [11], where in our example A = J 2 and all other quantities are as previously defined.Therefore, using the cycles and factors we identified in the Brusselator, we can compute the a 1 coefficient By computing a Sturm chain from the generic polynomial p(x) = 2x + a 1 we find that a 1 < 0 for the existence of an extremum of R(ω 2 ), thus, we re-derived condition (15).We simulated the trajectory of the stochastic Brusselator in the parameter regime which satisfies condition (15) using Gillespie's direct method [22], Figure 3, and plotted the power spectrum averaged over 500 repetitions.Our results can be found in Figure 4 and show good agreement with our theoretical prediction.
Conclusions.-Resonance in stochastic interaction networks is a well reported phenomenon and a prominent example of how internal stochasticity can lead to oscillatory behaviour.A vital tool to investigate stochastic oscillations is the power spectrum which is traditionally calculated from the Langevin equation.Current methods, however, require detailed knowledge of the underlying stochastic process which can be troublesome to calculate.In this letter we showed how resonance can be understood as a network property, independent of the noise correlations involved.We used Sturm chains to count the number of resonant modes and outlined a graph based method to determine parameter ranges in which stochas-tic oscillations occur.Future work will seek to extend the application of graph based methods to stochastic spatial systems such as stochastic Turing patterns in interaction networks.
where rem(•, •) is the remainder of the polynomial long division.Sturm's theorem proceeds by considering the signs of the Sturm chain p 0 , p 1 , • • • , p m evaluated at the points a and b.Similarly to Descartes' rule the number of sign changes of p 0 (a), p 1 (a), • • • , p m (a) and p 0 (b), p 1 (b), • • • , p m (b) is counted which we denote as σ(a) and σ(b).The number of distinct real roots is simply σ(a) − σ(b).Letting a = 0 and b = ∞ gives the number of all positive, distinct, real roots.For small networks, especially the case n = 2, the number of real roots follows trivially from the quadratic formula and det(A + xI) = x 2 + Tr(A)x + det(A), where T r(A) is the trace.When turning to larger networks, however, Sturm chains become an invaluable tool.Returning to the Brusselator, R(ω 2 ) is given by
vi 1 vi k .The Brusselator graph in Figure 2 has one cycle of length two with vertices c 1 := {v 1 , v 2 } and two cycles of length one given by the self loops on vertices v 1 and v 2 .A factor f k of degree k of G is a collection of pairwise disjoint cycles covering k distinct vertices with |f k | denoting the number of cycles in f k .The Brusselator has two factors of degree two f 2 = {{c 1 }} and f 2 = {{v 1 }, {v 2 }} and two factors of degree one which are identical to the cycles of length one.

FIG. 3 .
FIG. 3. A trajectory of the Brusselator with paramter valuesA = 1, B = 1.2.The smooth green line (colour online) is the solution of the ODE system (5) and the oscillating trajectory is the stochastic trajectory.