Graph-Facilitated Resonant Mode Counting in Stochastic Interaction Networks

TL;DR

This paper introduces a graph-theoretic and real root counting approach to efficiently identify stochastic oscillations and resonant modes in complex networks, simplifying analysis beyond small systems.

Contribution

It presents a novel method combining graph theory and root counting algorithms to determine stochastic resonance properties in large networks.

Findings

Resonant modes depend on the squared Jacobian matrix $J^2$.

Graph-theoretic tools simplify analysis of stochastic behavior in large networks.

Chemical reaction networks with multiple resonant modes can be identified easily.

Abstract

Oscillations in a stochastic dynamical system, whose deterministic counterpart has a stable steady state, are a widely reported phenomenon. Traditional methods of finding parameter regimes for stochastically-driven resonances are, however, cumbersome for any but the smallest networks. In this letter we show by example of the Brusselator how to use real root counting algorithms and graph theoretic tools to efficiently determine the number of resonant modes and parameter ranges for stochastic oscillations. We argue that stochastic resonance is a network property by showing that resonant modes only depend on the squared Jacobian matrix $J^2$ , unlike deterministic oscillations which are determined by $J$. By using graph theoretic tools, analysis of stochastic behaviour for larger networks is simplified and chemical reaction networks with multiple resonant modes can be identified easily.

Equations

$$A$$

$$2X+Y$$

$$B+X$$

$$X$$

$$T(X+1,Y|X,Y)$$

$$T(X-1,Y|X,Y)$$

$$T(X+1,Y-1|X,Y)$$

$$T(X-1,Y+1|X,Y)$$

$$E_X^{\pm }f(X,Y,t)$$

$$E_Y^{\pm }f(X,Y,t)$$

$$\frac{dP(X,Y,t)}{dt}$$

$$+(E_X^{-}-1)T(X+1,Y|X,Y)P(X,Y,t)$$

$$+(E_X^{+}{E}_Y^{-}-1)T(X-1,Y+1|X,Y)P(X,Y,t)$$

$$+(E_X^{-}{E}_Y^{+}-1)T(X+1,Y-1|X,Y)P(X,Y,t).$$

$$\frac{X}{\Omega }=u_1(t)+\frac{x_1}{\sqrt{\Omega }},$$

$${\dot{u}}_1$$

$${\dot{u}}_2$$

$$\dot{\mathbf{x}}=J\mathbf{x}+\mathbf{\lambda },$$

$$J=\left(\begin{array}{cc}B-1& A^2\\ -B& -A^2\end{array}\right),$$

$$P_k({\omega }^2)$$

$$R({\omega }^2)$$

$$Q_k({\omega }^2)$$

$$\frac{\text{d}R({\omega }^2)}{\text{d}({\omega }^2)}=R^{\prime }({\omega }^2)=0$$

$$p_0$$

$$p_1$$

$$p_2$$

$$⋮$$

$$p_i$$

$$⋮$$

$$0$$

$$R({\omega }^2)={\omega }^4+[(B-1{)}^2+A^2(A^2-2B)]{\omega }^2+A^4$$

$$p_0$$

$$p_1$$

$$p_2$$

$$\sigma (0)$$

$$\sigma (\infty )=0.$$

$$B<1+A^2\text{ (steady state condition),}$$

$$(B-1{)}^2+A^2(A^2-2B)<0\text{ (peak condition).}$$

$$a_{n-k}=\sum _{f_k\in \mathcal{G}}(-1{)}^{|f_k|+n+k}\prod _{c\in f_k}A[c]k=1,\cdots ,n$$

$$a_1=(J^2{)}_{v_1{v}_1}+(J^2{)}_{v_2{v}_2}=(B-1{)}^2+A^2(A^2-2B).$$

Full text

*Introduction.—*Interaction networks are ubiquitous in biological physics and mathematics Murray2002 ; Murray2008 ; Kampen2007 , from predator-prey models McKane2005 ; Wang1997 ; Yoshida2003 ; Lugo2008 ; Lotka1925 to the vast field of chemical reaction networks Feinberg1987 ; Feinberg1988 ; Mincheva2007 ; Michaelis1913 . Previous research highlighted how resonant amplification of noise in stochastic interaction networks can lead to behaviour not anticipated from deterministic ordinary differential equation (ODE) models, in particular the emergence of cyclic behaviour in stochastic models where the deterministic counterpart does not show a Hopf bifurcation McKane2005 .

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 Woolley2011 ; Kampen2007 ; Simpson2004 , 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 Kampen2007 . 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 Guldberg1864 ; Feinberg1987 . 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 Feinberg1987 ; Feinberg1988 ; Mincheva2007 . In Feinberg1987 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 multistationarity 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 Mincheva2007 . Graph theoretical methods provide the additional advantage over the approach in Feinberg1987 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 Meng2015 .

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.—*The Brusselator is a model for an autocatalytic reaction such as the Belousov-Zhabotinsky reaction nicolis1977 and follows the reaction scheme

[TABLE]

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 Anderson2010 we can determine the reaction rates

[TABLE]

where $\Omega$ represents the total volume of the system. To formulate the chemical master equation we define the operators $E_{X}^{\pm}$ and $E_{Y}^{\pm}$ on functions of $X$, $Y$ and time $t$ as

[TABLE]

The definition in (2) allows us to write the master equation for the Brusselator in the compact form

[TABLE]

For large $\Omega$ a Van Kampen expansion Kampen2007 of equation (3) is of the form

[TABLE]

with $x_{1}$ as a new stochastic variable. A similar expansion for $Y$ yields the deterministic equations for the Brusselator

[TABLE]

It is well known that the Brusselator exhibits a supercritical Hopf bifurcation when $A^{2}+1=B$ nicolis1977 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 $\Omega$ is large and at the next order we obtain a Fokker-Planck equation Kampen2007 . At steady state it is, however, simpler to use the equivalent representation of a chemical Langevin equation Kampen2007 ; McKane2005

[TABLE]

where bold quantities represent vectors, $J$ is the Jacobian of (5) evaluated at the fixed point,

[TABLE]

and $\bm{\lambda}$ is a vector of Gaussian Markov processes. Equation (6) determines the stochastic behaviour of the Brusselator at large, but finite $\Omega$.

A useful tool to find oscillations in stochastic trajectories is the power spectrum $P_{k}(\omega^{2})=\langle|\hat{x}_{k}|^{2}\rangle$ where $\hat{x}_{k}$ is the Fourier transform of the $k^{\text{th}}$ element of (6) and $\langle\cdot\rangle$ denotes the average over a number of realisations Woolley2011 . The general form of the power spectrum of the $k^{\text{th}}$ species of any interaction network whose stochastic behaviour can be described by equation (6) is

[TABLE]

where $I$ is the identity matrix, $\text{adj}(\cdot)$ is the adjugate matrix, $\text{det}(\cdot)$ is the determinant and $\langle\cdot\rangle$ denotes the average. $R(\omega^{2})$ and $Q_{k}(\omega^{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(\omega^{2})$ reduces to the characteristic polynomial of $J^{2}$ if we let $\omega^{2}=-\lambda$. 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}(\omega^{2})$ or minima of $R(\omega^{2})$ or both. In analogy with the damped harmonic oscillator we define $\omega_{R}$ as a resonance frequency or resonant mode such that $R(\omega_{R}^{2})$ 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 $\langle\lambda_{i}\lambda_{j}\rangle$, 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.—*We now 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(\omega^{2})$ has a minimum which translates into the condition

[TABLE]

and, since the angular frequency $\omega$ 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’ struik2014 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 Sturm1829 . 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)$

[TABLE]

where $\text{rem}(\cdot,\cdot)$ is the remainder of the polynomial long division. Sturm’s theorem proceeds by considering the signs of the Sturm chain $p_{0},p_{1},\cdots,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),\cdots,p_{m}(a)$ and $p_{0}(b),p_{1}(b),\cdots,p_{m}(b)$ is counted which we denote as $\sigma(a)$ and $\sigma(b)$. The number of distinct real roots is simply $\sigma(a)-\sigma(b)$. Letting $a=0$ and $b=\infty$ 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 $\text{det}(A+xI)=x^{2}+\text{Tr}(A)x+\text{det}(A)$, where $Tr(A)$ is the trace. When turning to larger networks, however, Sturm chains become an invaluable tool.

Returning to the Brusselator, $R(\omega^{2})$ is given by

[TABLE]

for which we can build the Sturm chain

[TABLE]

hence,

[TABLE]

Hence, necessary and sufficient conditions for the Brusselator to show stochastic oscillations are

[TABLE]

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(\omega^{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(\omega^{2})$ based on Mincheva2007 .

*Graph theoretic formula for the coefficients of a characteristic polynomial.—*Paper Mincheva2007 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. Maybee2006 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 $\mathcal{G}$. We use a vertex set $V(\mathcal{G})=\{1,\cdots,n\}$ for an $n$ species interaction network. There is an edge from vertex $i$ to vertex $j$ if $J^{2}_{ji}\neq 0$. The convention used in Maybee2006 ; Mincheva2007 is to only draw self loops if $A_{ii}>0$, however, for convenience, we will always draw a self loop if $J^{2}_{ii}\neq 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 $\mathcal{G}$ as a series of distinct vertices $\{v_{i_{1}},\cdots,v_{i_{k}}\}$ connected by edges $v_{i_{1}}v_{i_{2}},v_{i_{2}}v_{i_{3}},\cdots,v_{i_{k}}v_{i_{1}}$. For a cycle $c$ we denote $J^{2}[c]=(J^{2})_{v_{i_{2}}v_{i_{1}}}(J^{2})_{v_{i_{3}}v_{i_{2}}}\cdots(J^{2})_{v_{i_{1}}v_{i_{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 $\mathcal{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}^{\prime}=\{\{v_{1}\},\{v_{2}\}\}$ and two factors of degree one which are identical to the cycles of length one.

Consider the characteristic polynomial $p(x)=\sum_{i=0}^{n}a_{i}x^{i}$ of a matrix $A$. We can now apply a graph theoretic formula for the coefficients $a_{i}$, derived in Maybee2006 and applied to interaction networks in Mincheva2007 ,

[TABLE]

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

[TABLE]

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(\omega^{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 Gillespie1976 , 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 stochastic 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.

We would like to thank Eamonn Gaffney for helpful discussions.