Non-normality and non-monotonic dynamics in complex reaction networks
Zachary G. Nicolaou, Takashi Nishikawa, Schuyler B. Nicholson, Jason R. Green, Adilson E. Motter
NNon-normality and non-monotonic dynamics in complex reaction networks
Zachary G. Nicolaou, Takashi Nishikawa,
1, 2
Schuyler B. Nicholson,
3, 4
Jason R. Green,
3, 4, 5 and Adilson E. Motter
1, 2 Department of Physics and Astronomy, Northwestern University, Evanston, IL 60208, USA Northwestern Institute on Complex Systems, Northwestern University, Evanston, IL 60208, USA Department of Chemistry, University of Massachusetts Boston, Boston, MA 02125, USA Center for Quantum and Nonequilibrium Systems, University of Massachusetts Boston, Boston, MA 02125, USA Department of Physics, University of Massachusetts Boston, Boston, Massachusetts 02125, USA
Complex chemical reaction networks, which underlie many industrial and biological processes,often exhibit non-monotonic changes in chemical species concentrations, typically described usingnonlinear models. Such non-monotonic dynamics are in principle possible even in linear models if thematrices defining the models are non-normal, as characterized by a necessarily non-orthogonal set ofeigenvectors. However, the extent to which non-normality is responsible for non-monotonic behaviorremains an open question. Here, using a master equation to model the reaction dynamics, wederive a general condition for observing non-monotonic dynamics of individual species, establishingthat non-normality promotes non-monotonicity but is not a requirement for it. In contrast, weshow that non-normality is a requirement for non-monotonic dynamics to be observed in the R´enyientropy. Using hydrogen combustion as an example application, we demonstrate that non-monotonicdynamics under experimental conditions are supported by a linear chain of connected components,in contrast with the dominance of a single giant component observed in typical random reactionnetworks. The exact linearity of the master equation enables development of rigorous theory andsimulations for dynamical networks of unprecedented size (approaching 10 dynamical variables,even for a network of only 20 reactions and involving less than 100 atoms). Our conclusions areexpected to hold for other combustion processes, and the general theory we develop is applicable toall chemical reaction networks, including biological ones. Matrix non-normality is perhaps best known for its rolein a counter-intuitive form of nonlinear instability [1].Even when a fixed point is linearly stable in a nonlinearsystem described by ordinary differential equations, if thecorresponding Jacobian matrix is non-normal, a smallbut finite perturbation can transiently grow beyond thevalidity of the linear approximation and enter into thenonlinear regime, preventing the perturbation from de-caying to zero. The discovery of this phenomenon hasled to a thorough study of the spectral properties of non-normal matrices in the context of transient dynamics [2];it has also inspired recent works on implications of non-normality for network and spatiotemporal dynamics [3–17]. Given the common perception that linear dynam-ics are fully understood, the possibility of such transientgrowth offers interesting alternative interpretations forbehavior usually attributed to nonlinearity, such as igni-tion dynamics in combustion and temporary activationof biochemical signals.In network systems, however, even at the level of lin-ear dynamics, fundamental questions remain open con-cerning such transient growth—or, more generally, non-monotonic dynamics. How prevalent is non-normality indynamical networks and how often does it lead to non-monotonic dynamics? While non-normality is known tobe widespread among matrices encoding network struc-tures [12], the question is open for matrices representingdynamical interactions, which have direct implications onnon-monotonic dynamics. Since non-monotonicity couldbe observed for one variable but not for others within thesame system, how can we determine from the initial con-ditions whether a given variable will show non-monotonic behavior? Beyond the known tendency of non-normalityto be correlated with local and global directionality ofthe network [5–7, 12, 17], what other connectivity struc-tures have implications for non-normality and/or non-monotonic dynamics?In this article, we address these questions using thechemical master equation [18, 19], which describes thedynamics of a chemical reaction network in terms of atime-evolving probability distribution of the network’sstate. Master equations play an important role in sta-tistical physics [20, 21] and have been used to modelphysical and chemical processes in various contexts (e.g.,Refs. [22, 23]), including processes on networks (e.g.,Refs. [24–26]). The chemical master equation has the ad-vantage of being exactly linear and amenable to rigorousanalysis, while still reflecting the dynamics of species con-centrations, which are most often alternatively modeledby nonlinear reaction rate equations in the limit of largenumber of molecules. This connection between linearand nonlinear models is possible because, when the num-ber of molecules is large, the nonlinear dynamics in thefinite-dimensional space of species concentrations can belifted to linear dynamics in an infinite-dimensional space(which is akin to how the linear, infinite-dimensionalKoopman operator can be used to study nonlinear, butfinite-dimensional, dynamical systems [27]). Thus, whilethe individual interactions between the species concen-tration dynamics are inherently nonlinear, the interac-tions between the probability dynamics of different statesin the chemical master equation are strictly linear. Akey to our approach, particularly for reactions involvinga finite number of molecules, is to represent these linear a r X i v : . [ n li n . AO ] A ug p r o j e t i o n y x (a) Normal dynamics xx x d e r i v a t i v e a t t = t -202-101-1 0 1 0 2 4 0 2 4 0 2 4-1 0 1 -1 0 1 -1 0 1 t t t (b) Non-normal dynamics <0 v u vuv u vu FIG. 1. Non-monotonicity of normal and non-normal dynamics. In each column, the top panel shows trajectories of a two-dimensional linear system (gray curves in the background and one curve highlighted in green), while the bottom one shows theprojection of these trajectories onto a one-dimensional subspace (black solid line in the top panel) parallel to a given vector u (purple arrow) as a function of time. In the top panel, the dashed lines indicate the two eigenvectors of the system, and the redcolor intensity at each point encodes the rate at which the projection of the trajectory starting from that point initially movesaway from zero (i.e., the derivative of the projection at t = 0, multiplied by the sign of the projection). The vector v shown inpurple is orthogonal to a boundary line of the red region (the line on which the derivative of the projection is zero at t = 0).(a) Normal system, ˙ x = − x − y , ˙ y = − x − y , whose eigenvectors are orthogonal. The dynamics are projected onto a linewith a 55 ◦ slope (left column) and a vertical line (right column). The system exhibits weakly non-monotonic trajectories in bothcases. (b) Non-normal system, ˙ x = 3 x − y , ˙ y = 8 x − y , whose eigenvectors are far from being orthogonal. The trajectories areprojected on to the same lines as in (a). In contrast to the normal system in (a), this non-normal system exhibits a much morepronounced transient growth. While the two representative projection angles are used here to illustrate the contrast betweenthe two systems, the entire range of possible angles are shown in an animated version of this figure provided as SupplementalMaterial (click here for the video). The animation shows that the ranges of initial conditions and projection angles leading tonon-monotonicity is wider for the non-normal system than for the normal system. interactions by a directed weighted network of state-to-state transitions.Here, we study what the structure of such a networkcan tell us about the chemical process it encodes. In par-ticular, we focus on the consequences of non-normality,strongly connected components, and reaction irreversibil-ity underlying the local directionality of network links,examining their implications for non-monotonic dynam-ics. As a representative example of real complex chem-ical reactions, we consider hydrogen combustion, forwhich a dataset is available on the experimentally de-termined reactions, species, and rate constants [28–31].In the master-equation representation, the network growsrapidly in size with the number of atoms, reaching tens ofthousands of nodes for fewer than a hundred atoms. Toaddress the computational challenges of constructing thenetwork and solving the master equation, we developeda thresholding technique that substantially reduces thenetwork size without significantly affecting the accuracy of system trajectory calculations. This technique is in-corporated in our open-source toolbox for automaticallyconstructing a network from a given reaction dataset [32].Equipped with this toolbox, we use hydrogen combustionas a representative example to demonstrate that combus-tion networks are indeed non-normal, with some eigen-vectors having small angles between them (and hence farfrom being orthogonal), leading to non-monotonic behav-ior in the number of molecules of intermediate chemicalspecies.More generally, we derive a rigorous, geometrically in-terpretable condition on the initial probability distribu-tion under which the time evolution of a given set ofchemical species is non-monotonic. In particular, thecondition indicates that the dynamical non-monotonicitydepends both on the initial distribution and the se-lected set of species. The condition also shows that non-normality is not strictly required for non-monotonicity.Indeed, for a generic system, be it normal or non-normal, → O H + O → OHH + O → H + OHHO + O → O + OHH O + O → HO + OHH + O → HO H + → HO + O H O + H + O → H O + HO Ar + H + O → Ar + HO H + O → O + OH2H → H + H → H O + → H O + H H + OH → H OH + HO → H O + OH + HO → H + O H + HO → O + H → H + HO H O + H → H O + OHH + OH → H O + H2OH → H O → H O + OHO + OH → H O + O H O + OH → H O + HO H O + OH → H O + HO → H O + O → H O + O HO + OH → H O + O H H O O HO HOOHArH O (a) (b) (c) +H-H -OH+H O -H +O-O +H O+H-2H +O-O -OH+2H O+H-2H -O -OH+H O+H O +H-2H -O +OH+H O-H +HO -O -OH+H O 25001 5000 7500 10000125005000100007500 ×10 i j W ij (sec -1 )20 H
10 O FIG. 2. Network representations of the H /O combustion process. (a) Conventional bipartite graph representation, in whichthe nodes representing reactions (black dots) have incoming links (blue arrows) from the nodes representing the reactant speciesand outgoing links (red arrows) to those representing the product species. (b) Master-equation representation, in which eachnode represents a state defined by the number of molecules of all species and a weighted link the rate of transition betweentwo states. The network shown encompasses the 11 ,
129 states accessible from the initial state (orange node on the left) with20 H , 10 O , 1 OH, and 200 Ar molecules. The green nodes on the right indicate the states with strictly positive probabilityof occupancy in the steady-state distribution of the chemical master equation, Eq. (1). All the other states are indicated byblack dots, and the (reversible) transitions between the nodes are represented by links. The inset shows a zoom-up of the initialstate and its neighboring states, which are labeled with the changes in species composition relative to the initial state. Thelink strengths in the inset are proportional to the transition rates W ij . (c) Transition rate matrix W for Eq. (1). The statesare indexed by the order of discovery in the state enumeration algorithm based on depth-first search, starting from the initialstate i = 1. there is always a projection of the solution space for themaster equation that leads to non-monotonicity for someinitial distribution. Thus, the key question is whether thedynamics of a given set of species are determined by sucha projection. Non-normal systems, however, distinguishthemselves from their normal counterparts in that non-monotonicity is more prevalent and more pronounced (asillustrated in Fig. 1 using a simple two-dimensional sys-tem). In addition, we establish that the relation be-tween non-normality and non-monotonicity can be ex-pressed in terms of the R´enyi entropy [33], in analogywith the known relation between (the violation of) themolecular chaos assumption in the H-theorem and non-monotonicity of the Shannon entropy [34]. Furthermore,capturing the local link directionality and decomposingthe network into connected components, we reveal theglobal directionality of the network in the form of a di-rected acyclic graph (DAG) linking these connected com-ponents. Starting from an initial state most natural forthe combustion process, the system traverses a linearchain of irreversible steps within the DAG structure. Bycomparing with multiple classes of random networks, weconclude that the existence of this linear-chain structure must be attributed to non-random nature of the real com-bustion networks. I. MASTER-EQUATION FORMULATION
Given a set of N s chemical species and a set of N r reactions involving them, we consider the dynamics ofreactions in a mixture of these species. The state of thesystem is defined by the numbers of molecules of indi-vidual species. For a given state i , we denote by m ik the number of molecules of species k in that state. For agiven reaction n , we denote by ξ nk and ζ nk the numbersof molecules of species k that are reactants and prod-ucts, respectively. Under this reaction, state i with m ik molecules of each species k transitions to a state with m ik − ξ nk + ζ nk molecules of each species k . This al-lows us to represent the system as a network of states(nodes) connected by the state-to-state transitions in-duced by reactions (directed links). This representationis entirely different from (but related to) to the conven-tional representation of the species-reactions relation bya directed bipartite graph. For the combustion example network ( N = 347, 85 SCCs)8H network ( N = 1,042, 227 SCCs)10H network ( N = 2,621, 498 SCCs)12H network ( N = 5,799, 979 SCCs)14H network ( N = 11,631, 1,765 SCCs)16H network ( N = 21,717, 3,023 SCCs) -8 -6 -4 -2 t ( s e c ) FIG. 3. Networks of strongly connected components (SCCs), representing the structures of the weighted networks of statetransitions constructed from the 6 H , 8 H , . . . , 16 H initial states (with the numbers of O molecules and Ar atoms varyingproportionally). Each node represents an SCC, a set of states connected bidirectionally by paths of state transitions. The SCCsare identified in the network constructed by enumerating states without the (cid:15) -thresholding and then removing all transitionswhose rates are much smaller than the rates of the corresponding opposite transitions (see Sec. V for further details). The nodesize is proportional to the number of states in the corresponding SCC. Each directed link indicates that a transition is possiblefrom a state in one SCC to a state in the other SCC. The red star symbol indicates the SCC containing only the initial state.The few color-coded nodes connected by red links are the SCCs of the network obtained when applying the (cid:15) -thresholding,with the node color encoding the time at which the maximum is observed for the probability that the system’s state is inthe corresponding SCC. These SCCs are observed to form a directed linear chain in each case (when excluding SCCs withmaximum probability < − ). we describe in detail below, the bipartite representationis shown in Fig. 2(a) and compared to our representationin Fig. 2(b). Assuming that the molecules are well-mixedwithin a fixed volume, the dynamics can be describedprobabilistically by the chemical master equation [19],which determines the time evolution of P i = P i ( t ), theprobability that the system is in state i at time t :d P i d t = (cid:88) j W ij P j , (1)where W ij ≥ j tostate i given by W ij = (cid:88) n κ nj (cid:89) k m jk !( m jk − ξ nk )! ξ nk ! (cid:0) δ ij (cid:48) ( n ) − δ ij (cid:1) , (2) κ nj is a constant characterizing the rate at which reaction n occurs when the system is in state j , the notation δ ij is used for the Kronecker delta, and j (cid:48) ( n ) is the state towhich the system transitions from state j under reaction n . Since the rate of a reaction depends on the numbers ofmolecules of the reactants but not on those of the prod-ucts, W ij involves ξ nk but not ζ nk . Equation (2) impliesthat the transition rates W ij are independent of the prob-abilities P j , making Eq. (1) strictly linear. The systemtrajectory starting from a specific state, which we labelwith i = 1, can be determined by solving Eq. (1) with P (0) = 1 and P i (0) = 0 for all i >
1. The dynamics canthus be regarded as a Markov process over the networkof states driven by transition rates W ij , which can beregarded as the weight of the link from node j to i .As a specific example of chemical reaction dynamics,we consider combustion in a gas composed of hydrogen,oxygen, and a neutral buffer of argon, under a constant-volume, constant-energy condition. The argon buffer isincluded to absorb most of the energy released by thereactions, so that the evolving temperature of the gas re-mains within the acceptable range for the kinetic model(to be described below) while conserving energy. In ad-dition to the reactants and products of the overall com-bustion reaction, 2 H + O −−→ O, various inter-mediate, short-lived, radical species are created duringthis process (H, O, OH, HO , and H O ). From a giveninitial state of this H /O combustion process with spec-ified numbers of molecules of each species, the numberof states that are accessible through a sequence of reac-tions is finite, and we denote this number by N . Thisis because the number of each atomic species in the gasis conserved during each reaction event. Assuming thathydrodynamic effects are negligible (i.e., the gas is well-mixed, and the kinetic energy released by reactions ther-malizes quickly), each accessible state i has well-definedtemperature T i and pressure p i (which can vary with i due to the energy released or absorbed by the reactions).Because the temperature and pressure should not fall un-realistically low, the number of accessible states is furtherlimited. % R edu c ed (b)
20 30 40 50 60O and H atoms65 N O and H atoms N ~ ( a t o m s ) . (a) N ~ ( a t o m s ) . FIG. 4. Growing size and complexity of the H /O combus-tion networks with increasing number of atoms. (a) Log-logplot of network size N against the number of atoms with(solid circles) and without (open circles) the (cid:15) -thresholdingapplied during state enumeration. In each case, N appears togrow faster than the power law indicated by the dashed line.(b) Percentage reduction in N achieved by the (cid:15) -thresholding,as a function of the number of atoms. To map out the network of states, we implemented analgorithm based on a depth-first search to identify allstates accessible from a given initial state through thenetwork of possible transitions. Our implementation [32]uses Cantera [35] and can be applied to any reaction net-work data in the Cantera or ChemKin format [36]. Thenumber of accessible states can be very large even fora gas with a modest number of atoms. To see this, weconsider the initial state with 2 m molecules of H , m molecules of O , a single molecule of OH (included tostart the chain reaction process), and 20 m atoms of Arat 1 atm and 1 ,
000 K. We refer to this as the 2 m H initial state and label it with i = 1 throughout the ar-ticle. For this initial condition, we additionally imposea minimum temperature of 200 K and a minimum pres-sure of 0 .
01 atm for any state to be included. Note thatthe 2 m H initial state is stoichiometrically balanced, sothat all hydrogen in the mixture will be consumed dur-ing combustion. The number of states accessible from the6 H initial state (i.e., m = 3) is just N = 347, and thecorresponding network has a modest complexity (Fig. 3,the first network at the top). We note that this numberis significantly smaller than the number of all possiblestates with the same numbers of H, O, and Ar atoms(842 states, as determined by direct enumeration) dueto the accessibility requirement and the thermodynamicconstraints mentioned above. The number of accessi-ble states quickly grows as we scale up the number ofmolecules proportionally (Fig. 3), reaching N = 21 , initial state, with just 6 m + 2 = 50 atoms(counting only O and H, since Ar serves only as an energybuffer and does not actively participate in the reactions).In fact, the number of states appears to grow with thenumber of atoms slightly faster than a power law withexponent 4 . W ij , we need the con-stants κ nj in Eq. (2), which in this case is given by κ nj = κ n ( T j , p j )( N A V ) o n − , (3)where o n ≡ (cid:80) k ξ nk is the molecularity, N A is Avogadro’snumber, V is the system volume, and κ n ( T j , p j ) are therate constants defining the kinetic model (which for mostreactions has a temperature dependence of the modifiedArrhenius form, but may also involve efficiency coeffi-cients for a third-body reaction and a pressure depen-dence of a falloff reaction).Since the transition rates W ij span multiple orders ofmagnitude, the flow of probability under Eq. (1) is typ-ically limited to a small portion of the network. Thisobservation can be used upfront to reduce the size ofthe network, and thus the computational burden, with-out significantly affecting the dynamics. For this pur-pose, we apply thresholding in our state enumeration al-gorithm; if there is a transition from state j to a newstate i , then state i is added only if the rate of that tran-sition is not negligible compared to the rate of the oppo-site transition, i.e., if W ij > (cid:15)W ji , where (cid:15) is a (small)threshold [38]. Thus, if there are multiple states j fromwhich the system transitions to state i , we add state i ifand only if at least one of these transitions carries non-negligible flow of probability into state i . Once a full listof states is obtained, both W ij and W ji are computedfor each pair of enumerated states i and j . Unless oth-erwise noted, we employ (cid:15) = 10 − in the remainder ofthe paper, which is small enough for accurate computa-tion of system trajectories (see reference comment [39]for details). For example, for the 20 H initial state, this (cid:15) -thresholding reduces the size of the network by morethan 82% ( N = 63 , → , network and the corresponding 11 , W ≡ ( W ij ) ≤ i,j ≤ N are visualized inFig. 2(b) and (c), respectively. The (cid:15) -thresholding con-sistently yields significant reduction in N and in the com-plexity of the network, as seen by comparing the open andfilled circles in Fig. 4(a). Furthermore, the percentage re-duction in N achieved by this procedure grows with thenumber of atoms in the system, as shown in Fig. 4(b).With this reduction technique, we were able to considernetwork sizes as large as N = 74 , initial state with 92 atoms (counting O and H). II. SPECTRAL ANALYSIS
The linearity of Eq. (1) implies that all informationabout the dynamics is encoded in the eigenvalues andeigenvectors of the matrix W . We first note that thestructure of the matrix guarantees the sum over eachcolumn to be zero, i.e., (cid:80) i W ij = 0 for all j . Usingthis and applying the Gershgorin Circle Theorem [40],it follows that zero is an eigenvalue of W (which canbe degenerate), and that all the other eigenvalues have strictly negative real parts. Moreover, it can be shownthat each right eigenvector associated with the zero eigen-value can be normalized so that its components are allnon-negative and sum to unity (see Appendix). Eachsuch normalized eigenvector, or a convex combination ofmultiple such vectors, is thus a steady-state probabilitydistribution for Eq. (1).If the network of states is strongly connected, orequivalently if the matrix W is irreducible, the Perron-Frobenius Theorem [40] implies that the zero eigenvalueis actually non-degenerate and that the components ofits right eigenvector are all strictly positive, leading toa unique steady-state distribution with P i > i . This holds true for the H /O combustion processconsidered here, since all reactions in that process arereversible, making the networks of states undirected andthus strongly connected. Note, however, that many ofthe reactions have tiny reverse reaction rate, which leadsto many states with low probability P i , as we will showbelow.The temporal evolution of the system under Eq. (1)can be decomposed into eigenmodes. Assuming that thenetwork is strongly connected, we denote by χ i the nor-malized right eigenvector associated with the zero eigen-value, i.e., the unique steady-state distribution. For theremaining N − λ n to denote the n th eigenvalue (in an arbitrary order) and χ ni to de-note the i th component of the corresponding right eigen-vector normalized to unit length in 2-norm. Given theinitial probability distribution P i (0), the deviation fromthe steady-state distribution, Q i ( t ) ≡ P i ( t ) − χ i , canbe expressed in terms of the (generally non-orthogonal)projection onto the eigenbasis, Q i ( t ) = (cid:88) n> α n e λ n t χ ni , where α n ≡ (cid:88) i η ni P i (0) , (4)and η ni is the i th component of the left eigenvector as-sociated with λ n . Here, the left eigenvectors are normal-ized so as to satisfy (cid:80) i η ni χ n (cid:48) i = δ nn (cid:48) , where we recallthat δ nn (cid:48) denotes the Kronecker delta. In this eigen-decomposition, all terms decay monotonically in time be-cause λ n < n >
0, implying that P i ( t ) → χ i as t → ∞ , i.e., the system converges to the steady-statedistribution. The contributions of these eigenmodes arequantified by the mode strength coefficients α n , withtheir decay characterized by the timescale parameters τ n ≡ − / Re( λ n ) > /O combustion pro-cess is apparent in the distributions of α n and τ n shownin Fig. 5(a), spanning many orders of magnitude. Onthe one hand, the vast majority of the eigenmodes decayvery fast as τ n is tiny and thus contribute very little tothe overall dynamics. Indeed, the probability P i ( t ) com-puted using only the 150 modes highlighted in blue inFig. 5(a) stays within 3 × − of that computed usingall eigenmodes for all i and all t , while P i ( t ) computedusing only the 1 ,
000 eigenmodes with largest τ n stayswithin 2 × − of that computed using all eigenmodes (a) -5 -10 -15 -6 -4 -2 (sec) -8 π /2 π /8 π /4 π /8 (c) -8 -6 -4 -2 (d) (sec) t magnitude of vector sum (b) H OH O FIG. 5. Spectrum of the 20 H network in Fig. 2(b). (a) Projection α n onto eigenvectors vs. timescale τ n associated withthe corresponding eigenvalues. Highlighted in blue are the 150 eigenvalues with largest α n among the 1 ,
000 eigenvalues withthe largest τ n . Histograms for α n and τ n are shown on the right and the top of the panel, respectively. In each case, thehistogram is shown separately for the blue and gray eigenvalues. (b) Illustration of non-monotonic dynamics due to a pairof nearly parallel eigenvectors χ n and χ n (cid:48) (green arrows). When the eigenvectors have large coefficients ( α n , α n (cid:48) (cid:29)
1) andassociated with very different timescales ( τ n (cid:29) τ n (cid:48) ), the sum of the corresponding terms in Eq. (4) (blue arrows) initially growsbut eventually shrinks to zero. (c) Angle θ nn (cid:48) between the n th and n (cid:48) th eigenvectors among those corresponding to the bluedots in (a), sorted by descending magnitude of their projection coefficients α n and α n (cid:48) . (d) Non-monotonic decay of (cid:80) i Q i forthe network in Fig. 2(b). For reference, we also plot the expected number of molecules of reactants (H and O ) and product(H O) as functions of time, showing the progression of the combustion process. The red shading indicates the interval of theignition event (defined in the text). (code for efficient trajectory computation available in ouronline repository [32]). On the other hand, there areseveral modes with slow timescales and significant α n .These features are shared with sloppy models [41, 42],which have modes with strength and timescale varyingover many orders of magnitude but can describe the pro-cess accurately with only a handful of these modes.However, if the matrix W is non-normal (i.e.,if the matrix normality condition (cid:80) k W ik W jk = (cid:80) k W ki W kj , ∀ i, j is violated), the monotonic decay of in-dividual eigenmodes does not tell the whole story. Fora normal W , the right eigenvectors are all orthogonal toeach other. It then follows from the eigen-decompositionin Eq. (4) that Q i ( t ) decays monotonically in 2-norm,i.e., (cid:80) i Q i → t → ∞ . For a non-normal W ,the eigenvectors are not necessarily orthogonal to eachother and can be nearly parallel (or even completely par-allel, which is equivalent to having a degenerate eigen-vector). Consider a system whose n th and n (cid:48) th eigen-vectors are nearly parallel, and suppose that the coeffi- cients α n and α n (cid:48) have opposite signs with large magni-tudes, and that the timescales τ n and τ n (cid:48) are very dif-ferent. Then, the sum of the corresponding terms inEq. (4) will exhibit highly non-monotonic dynamics, ini-tially growing to a large size before eventually shrink-ing to zero, as illustrated in Fig. 5(b). For the H /O combustion network of Fig. 2, the matrix W is highlynon-normal, with many eigenvectors having small anglesbetween them (Fig. 5(c)), and we indeed observe (cid:80) i Q i to behave non-monotonically (Fig. 5(d), blue curve).The non-monotonic dynamics emerge before the ignitionevent, which starts at t ≈ . × − sec and ends at t ≈ . × − sec (Fig. 5(d), red shading). Here, theignition event is defined as the interval in which between5% and 95% of the total temperature change is observed.The non-monotonicity continues during the conversion ofH and O to H O, until about half of the fuel is con-sumed at t ≈ − sec. III. NON-MONOTONIC ENTROPY DYNAMICS
The non-monotonic decay described above for non-normal W is closely related to the so-called R´enyi en-tropy H ≡ − log( (cid:80) i P i ) [33]. Like the Shannon entropy H ≡ − (cid:80) i P i log( P i ), the R´enyi entropy quantifies thespread of the probability distribution and can thus beregarded as a measure of the uncertainty about the stateof the system. Just as Boltzmann’s H-theorem guaran-tees monotonic evolution of the Shannon entropy when amolecular chaos assumption is satisfied (which would en-tail W ij = W ji in this setting) [34], the R´enyi entropy ismonotonic when the normality condition (cid:80) k W ik W jk = (cid:80) k W ki W kj is satisfied, since the orthogonality of theeigenvectors guarantees monotonic decay [43]. To illus-trate this, we first note that the following formula can bederived: d H d t = − · (cid:80) ij P i S ij P j (cid:80) i P i , (5)where S ≡ ( S ij ) ≤ i,j ≤ N , S ij ≡ ( W ij + W ji ) is the sym-metric part of W . This equation establishes a connec-tion between information theory and dynamical systemstheory: the rate of change of the R´enyi entropy is di-rectly proportional to the Rayleigh quotient of S , whichis of interest in non-normal growth analysis [2, 3]. Thisconnection reflects the fact that the R´enyi entropy is afunction of the 2-norm of the probability distribution P i .In particular, Eq. (5) implies that the minimum possi-ble d H / d t is determined by the largest eigenvalue of S ,which is non-positive and known as the reactivity (up tothe factor of − W is normal, it is known that the re-activity equals the largest eigenvalue of W , which is zero.This implies that d H / d t ≥
0, i.e., the R´enyi entropy canonly increase during the system’s evolution towards thesteady-state distribution. If W is non-normal, the re-activity can be strictly positive, meaning that the R´enyientropy decreases. The leading eigenvector of S gives thedistribution maximizing the rate of decrease.For the H /O combustion network of Fig. 2, the reac-tivity is indeed strictly positive (Fig. 6(a), larger purpledot), indicating the existence of a probability distribu-tion for which the system exhibits the maximum pos-sible d H / d t (Fig. 6(b), smaller purple dots). Such adistribution can be constructed by properly normaliz-ing the corresponding eigenvector of S , which is guar-anteed to be possible by the Peron-Frobenius theorem.The system evolution starting from that distribution in-deed exhibits the predicted rate of entropy decrease,d H / d t = − . × sec − (Fig. 6(c), purple dot),reflecting the quick focusing of the probability initiallyspread over many states onto a small number of states.This distribution, however, is not the only one withd H / d t <
0, as there are many other strictly positiveeigenvalues for S (Fig. 6(a)), from which distributionssharing the same property can be constructed. The de-creasing H (albeit at a much slower rate than the maxi-mum) is also observed during the system evolution start- × -7 = sec × = sec -1 × = sec -1 R é n y i e n t r o p y -6 -8 -4 -2 (c) (sec) -2 -4 -6 (b) +20H +OHsteady-statemax × -2-4-60 10000 (a) × FIG. 6. Non-monotonic entropy dynamics for the 20 H net-work of Fig. 2(b). (a) Spectra of the non-normal matrix W and its symmetric part S . While the leading eigenvalue of W is zero (larger orange dot), the leading eigenvalue of S isstrictly positive (larger purple dot), indicating that the R´enyientropy H can decrease, with the fastest rate | d H / d t | de-termined by that eigenvalue. (b) Steady-state distribution(orange) and the distribution with the maximum | d H / d t | (purple), constructed by normalizing the leading eigenvectorsof W and S , respectively. Also shown is the single-state distri-bution corresponding to the initial condition used to constructthe network (green). (c) Evolution of H starting from thethree distributions in (b). The green trajectory achieves itsmaximum | d H / d t | at t = 9 . × − sec (green dot), whilethe purple trajectory (shifted forward in time to facilitatecomparison) initially exhibit two orders of magnitude fasterdecrease in H (purple dot). The red shading indicates thesame interval of ignition event shown in Fig. 5(d). ing from the 20 H initial state, with several periods ofd H / d t < t = 9 . × − sec, before the ignition event, indi-cated by the red shading in Fig. 6(c). This suggests theexistence of bottlenecks in the network of states, wherenon-normality directs the flow of probability to accumu-late on a small number of states during the pre-ignitiondynamics and sets up the stage for the explosive ignitiondynamics. IV. CONDITION FOR NON-MONOTONICSPECIES DYNAMICS
We now derive an analytical, geometrically inter-pretable condition for individual species to exhibit non- uv uv (b)
Normal Non-normal X (a) t (sec) -8 -6 -4 -2 N = N = N = N = N = N = N = N = ∞ FIG. 7. Non-monotonic chemical species dynamics. (a) Ex-pected fraction of molecules that are radicals in the H /O combustion process, as a function of time for several differ-ent network sizes N . The seven solid curves correspond tothe 6 H , 10 H , 14 H , . . . , 30 H networks, whose sizes N areshown in the plot. The dashed curve represents the continuumlimit N = ∞ . (b) ( N − P i for N = 3. The dashedlines, orthogonal to the vectors u and v , divide the simplexinto quadrants, with their intersection corresponding to thesteady-state distribution χ i (which would be at the center ofthe triangle if W is normal. Trajectories starting in a grayquadrant exhibits non-monotonic dynamics. monotonic dynamics. Given a subset of species X , let X = X ( t ) denote the expected fraction of moleculesthat are of the species in X at time t , relative to itsvalue for the steady-state distribution. This quantitycan be expressed as X ( t ) = (cid:80) i M i Q i , where M i ≡ (cid:80) k ∈X m ik / (cid:80) k m ik is the normalized counts of moleculesof the species in X for state i . Using this formula, thetime evolution of the species can be calculated from thesolutions of Eq. (1). For example, the evolution of theradicals for the H /O combustion process computed thisway is shown for different initial number of molecules(and thus different N ) in Fig. 7(a). As N increases,the trajectory appears to approach the continuum limitdetermined by Cantera. Thus, our results suggest thatthe sharp peaks of radical concentrations associated withignition observed in the continuum limit have origin innon-monotonic dynamics promoted by the non-normalityof the transition rate matrix W in the chemical masterequation.Since the probabilities P i are all strictly positive andsum to unity, the trajectory of the corresponding pointunder Eq. (1) is limited to the set { ( P , . . . , P N ) | P i ≥ , ∀ i and (cid:80) i P i = 1 } , which is an ( N − N = 3, in which case it is a triangle. Eachvertex of the simplex corresponds to the distribution withall the probability concentrated on a single state. Eachface (or edge for N = 3) corresponds to distributionswith the probability spread over multiple states. Theinterior corresponds to distributions with non-zero prob-ability for all states. Trajectories determined by Eq. (1)travel within this simplex and eventually approach thepoint corresponding to the steady-state distribution χ i .If the matrix W is normal, this point would be at the cen-ter of the simplex, since the vector of all ones is a righteigenvector of W (in addition to being a left eigenvector),which implies that the steady-state distribution is uni-form. However, if W is non-normal, the point could beanywhere in the simplex and can be close to the bound-aries of the simplex, since χ i can be non-uniform. In-deed, for the H /O combustion network of Fig. 2, itlies very close to a 24-dimensional hyper-edge, with a 2-norm distance ≈ . × − ( ≈ .
01% of the distancefrom the center of the 11 , .
99% of the probability is concentrated on the corre-sponding 24 states (mostly composed of H O, with justa few molecules of H , O , and other radicals).To derive the condition for non-monotonic dynamics,we first note that X converges to zero as t → ∞ , since X is defined relative to the steady-state distribution. Thus,the condition for X to initially move away from zero be-fore converging to zero is that X and its derivative d X/ d t ,given by d X/ d t = (cid:80) ij M i W ij Q j , have the same sign.To rewrite this condition in a geometrically interpretableform, we define vectors u and v as those with the i thcomponent u i ≡ (cid:80) j M j O ji and v i ≡ (cid:80) jk M j W jk O ki ,respectively, where O ij ≡ δ ij − /N . The vector u canbe interpreted as the projection of the vector m ≡ ( M i )onto the simplex. Likewise, v is the projection of thevector m (cid:48) ≡ ( M (cid:48) k ), M (cid:48) k ≡ (cid:80) j M j W jk , onto the sim-plex. Using these vectors, we can rewrite X and d X/ d t as X = (cid:80) i u i Q i and d X/ d t = (cid:80) i v i Q i . The simplexcan be divided into “quadrants” by two hyperplanes:one orthogonal to u (given by (cid:80) i u i Q i = 0) and theother orthogonal to v (given by (cid:80) i v i Q i = 0). Both hy-perplanes contain the point Q i = 0, which correspondsto the steady-state distribution χ i . Then, a geometriccondition for exhibiting non-monotonic dynamics aftertime t is that the point Q i ( t ) lies in either the quad-rant (cid:8)(cid:80) i u i Q i > (cid:80) i v i Q i > (cid:9) or the quadrant (cid:8)(cid:80) i u i Q i < (cid:80) i v i Q i < (cid:9) . These quadrants areshaded in gray in Fig. 7(b). The vectors u and v , as wellas the quadrants of non-monotonicity, can also be definedand are shown in Fig. 1 for normal and non-normal two-dimensional systems. (We note that u and v for thesesystems do not require the projection by O ij , since theirstates are not constrained to a simplex.)The condition just derived shows that, for both nor-0 = {5H , 2O , 3H, 2HO , H O}= {4H , 3H, 3O, 2HO , H O} = {2H , 4H, 4O, OH, 4H O}= {5H , H, 2O, 5OH, HO } -8 -6 -4 -2 t ( s e c ) FIG. 8. Flow of probability over the DAG structure in the H /O combustion network constructed from the 8 H initial state.In each panel, the network is visualized as in Fig. 3, except that the initial probability distribution is fully concentrated on arandomly chosen state (red star symbol). During system evolution, probability branches out and flows downward along variouspaths (red arrows; showing only those with probabilities < − ), eventually converging to the SCC at the bottom (greencircle), which supports the steady-state distribution. Note that the four randomly chosen initial states are different from the8 H initial state used to construct the network. mal and non-normal W , there are initial distributionsleading to non-monotonic behavior of X . This meansthat, even though the 2-norm (cid:80) i Q i converges to zeromonotonically for normal W (as mentioned earlier), theconvergence of the projection X = (cid:80) i u i Q i may be non-monotonic. We thus conclude that non-normality of W isnot a requirement for non-monotonicity of X . However,non-normal W is different from normal W in that thesteady-state distribution tends to lie near the boundariesof the simplex. This implies that the quadrants of non-monotonicity occupy larger fraction of the simplex thanfor normal W , indicating that non-normal W is morelikely to exhibit non-monotonic dynamics. Moreover, thenon-monotonicity tends to be more pronounced for non-normal W , since the further away the initial point isfrom the hyperplane orthogonal to v , the larger is therate at which X moves away from zero (recalling thatd X/ d t = (cid:80) i v i Q i is the distance from the hyperplane).For our H /O combustion example, we indeed observe asharp transient increase in the expected number of radi-cal molecules, as we saw in Fig. 7(a). V. CONNECTED COMPONENT ANALYSIS
In general, the network of states representing the sys-tem (1) is directed and can be decomposed into stronglyconnected components (SCC), where an SCC is definedas a subset of nodes in which a directed path exists be-tween every node pair in both directions. This allows fora coarse-grained representation of the network by a (dif-ferent) network of SCCs. In this representation, whichwe refer to as the
SCC network , each SCC is regardedas a node, and a directed link is drawn from one SCC toanother if, in the original network, there is a directed linkfrom some node in the first SCC to a node in the secondSCC. It follows from a well-known fact from graph the-ory that the SCC network must be a DAG, i.e., it doesnot contain any closed loop.For the H /O combustion, the whole network is actu-ally strongly connected (since all the reactions are mod-eled as a reversible one) and thus has just one SCC. How-ever, a non-trivial DAG structure emerges when we takeinto account the link weights W ij , many of which arenegligibly small compared to the weights W ji of the re-ciprocal links. This is because the reaction giving rise toa link is often nearly irreversible under the given condi-tion in the sense that its rate in one direction is orders of1magnitude larger than in the other direction. To capturethis link directionality, we remove all transitions fromstate j to state i satisfying W ij < µW ji for a given smallconstant µ >
0. We note that this µ -thresholding pruneslinks representing negligibly rare transitions, while the (cid:15) -thresholding introduced above prunes nodes represent-ing rarely visited states (together with the links involv-ing those nodes). Once the µ -thresholding is applied, wedecompose the resulting network into SCCs, revealinga DAG structure. Throughout this article, we use µ =2 × − . With this value, we find that the µ -thresholdingremoves much fewer links than the (cid:15) -thresholding (4 . .
3% of all links in the 20 H network con-structed with no (cid:15) - or µ -thresholding) and does not signif-icantly affect the dynamics (with maximum error in P i ( t )less than 5 × − for all i and all t ). Because the SCCnetwork has a DAG structure, one can arrange the SCCsin horizontal layers in such a way that the directions ofall links, and thus flows of probability under Eq. (1), aredownwards. This layered SCC arrangement is used inFig. 3 for the 6 H , 8 H , . . . , 16 H networks constructedwithout the (cid:15) -thresholding but with µ -thresholding. TheDAG structure dictates the downward flow of probabilityfrom an arbitrary initial distribution to the final steady-state distribution contained in the bottom SCC, as shownin Fig. 8 for randomly chosen initial state in the 8 H network. Code for computing the SCCs, as well as datafiles describing individual states, is available in our onlinerepository [32].We observe that the initial condition used to constructthe network tends to belong to an SCC located towardthe bottom of the layered DAG structure (the red stars inFig. 3). If we focus now on the SCCs that are accessiblefrom that initial state, we find that these core SCCs forma linear chain leading to the SCC containing the steady-state distribution (color-coded circles in Fig. 3, excludingSCCs with P i < − for all t ). The same linear chainis obtained also if we take the network constructed withthe (cid:15) -thresholding and then remove insignificant reversetransitions using the µ -thresholding. While the expectednumber of the reactants H and O decay monotoni-cally and the product H O increase monotonically alongthe chain, the numbers of radical molecules change non-monotonically, as shown in Fig. 9 for the 20 H network.While this coarse-grained view of the network summa-rizes the dynamics with a simple linear chain, complexnetwork structures exist both within and between the in-dividual SCCs in the chain, as shown in Fig. 10. Thisfull network is the same as that in Fig. 2(b), except thatit is visualized using the linear chain structure in Fig. 9.The DAG structure of the SCC network has a directimplication for the steady-state distribution χ i . Indeed,we can show that P i = χ i is strictly positive only fornodes i belonging to the SCCs without any outgoinglinks, which can always be placed at the bottom of thelayered arrangement (see Appendix for a rigorous proof).Since the number and the sizes of such SCCs are oftensmall, this implies that the steady-state distribution is -8 -6 -4 -2 t (sec)0 e x p e c t e d nu m b e r o f m o l e c u l e s c o r e S CC s H OH O HOH O H O HO FIG. 9. Molecular composition profiles in the linear chain ofcore SCCs in the 20 H network in Fig. 2(b). For a given SCC,we show the expected number of each chemical species in astate belonging to that SCC at the time of maximum prob-ability (encoded as the node color). We observe that, whilethe numbers of reactant molecules (H and O ) monotonicallydecrease and the number of product molecules (H O) mono-tonically increases along the linear chain, some of the radicals(H and O) have peaks for intermediate SCCs. localized. For the networks in Fig. 3, the distribution ishighly localized, with all probability concentrated in thegreen SCC at the bottom of each network (and furtherlocalized within that SCC, as illustrated in Fig. 10 for the20 H network). Moreover, we observe a similar localiza-tion even when the µ -thresholding is not applied and thewhole network forms a single SCC, as we saw in the pre-vious section, reflecting the fact that the dynamics areessentially unaltered by the thresholding. VI. RANDOM NETWORKS
To interpret the linear chain structure shown in Figs. 9and 10 for the ( µ -thresholded) 20 H network, we nowcompare its topological properties to those of generic di-rected random networks. We first consider fully randomnetworks with the same numbers of nodes N = 11 , M = 150 , N nodes with very high probability (96 . ,
000 realizations), and the average sizeof the giant SCC is ≈ , . ± .
18. The giant SCCappears because the average degree d ≡ M/N ≈ . d = 1 [44, 45].In 1 ,
000 realizations, any node not in the giant SCC ei-ther had outgoing links only (0 . ± .
12 nodes on aver-age) or incoming links only (0 . ± .
13 nodes). Thus,the linear-chain structure of the SCC network (the color-2 -8 -6 -4 -2 t ( s e c ) FIG. 10. Detailed connectivity structure within and between the core SCCs in the linear chain from Fig. 9 (reproduced hereat the bottom). The size of each node i is proportional to the maximum of P i ( t ) over t , and the color coding for the nodes isthe same as in Fig. 9. Nodes with negligibly small probability ( P i ( t ) < − for all t ) are excluded to avoid overcrowding. coded SCCs in Fig. 9) is not typical of a network withthe same numbers of nodes and links. One possible rea-son for this is that the in- and out-degree distributionsof the combustion network are different from the Pois-son distributions expected for the random counterpart(Fig. 11(a)). However, even if we consider random net-works with the same expected in- and out-degree dis-tributions [46] as the combustion network, the networkstill typically consists almost entirely of the giant SCC(having 11 , . ± .
01 nodes on average, when esti-mated from 100 realizations). Moreover, the probabil-ity that a directed link is reciprocated is relatively high( ≈ .
36) in the combustion network, and this is far fromtypical for the random networks. Indeed, for either ofthe two random models, the probability that a directedlink is reciprocated is p / [2 p (1 − p )] ≈ . p = M/ [ N ( N − /O combustion process), we now consider the class of randomchemical reaction networks defined for a given number N r of randomly chosen reactions involving a given set of N s species. For simplicity, we consider only bi-species,bi-molecular reactions of the form S + S → S + S in which the four species involved are distinct and cho-sen randomly from the N s species. For each reactionchosen, we include the reverse reaction with probabil-ity R (which can be tuned to reproduce the probabilityof reverse transition observed for the H /O combustionnetwork). For a given set of reactions and given the to-tal number of molecules N m (which is conserved by anybi-species, bi-molecular reaction), we construct the corre-sponding matrix W through Eq. (2), which in this case reduces to W ij = κ n ( T j ,p j ) N A V m jk m jk for i (cid:54) = j , wherereaction n involves two reactants S k and S k and trans-forms state j to state i . Note that the factor κ n ( T j ,p j ) N A V can be omitted for this analysis as it does not affectthe topological properties of the network. For N s = 8, N r = 33, N m = 9, and R = 0 .
39 (leading to the prob-ability of reverse transition ≈ . N = 11 , M = 156 , . ± , .
33 (yieldingthe average degree of 13 . ± . , . ± .
94, representing ≈ .
3% of the nodes(Fig. 11(b)). Similarly to the purely random networksconsidered above, this class of random networks also un-dergoes a percolation transition close to average degree= 1 (Fig. 11(c), blue curve). The number of nodes withoutgoing (14 . ± .
93) or incoming (33 . ± .
12) linksonly constitutes a tiny fraction of the network (thoughsignificantly larger than for fully random networks). Inparticular, the nodes with only incoming links, which in-dividually form single-node SCCs, corresponds to a tinyfraction of the states on which the steady-state probabil-ity concentrates. This fraction becomes even smaller asthe average degree increases (Fig. 11(c), orange curve).
VII. CONCLUSIONS
The theory developed here reveals two distinct mech-anisms through which non-monotonic dynamics emergefrom non-normality of the network. First, the non-normality moves the reactivity (the largest eigenvalue of S ) away from zero, enabling the decrease of the R´enyientropy and inducing non-monotonic entropy dynamics.Second, non-normality allows the steady-state distribu-3 in-degreeout-degreePoisson (a) (b)degree p r o b a b ili t y mean degree20 H network(c) r e l a t i v e s i z e giant SCC P i > 0 in steady stateoutgoing links onlyincoming links only disconnected giant SCC FIG. 11. Topological comparison between the 20 H network in Fig. 2(b) and various classes of random networks. (a) In-and out-degree distributions of the combustion network (blue and orange histograms, respectively), compared to the Poissondistribution (black curve). (b) Typical SCC structure of a random chemical reaction network with topological parametersmatching those of the combustion network. (c) Relative size of the giant SCC (blue) and relative number of nodes with P i > tion χ i to be non-uniform and often highly localized(as we observed for the hydrogen combustion examplein Sec. IV and attributed to the DAG structure of theSCC network in Sec. V), placing it off-center and closeto a hyper-edge of the high-dimensional simplex of prob-ability distributions. This tends to widen the quadrantsof growing deviation in the non-monotonicity conditionwe derived and thus increases the likelihood of observingan initial swing of the species concentrations away fromtheir steady-state values.These mechanisms for non-monotonic dynamics are in-timately related to the fact that many chemical reactions,particularly in combustion processes, are irreversible ornearly irreversible. In our formulation, reaction irre-versibility translates to local directionality of the corre-sponding links in the network. Since non-normality re-quires the network to be directed [5–7, 12, 17], some of thereactions need to be nearly irreversible for the two mecha-nisms mentioned above to induce non-monotonic dynam-ics. By removing the rare reverse transitions associatedwith nearly irreversible reactions, we reveal global di-rectionality in the hydrogen combustion networks in theform of a linear chain of strongly connected components.This directionality underlies the extreme localization ofsteady-state probability that promotes non-monotonicdynamics through the second mechanism above.These two mechanisms for non-monotonicity requiresnon-normality, which in turn requires local link direction-ality. However, our non-monotonicity condition clarifiesthat normal networks, even those with undirected links,can exhibit non-monotonic dynamics (but typically lesspronounced than for non-normal networks), dependingon the initial distributions of states and the choice ofthe observed quantity. The observed quantities we con-sidered here are similar to a weighted 1-norm and thusfundamentally different from the 2-norm of the devia- tion from the steady-state distribution, which can exhibitnon-monotonicity only for non-normal networks.The structure found in random networks is different(albeit still globally directional) and characterized by thedominance of a single giant connected component andeven more extreme localization of the steady-state prob-ability. This holds true also for a class of random reac-tion networks that respects the network structural con-straints imposed by the chemical master equation. Thus,the emergence of a linear chain in hydrogen combus-tion networks must be due to other factors not capturedby the random selection of reactions. One possibility isthat, when the product molecules are thermodynamicallymore stable than the reactants, this bias constrains theconnected component structure. While a model of ran-dom reaction networks accounting for this effect couldexplain the type of global directionality observed in hy-drogen combustion networks, different factors could leadto different global structures in other types of chemicalreaction networks, such as biological ones [47]. Manyintracellular biochemical reaction networks that exhibittransient dynamics, such as signaling networks, are ex-pected to have a high level of directionality. The non-monotonicity condition we derived here and our SCC-based analysis of the network’s global directionality areapplicable beyond those networks we explicitly consid-ered and lay a foundation for future research on non-normality and non-monotonic dynamics in complex re-action networks in general. Our approach, which focuseson transient dynamics far from equilibrium and can beextended to externally driven systems using time-varyingmaster equations, may extend to other network systemsand also contribute to the ongoing development of non-equilibrium statistical mechanics.4 ACKNOWLEDGEMENTS
This work was supported by MURI Grant No.W911NF-14-1-0359 and ARO Grant No. W911NF-19-1-0383.
APPENDIX
We prove that any right eigenvector of W associatedwith the (possible degenerate) zero eigenvalue can be cho-sen so that its components are all non-negative and sumto unity. We also prove that the components are zeroexcept for those corresponding to the nodes in the SCCsthat have no outgoing links. The latter implies that,for the steady-state distribution, all states outside theseSCCs have P i > c denote the number of SCCs with no outgoinglinks and c (cid:48) the number of the other SCCs. We first notethat the DAG structure can be used to re-index the nodesand make W a block lower-triangular matrix, with eachdiagonal block corresponding to an SCC. The matrix W then takes the following form: W = W (cid:48) ∗ . . . ∗ ∗ W (cid:48) c (cid:48) ∗ ∗ ∗ W ∗ ∗ ∗ . . . ∗ ∗ ∗ W c , (6)where an empty space indicates a zero block, and a starsymbol indicates that the block can have non-zero com-ponents. The c blocks corresponding to the SCCs with-out outgoing links necessarily appear as the last ones,denoted by W , . . . , W c , which form a block diagonalsubmatrix in W because there cannot be links connect-ing these SCCs by definition. Since the sum of com-ponents along each column of these blocks is zero, each block W i has a zero eigenvalue along with the associ-ated right eigenvectors. Since each SCC is by definitionstrongly connected, and hence each W i is irreducible, itfollows from the Peron-Frobenius Theorem that the zeroeigenvalue is not degenerate and that the correspondingeigenvector can be chosen to have only non-negative com-ponents. We choose such an eigenvector, and we furthernormalize it, so the sum of the component equals one.Extending these vectors to the size of the whole matrix W (setting all added components to zero), we obtainright eigenvectors associated with the zero eigenvalue of W .To see that the c -dimensional subspace spanned bythese eigenvectors covers the entire eigenspace associatedwith the zero eigenvalue, consider a column vector of W intersecting with the diagonal block W (cid:48) i , i = 1 , . . . , c (cid:48) .The sum of the off-diagonal components of W (cid:48) i thatappear in this vector is strictly less than the diagonalcomponent, since the sum of all the vector componentsis zero and all the off-diagonal components are non-negative. This implies that the Gershgorin circle cor-responding to this column of W (cid:48) i is to the left of theorigin of the complex plane and at a finite distance awayfrom the origin. Since this applies to all the columnsof W (cid:48) i for any i = 1 , . . . , c (cid:48) , zero cannot be an eigen-value of any of W (cid:48) , . . . , W (cid:48) c (cid:48) . Thus, all repetitions of thezero eigenvalue of W must come from the zero eigenval-ues of W , . . . , W c , and the c -dimensional subspace con-structed above is indeed the eigenspace of W correspond-ing to eigenvalue zero. We note that the eigenvectors of W spanning the eigenspace were chosen to have compo-nents that are all non-negative and are zero outside theSCCs with no outgoing links. Thus, any eigenvector of W associated with the zero eigenvalue shares the sameproperty and can thus be normalized so that its com-ponents sum to unity, giving the steady-state distribu-tion. Consequently, we have P i > [1] L. N. Trefethen, A. E. Trefethen, S. C. Reddy, and T.A. Driscoll, Hydrodynamic stability without eigenvalues,Science , 578 (1993).[2] L. N. Trefethen and M. Embree, Spectra and Pseudospec-tra: the Behavior of Nonnormal Matrices and Operators (Princeton University Press, Princeton, 2005).[3] M. G. Neubert and H. Caswell, Alternatives to resiliencefor measuring the responses of ecological systems to per-turbations, Ecology , 653–665 (1997).[4] T. Nishikawa and A. E. Motter, Synchronization is op-timal in nondiagonalizable networks, Phys. Rev. E, ,065106(R) (2006).[5] T. Nishikawa and A. E. Motter, Maximum performanceat minimum cost in network synchronization, Physica D , 77 (2006).[6] B. Ravoori, A. B. Cohen, J. Sun, A. E. Motter, T. E.Murphy, and R. Roy, Robustness of optimal synchro-nization in real networks, Phys. Rev. Lett. , 034102(2011).[7] G. Hennequin, T. P. Vogels, and W. Gerstner, Non-normal amplification in random balanced neuronal net-works, Phys. Rev. E , 011909 (2012).[8] S. Tang and S. Allesina, Reactivity and stability of largeecosystems, Front. Ecol. Evol. , 21 (2014).[9] T. Biancalani, F. Jafarpour, N. Goldenfeld, Giant Am-plification of Noise in Fluctuation-Induced Pattern For-mation, Phys. Rev. Lett. , 018101 (2017).[10] V. Klika, Significance of non-normality-induced patterns: Transient growth versus asymptotic stability, Chaos ,073120 (2017).[11] M. Asllani and T. Carletti, Topological resilience in non-normal networked systems, Phys. Rev. E , 042302(2018).[12] M. Asllani, R. Lambiotte, and T. Carletti, Structure anddynamical behavior of non-normal networks, Sci. Adv. ,eaau9403 (2018).[13] S. Nicoletti, D. Fanelli, N. Zagli, M. Asllani, G. Battis-telli, T. Carletti, L. Chisci, G. Innocenti, and R. Livi,Resilience for stochastic systems interacting via a quasi-degenerate network, Chaos , 083123 (2019).[14] R. Muolo, M. Asllani, D. Fanelli, P. K. Maini, and T. Car-letti, Patterns of non-normality in networked systems, J.Theor. Biol. , 81 (2019).[15] W. Tarnowski, I. Neri, and P. Vivo, Universal transientbehavior in large dynamical systems on networks, Phys.Rev. Res. , 023333 (2020).[16] G. Baggio, V. Rutten, G. Hennequin, and S. Zampieri,Efficient communication over complex dynamical net-works: The role of matrix non-normality, Sci. Adv. ,eaba2282 (2020).[17] S. Johnson, Digraphs are different: Why directional-ity matters in complex systems, J. Phys. Complexity ,015003 (2020).[18] D. A. McQuarrie, Stochastic approach to chemical kinet-ics, J. Appl. Probab. , 413 (1967).[19] D. T. Gillespie, A rigorous derivation of the chemicalmaster equation, Physica A , 404 (1992).[20] L. E. Reichl, A Modern Course in Statistical Physics , 4thedition (Wiley–VCH, Weinheim, Germany, 2016).[21] P. L. Krapivsky, S. Redner, and E. Ben-Naim,
A KineticView of Statistical Physics (Cambridge University Press,Cambridge, UK, 2010).[22] V. Seshadri, B. J. West, K. Lindenberg, Analytic theoryof extrema. III. Results for master equations with ap-plication to unimolecular decomposition, J. Chem. Phys. , 1145 (1980).[23] C. Jarzynski, Equilibrium free-energy differences fromnonequilibrium measurements: A master-equation ap-proach, Phys. Rev. E , 5018 (1997).[24] R. Albert and A.-L. Barab´asi, Statistical mechanics ofcomplex networks, Rev. Mod. Phys. , 47 (2002).[25] R. Pastor-Satorras and A. Vespignani, Epidemic Spread-ing In Scale-free Networks, Phys. Rev. Lett. , 3200(2001).[26] T. Hoffmann, M. A. Porter, R. Lambiotte, Generalizedmaster equations for non-Poisson dynamics on networks,Phys. Rev. E , 046102 (2012).[27] I. Mezi´c, Spectrum of the Koopman Operator, Spectralexpansions in functional spaces, and state-space geome-try, J. Nonlinear Sci. (2019); I. Mezi´c, Spectrum of theKoopman Operator, Spectral Expansions in FunctionalSpaces, and State Space Geometry, arXiv:1702.07597.[28] J. Li, Z. Zhao, A. Kazakov, and F. L. Dryer, An up-dated comprehensive kinetic model of hydrogen combus-tion, Int. J. Chem. Kinet. , 566 (2004).[29] M. ´O Conaire, H. J. Curran, J. M. Simmie, W. J. Pitz,and C. K. Westbrook, A comprehensive modeling studyof hydrogen oxidation, Int. J. Chem. Kinet. , 603(2004).[30] D. Baulch et al. , Evaluated kinetic data for combustionmodeling: supplement II, J. Phys. Chem. Ref. Data , 757 (2005).[31] A. A. Konnov, Yet another kinetic mechanism for hydro-gen combustion, Combust. Flame , 14 (2019).[32] See our Github repository at https://github.com/znicolaou/ratematrix for data sets, as well as code togenerate the state space, the ensemble evolution, and theconnected components for any reaction network data filein the Chemkin or Cantera format.[33] A. R´enyi. On measures of entropy and information, In Proceedings of the Fourth Berkeley Symposium on Math-ematical Statistics and Probability , vol. 1 (The Regentsof the University of California, Oakland, 1961).[34] R. C. Tolman,
The Principles of Statistical Mechanics (Dover Publications, New York, 2010).[35] D. G. Goodwin, H. K. Moffat, and R. L. Speth, Cantera:An object-oriented software toolkit for chemical kinetics,thermodynamics, and transport processes, , 2017.[36] R. J. Kee, J. A. Miller, and T. H. Jefferson, CHEMKIN:A general-purpose, problem-independent, transportable,FORTRAN chemical kinetics code package, No. SAND–80-8003, Sandia Labs., 1980.[37] The number N of states accessible from the 2 m H initialstate can be shown to be upper-bounded by a power lawwith an exponent of 8. To see this, we first note thatany accessible state has the same number of H and Oatoms as the initial state, which has (4 m + 1) H atomsand (2 m + 1) O atoms. Here, we do not need to considerthe Ar atoms since the Ar molecules do not change inthe process. Then, the maximum number of molecules ofspecies H , O , H, O, OH, H O, HO , and H O that anystate can have are 2 m , m , 4 m +1, 2 m +1, 2 m +1, 2 m , m ,and m , respectively, implying that N ≤ m · m · (4 m +1) · (2 m +1) · (2 m +1) · m · m · m = 4 m (2 m +1) (4 m +1) ∼ m .This bound, however, is typically far larger than N . Forexample, m = 3 yields N = 347 (cid:28) m (2 m + 1) (4 m +1) = 619 , m = 3, this bound is 842, as mentioned earlierin the text.[38] In contrast to existing methods (e.g., the finite stateprojection [48] and stochastic simulations [49, 50]), ourthresholding approach offers a scalable alternative thatcan both account for link weights W ij and calculate ac-curate direct solutions of the master equation.[39] We find that the maximum error in P i ( t ) caused by the (cid:15) -thresholding over all i and over all t between t = 10 − secand t = 10 sec (well after the combustion process isconsidered complete) is small and generally decreaseswith the network size. Already for the network of size N = 4 ,
184 built from the 16 H initial state (the largestnetwork for which the error could be computed), themaximum error is less than 2 × − .[40] R. A. Horn and C. R. Johnson, Matrix Analysis , 2ndedition (Cambridge University Press, New York, 2013).[41] M. K. Transtrum, B. B. Machta, K. S. Brown, B. C.Daniels, C. R. Myers, and J. P. Sethna, Perspective:Sloppiness and emergent theories in physics, biology, andbeyond, J. Chem. Phys. , 010901 (2015).[42] A. White, M. Tolman, H. D. Thames, H. R. Withers,K. A. Mason, and M. K. Transtrum, The limitationsof model-based experimental design and parameter es-timation in sloppy systems, PLOS Comput. Biol. ,e1005227 (2016). [43] Note that the Shannon and R´enyi entropies in this set-ting are distinct from the Boltzmann entropy, which in-cludes uncertainty about the positions and velocities ofthe molecules in addition to the chemical composition.While the Shannon and R´enyi entropies may exhibit non-monotonic evolution, the Boltzmann entropy strictly in-creases according to the second law of thermodynamics.[44] S. N. Dorogovtsev, J. F. F. Mendes, and A. N. Samukhin,Giant strongly connected component of directed net-works, Phys. Rev. E , 025101 (2001).[45] M. E. J. Newman, S. H. Strogatz, and D. J. Watts, Ran-dom graphs with arbitrary degree distributions and theirapplications, Phys. Rev. E , 026118 (2001).[46] F. Chung, L. Lu, and V. Vu, Spectra of random graphswith given expected degrees, Proc. Natl. Acad. Sci. U.S.A. , 6313 (2003).[47] K. Kaneko, Life: An Introduction to Complex SystemsBiology (Springer–Verlag, Berlin, Germany, 2006).[48] B. Munsky and M. Khammash, The finite state projec-tion algorithm for the solution of the chemical masterequation, J. Chem. Phys. , 044104 (2006).[49] L. B. Newcomb, M. Alaghemandi, and J. R. Green.Nonequilibrium phase coexistence and criticality near thesecond explosion limit of hydrogen combustion, J. Chem.Phys. , 034108 (2017).[50] L. B. Newcomb, M. E. Marucci, and J. R. Green. Explo-sion limits of hydrogen–oxygen mixtures from nonequilib-rium critical points, Phys. Chem. Chem. Phys.20