A closed form for Jacobian reconstruction from timeseries and its application as an early warning signal in network dynamics
Edmund Barter, Andreas Brechtel, Barbara Drossel, Thilo Gross
AA closed form for Jacobian reconstruction fromtimeseries and its application as an early warningsignal in network dynamics
Edmund Barter ∗ , Andreas Brechtel ∗ , Barbara Drossel and Thilo Gross University of Bristol, Department of Engineering Mathematics, MerchantVenturers Building, Woodland Road, Bristol BS8 1UB, UK TU Darmstadt, Fachbereich Physik, Hochschulstrasse 12, 64289 Darmstadt,Germany UC Davis, Department of Computer Science, 1 Shields Av. Davis Ca 95616, USA @ please address correspondence to [email protected] ∗ these authors contributed equally Abstract
The Jacobian matrix of a dynamical system describes its responseto perturbations. Conversely one can estimate the Jacobian matrix bycarefully monitoring how the system responds to environmental noise.Here we present a closed form analytical solution for the calculation ofa system’s Jacobian from a timeseries. Being able to access a system’sJacobian enables us to perform a broad range of mathematical analysesby which deeper insights into the system can be gained. Here we considerin particular the computation of the leading Jacobian eigenvalue as anearly warning signal for critical transition. To illustrate this approach weapply it to ecological meta-foodweb models, which are strongly nonlineardynamical multi-layer networks. Our analysis shows that accurate resultscan be obtained, although the data demand of the method is still high.
As humans we are dependent on the functioning of complex systems on manyscales, ranging from our own body with its interlinked metabolic, signallingand microbial networks, via technical and organizational systems such as powergrids and political systems, to the planetary-scale supply chains and the climatesystems. All of these are complex nonlinear many variable systems, and assuch at a risk of undergoing sudden, qualitative and potentially irreversibletransitions [1]. Over the past decades there has been a steadily growing interestin methods that provide early warning of such transitions [2, 3]. In particular1 a r X i v : . [ n li n . C D ] O c t here is a growing awareness that such methods are needed for systems that arespatial or network based and hence inhherently high-dimensional [4, 5, 6].The traditional approach to anticipating qualitative transitions in complexsystems is mechanistic modelling. Many well-understood systems can be mod-elled so precisely that the model can accurately predict the threshold parametervalues where transitions occur. For technical applications such as aircraft flightor stability of structures the model-based stability analysis is well establishedand in many cases part of a legally mandated licensing process [7]. However,model-based approaches tend to produce poor results in systems that are lesswell known as seemingly minor details of the model can sometimes drasticallyaffect transition points.It has long been known that critical transitions are generally preceded bycritical slowing down [8, 9]. This phenomenon leads to a distinctive increase inthe auto-correlation and cross-correlations of timeseries before a transition. Theadvantage of correlation-based warning signs is that detailed understanding ofthe system is not needed. Its disadvantage is that high frequency time seriesdata are necessary for robust warning, which imposes a strong, and for manyapplications prohibitive, constraint.The model-based and correlation-based approaches to predicting criticaltransitions can be seen as two extreme strategies. The former is entirely basedon structural knowledge of the system and uses real world data at most to fitmodel parameters. The latter eliminates the need for structural information, atthe cost of a high time series data demand [10, 11].In many potential applications there are aspects of the system that are wellunderstood because different variables are related via physical laws or subjectto logical constraints [11]. If such bits of ‘structural’ knowledge are available ina system it is desirable to exploit them for the construction of early warning sig-nals. However, the same systems may also contain aspects that are considerablyless well understood and hence make purely model based predictions unreliable.This defines the need for a middle way, where available structural informationon a system is used to reduce the data demand, while time series data is usedto close gaps in the understanding in areas where structural information is notavailable [2].A common approach to finding a middle way in the construction of earlywarning signs is to use data assimilation approaches to continuously improvea dynamical model [11]. A classical application of these approaches is the eu-trophication of shallow lakes [12, 13], for which good early warning signs canbe constructed using techniques such as Bayesian learning and Kalman filters[14, 15, 16]. However, a broad investigation of such approaches shows that theymay result in warning signs that are “faint and late”[11], while Boettiger andHastings point out advantages of a model-based approach [13].A different approach to critical transitions builds on linear stability analysis[17, 18]. In [17], the authors formulated a model that was sufficiently generalto encompass the whole class of conceivable models into which a given systemcould potentially fall. Using the so-called generalized modelling approach, theJacobian matrices were computed that govern the system’s response to pertur-2ations and from which the bifurcation points can be computed. Because theunderlying models contain unknown functional relationships the Jacobian ma-trices that are thus obtained still contain unknown parameters. The authorsused time series data to eliminate these remaining uncertainties. They showedthat in this way accurate early warning signals can be constructed that onlyrequire very limited timeseries data.The present paper takes the idea of reconstructing a systems Jacobian fromdata one step further by avoiding assumptions on the functional form of theinteraction terms. In a system that is subject to some noise, the cross correla-tions in timeseries encode very similar information to the Jacobian matrix of itsdeterministic backbone. Considered in isolation the cross correlations are notsufficient to reconstruct the full Jacobian, however a complete reconstruction ispossible if we have some additional structural information. In systems that canbe described as complex networks the network structure imposes constraints onwhich variables can interact directly, which in turn implies that some entries ofthe Jacobian must be zero. In a sufficiently sparse network, and particularlyin multi-layer networks, knowing these zeroes provides sufficient information toreconstruct the remainder of the Jacobian from time series correlations.For illustration we apply Jacobian reconstruction approach to an ecologi-cal meta-foodweb model, formulated as a dynamical multi-layer network. Bycomparing with the known ground truth of the model, i.e., its exact Jacobian,we show that the Jacobian can be reconstructed faithfully and demonstrate itsvalue as an early warning signal. We find that despite leveraging the structuralinformation, the amount of timeseries data required for accurate results is atpresent still prohibitively high for the ecological application. However we dis-cuss several avenues of future research that may reduce the data requirementsto a point where the method becomes widely applicable. In every dynamical system that is in the vicinity of some form of long-termbehaviour, the response of the system to small perturbations in the variablescan be captured by some matrix. In the simplest case the system is a system ofordinary differential equations (ODEs)˙ x = f ( x , p ) (1)in which a vector of variables x evolves in time in a way that is dependent ona set of parameters p . The simplest form of long-term behavior is rest in astationary state x ∗ ( p ), such that f ( x ∗ , p ) = 0 (2)The response to sufficiently small perturbations in the the steady state is thendescribed by the Jacobian matrix J , whose elements are computed as J ij = ∂∂x j ˙ x i (cid:12)(cid:12)(cid:12)(cid:12) x = x ∗ (3)3he steady state under consideration is stable when all eigenvalues of the Ja-cobian have negative real parts. A smooth change in the parameters p willgenerally cause the the eigenvalues to change smoothly. When such a changecauses one or more eigenvalues to acquire positive real parts, then a bifurcationoccurs in which the dynamics change qualitatively.Because the Jacobian is a real matrix its eigenvalues are real or form com-plex conjugate pairs. Hence there are two fundamental types of bifurcations. Inbifurcations of fold-type a single real eigenvalue acquires a positive real part asit passes through zero. Several different forms of this bifurcation are commonlyencountered in dynamical systems (fold bifurcation, pitchfork bifurcation, tran-scritical bifurcation), but generally speaking these are associated with a changein the number of steady states in the system or the exchange of stability prop-erties.In a Hopf bifurcation a complex conjugate eigenvalue pair acquires a positivereal part by crossing the imaginary axis of the complex plane. This bifurcationis generally associated with the onset of at least transient oscillations as thesystem departs from the steady state.Both types of bifurcations can occur in several different forms, some of whichcause only relatively mild non-catastrophic transitions (say, replacing stationarybehaviour with low amplitude oscillations), while others lead to a catastrophic(and potentially irreversible) departure from the steady state.One can distinguish between different types of bifurcations by computingnormal form parameters that are functions of higher derivatives of f . However,this is beyond the scope of the present paper.In spatially extended systems, which are defined on a continuous space or ona spatial network of discrete nodes, the fundamental bifurcations can come intwo flavors: A bifurcation may affect all points in space simultaneously in thesame way, or it may affect points in space differently leading to the formation ofspatial patterns. For clarity the bifurcations of the latter type are called Turingbifurcations (fold-like case) and wave instabilities (Hopf-like case).While not every bifurcation in a dynamical system is a critical transition,any bifurcation occurring in an important real world system is certainly a causefor concern. In this spirit our aim in the following is to reconstruct the Jacobianof a dynamical system from data in order to determine its leading eigenvalue. Ifthe real part of this eigenvalue approaches zero, we interpret this as a warningsignal for an impending bifurcation. Our aim in this section is to formulate a method that can reconstruct the Jaco-bian matrix of a dynamical system from time series. We do this by expandingon the work of Honerkamp [19], van Kampen [20], and Steuer et al. [21].We start from a stochastic timeseries that fluctuates around a steady state x ∗ of the underlying deterministic backbone of the system. As shown in [19, 20, 21]and reproduced in appendix B the Jacobian J close to x ∗ is related to the4ovariance matrix Γ of timeseries, and to the fluctuation matrix D of the noise,via the equation JΓ + ΓJ T = − D . (4)In the following, we show that we can use this relationship to compute J . Con-sider that a Jacobian matrix of linear dimension N contains N independententries. By contrast, Eq. (4) equates two symmetric matrices, and hence imposesonly N ( N + 1) / J . Therefore in any applicationwith multiple variables ( N >
1) the system is underdetermined such that wecannot recover the complete Jacobian purely from the time series.We can still recover the full Jacobian if we have additional information thatwe can leverage. Fortunately, in many applications some structural informationis easily accessible. In particular in large spatially complex systems or reason-ably sparse networks we know that certain variables cannot interact and hencethe corresponding elements of the Jacobian must be zero. This yields a set ofstructural constraints of the form J ij = 0 for certain pairs ( i, j ). If we can iden-tify at least N ( N − / N variables and a set of G additionalstructural constraints, with G ≥ N ( N − / (cid:18) J J J J (cid:19) (cid:18) Γ Γ Γ Γ (cid:19) + (cid:18) Γ Γ Γ Γ (cid:19) (cid:18) J J J J (cid:19) = − (cid:18) D D D D (cid:19) which implies the independent conditions2 J Γ + 2 J Γ = − D (5) J Γ + J Γ + J Γ + J Γ = − D (6)2 J Γ + 2 J Γ = − D (7)with the second condition applying to both off-diagonal terms. The left-handside of these equations is a linear system. Hence we can write the conditions inthe form B j = − d (8)where B is a matrix, j is a column vector that contains the entries of theJacobian, i.e. j = ( J , J , J , J ) T and d is the corresponding vector for D .For the two dimensional example this reads Γ Γ Γ Γ Γ Γ Γ J J J J = − D D D D (9)5he form of this equation suggests that we can solve it for j by multiplying B − from the left. However, we have to take care because B is not invertible becausethe two center rows are identical, which is a consequence of the missing infor-mation. We can fix this problem by using the additional constraints. Imposingstructural constraints on the two-variable system makes this example almostpointless, but, for the purpose of illustration, let us assume that we know thatvariable 1 cannot depend on variable 2 and hence J = 0. We can representthis constraint in the same matrix equation as the system by adding it as anadditional line, Γ Γ Γ Γ Γ Γ Γ J J J J = − D D D D (10)If this is the only constraint then we can now drop the third row of the matrixand the third entry of the vector on the right-hand side and solve for j by matrixinversion. In practice there are typically additional constraints and hence thesystem will be overdetermined. In this case we use least squares optimizationto find an approximate solution. As we will see below the form of Eq. (10) isvery convenient for finding the least squares solution. In particular it allows usto obtain a analytical closed form solution for j .Let us now generalize from the two-dimensional example to systems withan arbitrary number of variables. For this purpose we define the vectorizationoperator [22] vec( X ) = ( X , X , . . . X N , X , . . . ) T (11)We can now write Eq. (8) as B vec( J ) = − D ) (12)i.e. the vectorization of a matrix is the columns of the matrix stacked on top ofeach other. To find the general form of B we start from Eq. (4) and vectorizeboth sides, which yields vec( JΓ + ΓJ T ) = vec( − D ) . (13)Because vectorization is a linear operator we can pull the -2 out of the vector-ization on the right hand side and apply the vectorization separately to the twoterms on the left hand side, hencevec( JΓ ) + vec( ΓJ T ) = − D ) . (14)It is known [22] that for matrices X , Y , Z the following identity holds (cf. ap-pendix C): vec( XYZ ) = ( Z T ⊗ X ) vec( Y ) (15)6here ⊗ is the Kronecker product of matrices defined by X ⊗ Y = X Y X Y . . .X Y X Y . . . ... ... . . . (16)We now substitute Y = J , Z = Γ and X = I , where I is the identity matrix ofappropriate size. This yieldsvec( JΓ ) = vec( IJΓ ) = ( Γ ⊗ I ) vec( J ) (17)which brings the first term from Eq. (14) into the desired form. If we try thesame for the second term we findvec( ΓJ T ) = vec( ΓJ T I ) = ( I ⊗ Γ ) vec( J T ) (18)which is not quite the desired form because vectorization of the transpose of J appears rather than the vectorization of J itself. The vectorization of a ma-trix and the vectorization of its transpose are not identical but closely related.Consider that for our two-dimensional the two vectorizations are related byvec( J T ) = J J J J = C J J J J = C vec( J ) (19)where C = (20)is a permutation matrix. In the general case we can still writevec( J T ) = C vec( J ) (21)where C is now a permutation matrix of size N × N , one can construct thismatrix from blocks of size N × N as C = C C . . . C C . . . ... ... . . . (22)Here, C nm is a N × N -matrix defined by( C nm ) ij = δ im δ nj (23)where δ is the Kronecker delta. 7sing the commutation matrix we can write Eq. (18) asvec( ΓJ T ) = ( I ⊗ Γ ) vec( J T ) = ( I ⊗ Γ ) C vec( J ) (24)Substituting this relationship and Eq. (17) into Eq. (14) we get the desired form(( Γ ⊗ I ) + ( I ⊗ Γ ) C ) (cid:124) (cid:123)(cid:122) (cid:125) B vec( J ) (cid:124) (cid:123)(cid:122) (cid:125) j = − D ) (cid:124) (cid:123)(cid:122) (cid:125) d . (25)We now know how to construct the matrix B such that we can write the systemin the form B j = − d . (26)However, this system is still underdetermined. If we know that some elements of j must be zero we can represent this knowledge in a matrix U . For this purposewe can gather the respective elements of j in an ordered list U and then define U as an |U| × N matrix with U nm = (cid:26) j m is the n ’th entry of U j must be zero.In analogy to the small example we can now impose the additional conditionson the system by stacking them below the matrix B such that Eq. (26) becomes (cid:18) BU (cid:19) j = − (cid:18) d (cid:19) (28)where is a column vector containing |U| zeroes. The () that appear in thisequation should be read as a block-wise notation for matrices/vectors, whererows are stacked on top of each other in analogy to Eq. (10).To simplify the notation we introduceˆ B = (cid:18) BU (cid:19) ˆ d = (cid:18) d (cid:19) (29)which allows us to write the whole set of conditions once again in the formˆ B j = ˆ d (30)In practise this will now be an overdetermined system, such that no exact so-lution exists. However, finding an approximation that minimizes the squares ofthe deviations in each row is a well-known problem. The known solution [23, 24](cf. appendix D) for this problem is j = ( ˆ B T ˆ B ) − ˆ B T ˆ d (31)where the expression ( ˆ B T ˆ B ) − ˆ B T , the pseudoinverse of ˆ B ,appears.8he equation provides a closed form solution by which the Jacobian elementscan be computed from the covariance matrix of the timeseries, a known orestimated fluctuation matrix, and a set of additional structural constraints onthe Jacobian.Typically the least squares fit will not set the Jacobian elements governedby the structural constraints exactly to zero. One may therefore enforce theseknown zeros by setting the respective elements of the Jacobian explicitly to zeroafter the computation of Eq. (31) has finished. In numerical experiments, de-scribed below, we found that this improved the accuracy of Jacobian eigenvaluesestimated by this method.To summarize, the Jacobian J of an N -dimensional system close to a steadystate can be reconstructed as follows:1. Compute the co-variance matrix from the time series data, Γ ij = (cid:104) X i X j (cid:105) .2. Construct the diagonal fluctuations matrix D , and compute d = vec( D ).In some systems these fluctuations can be measured directly, otherwise areasonable approximation may be derived based on assumptions on theunderlying noise process [20].3. Construct the permutation matrix C according to Eqs. (20,21).4. Compute the matrix B = ( Γ ⊗ I + ( I ⊗ Γ ) C ) , where I is the N × N identity matrix.5. Define j = vec( J ) and identify at least N ( N − / j that mustbe zero due to structural constraints. Use these to construct the matrix U (see Eq. (27)). The column-dimension of U is the row dimension of j , androw dimension of U is identical to the number of structural constraints.The matrix has exactly one non-zero entry in each row such that U nm = 1if the n th structural condition reads j m = 0.6. Construct ˆ B = (cid:18) BU (cid:19) ˆ d = (cid:18) d (cid:19)
7. Compute j = ( ˆ B T ˆ B ) − ˆ B T ˆ d and recover the Jacobian J from its vectorization j .8. Set the elements of J governed by structural constraints explicitly to zero.9 Application to a meta-foodweb model
In the following we explore if the Jacobian matrix can be reconstructed suffi-ciently accurately to warn of impending critical transitions. For this purpose, weconstructed a test system of realistic complexity for which the Jacobian matrixis nevertheless known analytically, and we created noise timeseries.We used a meta-foodweb model already studied in [25, 26] (see appendixE for details). The model consists of a spatial network of P habitat patcheslinked by avenues of species dispersal (Fig. 1). Each patch harbours a complexfoodweb, consisting of S nodes that represent populations of different species,which are linked by predator-prey interactions. The dynamics of the system aregiven by a set of differential equations that govern the changes in variables dueto diffusive dispersal between patches and biological processes occurring withina patch (primary production, predator-prey interaction, natural mortality).The model contains several complicating factors encountered in real systems.The functions used to model the biological processes are strongly nonlinear.They capture various saturation effects and realistic responses to the availabilityof different food sources (prey switching). The dynamics of different speciesoccurs on different time scales according to biological scaling relationships whichrelate a species position in the foodweb to its expected biomass turnover rate[27]. Similar scaling relationships also govern the rate at which different speciesdisperse across the spatial network[25].We consider two different versions of the patch topology. The smaller of thetwo consists 6 patches which are connected in such a way that they form thesmallest completely asymmetric network (Fig. 1D). The larger one contains 15patches and was generated as a random geometric graph (Fig. 1E). It has a com-paratively large diameter and high clustering and contains several symmetriesthat are characteristic for this type of spatial networks [28]. For the food webswe used a predator-prey system consisting of two species (Fig. 1C), as well asa four-species system in the shape of the so-called intra-guild predation motif(Fig. 1B), a common and well-studied foodweb motif.Work by Nakao and Mikhailov [29] and the extension of their approachto meta-foodwebs [25] showed that network models behave analogously to dy-namical systems in continuous space. Hence the theory of pattern formationin partial-differential equation systems can transferred almost exactly to thesedynamical networks. This means that our test systems can exhibit instabili-ties that are best described as pattern-forming bifurcations, specifically Turingand Wave instabilities. Pattern forming instabilities in predator-prey systemsin continuous space were studied in detail in [30] which made it easy to locatethese bifurcations also in our multi-layer network predator-prey system (Fig. 2).Using a generalized modelling approach [31, 32, 33] we analytically computedthe Jacobian that describes the Jacobian matrices of the two example food webson arbitrary spatial topologies. We then picked specific realizations of models inwhich relevant bifurcations occurred. To produce noisy timeseries we simulatedthese models with added noise using the Euler-Maruyama method (see appendixA). 10igure 1: Schematic representation of the model system. We consider a mul-tilayer network where copies of an ecological food web exist in different geo-graphical patches (A). We consider an intra-guild predation food web with anadditional predator (B) and a predator-prey system (C). For the spacial net-work, we use the smallest completely asymmetric graph (D) and a larger randomtopology, generated as a random geometric graph(E). As a first test we generated simulated noisy timeseries containing 2 · datapoints from the two-species system on the six-patch topology. Parameter valuesfor these simulations were chosen to lie on four transects through the parametersspace that crossed bifurcation points. For each of these time series we thenreconstructed the Jacobian matrix and computed the leading eigenvalue. Beforethe bifurcation point the estimated eigenvalues are in very good agreement withthe known ground truth provided by the analytic eigenvalues (Fig. 3).As we follow the transects the real part of the leading eigenvalue crosseszero in a bifurcation. When the bifurcation occurs the reconstructed eigenvaluedeparts from the analytical value. This behavior is expected as the analyticalsolution continues to show the eigenvalues around the, now unstable, steady11 .4 0.5 0.6 0.7 0.8 0.9 1.0 φ γ stable AT B (a) (b) (c) (d) TuringHopfTrascritical
Figure 2: Bifurcation diagram of the predator-prey system on the six patchestopology. Stability of the state under consideration changes in response tochanges in the sensitivity of biomass production to producer biomass φ and thesensitivity of predation to prey biomass γ . The state under consideration is sta-ble in the top left area. Stability is lost when either of three bifurcations occur(Turing, Transcritical, Hopf). After the loss of stability the system approaches astate of homogeneous oscillations (A), a different homogeneous stationary state(B), or a state of stationary patterns (T). The bifurcation diagram was com-puted using the master stability function approach from [25] (see appendix). Itcorresponds directly to Fig. 1 from [30] which studies a predator-prey system incontinuous space. Lines (a-d) indicate the transects used for the correspondingsimulations in Fig. 3.state, whereas the reconstruction algorithm computes the leading eigenvalueassociated with the new dynamics, which has now departed from the previoussteady state.We note that in transect d) the reconstructed eigenvalue is almost exactlyzero in a wide region after the bifurcation. This happens because the systemapproaches a stable limit cycle for which the leading Lyapunov exponent iszero. The recovery of this zero by the algorithm provides some (unexpected)evidence suggesting that the method reveals some salient information even innon-stationary states.Careful examination of the transects (Fig. 3) suggests that the accuracyof reconstruction improves as we approach the bifurcation point. To explorethis further we generated 15 transects in the vicinity of bifurcation points and12igure 3: Comparison of the analytical ground truth for the leading eigenvaluewith an eigenvalue estimate from Jacobian reconstruction. The panels corre-spond to the 4 transects shown in Fig. 2. Estimates are in good agreementwith the analytical value until stability is lost and the timeseries depat fromthe steady state. [Parameters φ, γ : a) 0.75, 0.38 to 0.8, 0.33; b) 0.85, 0.46to 0.9, 0.41; c) 0.9, 0.51 to 0.95, 0.46; d) 0.95, 0.56 to 1.0, 0.51. Bifurcationsencountered are Turing (a,b) and Hopf (c,d) ]generated three sets of timeseries along every one of the transects (see appendixfor details). The results confirm that the accuracy of the estimate improves asthe bifurcation point is approached (Fig. 4).We now consider the case where parameters are slowly changing over a longsimulation run. Our aim is to estimate the leading eigenvalue of the Jacobianover time as this slow change in the system is taking place. For this purposewe apply the proposed method to reconstruct the Jacobian in sliding time win-13igure 4: The plot shows the error of the eigenvalue approximation in depen-dence to the analytic eigenvalue for the S = 2 species on N = 6 patch predatorprey system. The error is the difference of the estimated eigenvalue and theanalytic one. To gather the data 15 parameter trajectories covering transcriti-cal, Hopf and Turing bifurcations where sampled 3 times each. The errors aresmall for negative leading eigenvalues. The errors for greater analytic eigen-values where cut, because the estimation is not expected to work for positiveeigenvalues.dow of length τ . The results (Fig. 5) show that the estimates based on thesliding window are more noisy, undergoing visible fluctuations around the truevalue. Nevertheless the trend of the eigenvalue approaching zero is still clearlycaptured.To test the applicability to larger networks, we apply Jacobian reconstructionto a model system with 4 species on 15 patches. Even in this larger systemthe accuracy is still quite good but the estimated eigenvalue is systematicallyslightly less than the true value. For the case of fixed parameter values thedifference once again disappears as we approach the bifurcation. But for thecase of sliding parameter values a small difference remains. We suspect thatthis may be the combined effect of the non-autonomous nature of this systemsand the noise leading to bifurcation delay. This delay effect can be expectedto be more pronounced in the larger food web due to the presence of higher-level predator whose dynamics happen on correspondingly longer timescales. Ifthis is the case then the reconstructed Jacobian eigenvalue may actually offer abetter estimate of the relevant transition point than the analytical solution forthe system without noise. 14igure 5: Eigenvalue estimation in a system with continuously changing param-eters. We compare the estimated eigenvalue from Jacobian reconstruction in asliding window to the analytic ground truth for a predator-prey system on 6geographical patches close to a Turing bifurcation. The panels show the sameresult, however the estimate is placed at the end of the observation window(panel a) or in the center (panel b). For large φ the steady state described bythe analytical eigenvalues is unstable (greyed region) and hence the reconstruc-tion yields eigenvalues of a different state. The width of the sliding window isindicated on the left (shaded blue). [Parameters: φ, γ is changed from 0.70,0.35to 0.75,0.30]So far we studied systems that were designed using the generalized modellingapproach. We complement this by the analysis of a well established ecologi-cal model, the Rosenzweig-MacArthur predator-prey model [34] with quadraticmortality and diffusion on the 6-patch network (see appendix for details). Theresults of Jacobian reconstruction (Fig. 7) show that the leading eigenvaluecan be recovered with reasonable accuracy, again the accuracy of the estimateimproves significantly as the system approaches the bifurcation point. In this paper we expanded on previous work by Honerkamp, van Kampen, Steuerand others to derive a closed form expression for the reconstruction of Jacobianmatrices from time series data.For illustration we applied the mathematical formula to the an ecologicalmeta-foodweb model. This example illustrated that a relatively robust recon-struction of the leading eigenvalue of the Jacobian is possible even in a stronglynonlinear multi-layer network with dynamics on multiple timescales. However,the example also revealed that the required amount of data is still prohibitive15igure 6: Eigenvalue estimations for the system with 4 species on 15 patches.We consider loss of stability as a result of the change of two different param-eters, governing the nonlinearity of primary production, φ (panel a,c) and thenonlinearity of the functional response of attack rate to prey density, γ (b,d).The leading eigenvalue of the Jacobian was reconstructed from timeseries forfixed parameters (a,b) and parameter transsects (c,d). The steady state un-der consideration is unstable subsequent to a Hopf bifurcation (region shadedgrey). For the transects estimates are shown in the center of the sampling win-dow (window width is indicated by blue area). See appendix for parameters anddetails.for the ecological application.There is reason to believe that future research and particular a deeper math-ematical understanding of the Jacobian reconstruction can significantly reducethe required amount of data. Particularly interesting in this respect is the ob-16igure 7: Eigenvalue estimations in the Rosenzweig-MacArthur with quadraticmortality and diffusion on the six-patch network. The steady state under con-sideration loses stability due to a Fold (panels a,c) and Hopf (b,d) bifurcations.The leading eigenvalue was reconstructed from simulation runs with fixed pa-rameters (a,b) and slowly changing parameters (c,d). The accuracy of the esti-mates improves as the system approaches the bifurcation. (Estimated values areshown in the middle of the sampling window, indicated in blue. See appendixfor details)served increase in accuracy close to the bifurcation point. A promising goal forfuture exploration would be to understand how far this region of heightenedaccuracy extends. If a real world application under consideration is alreadyclose to bifurcation one might find that much less data is required to reach thedesired accuracy.Alternatively, it might be possible to reduce the data demand by optimizing17he sampling scheme. One can envision an iterative scheme, similar to [35],where a small number of samples is used to find an initial estimate of the Jaco-bian. Using the initial estimate one could then identify the relevant timescalesand important entries in the covariance matrix and optimize the sampling effortaccordingly.A third alternative may be to use additional information that may be avail-able. The advantage of our jacobian-based approach is that it can take advan-tage of additional knowledge on the systems that may be available in ecologicalapplications, such as closure exponents or the predator-dependence of the pre-dation rate. The previous example [17] suggests that the use of this informationmay reduce the data demand considerably.Meanwhile the method proposed here may be useful in fields where data ismore readily available, such as studies of metabolism, power grids, or economicdata. In the study of metabolism Jacobian reconstruction is already frequentlyused [36], for this application the present work yields an analytic closed formsolution to a problem that is so far solved by machine learning methods. Forpower grids, reconstructing Jacobians may be particularly interesting becauseit could yield deeper insights into the functioning of the system in additionto providing an early warning signal. For economics it may be particularlyinteresting that the Jacobian can be seen as a representation of causality in thesystem. Closed form jacobian reconstruction thus offers a way to infer causalityfrom correlation and is guaranteed to be exact in the large data limit. Animportant caveat is that the method requires some additional information toavoid underdeterminedness. However, as illustrated here, it is already sufficientif we can identify N ( N − / N variables of the system. In networks language, this is equivalent to sayingthat the mean degree z of the network of interactions must obey z < N − , (32)a condition that should be easy to satisfy in many applications and becomeseasier to meet in larger systems.In summary, we find that Jacobian reconstruction is a promising approachto the analysis of complex systems near critical transitions, although the datarequirements presently still limit its applicability. We expect that the closedform solution derived in the present paper inspires future mathematical workto alleviate these requirements. Acknowledgements
The authors thank Caio Lucidius N. Azevedo and Nicolas Verschueren van Reesfor valuable pointers to the algebra literature. This work was supported byDeutsche Forschungsgemeinschaft Research Unit FOR 1742, under grant con-tract number Dr300/13-2 and by the EPSRC (EP/N034384/1).18 eferences [1] M. Scheffer,
Critical transitions in nature and society . Princeton: PrincetonUniversity Press, 2009.[2] M. Scheffer, S. R. Carpenter, T. M. Lenton, J. Bascompte, W. Brock,V. Dakos, J. van de Koppel, I. A. van de Leemput, S. A. Levin, E. H. vanNes, M. Pascual, and J. Vandermeer, “Anticipating critical transitions,”
Science , vol. 338, no. 6105, pp. 344–348, 2012.[3] C. Kuehn, G. Zschaler, and T. Gross, “Early warning signs for saddle-escape transitions in complex networks,”
Scientific Reports , vol. 5, no. 1,p. 13190, 2015.[4] S. K´efi, V. Guttal, W. A. Brock, S. R. Carpenter, A. M. Ellison, V. N.Livina, D. A. Seekell, M. Scheffer, E. H. van Nes, and V. Dakos, “Earlywarning signals of ecological transitions: Methods for spatial patterns,”
PLoS ONE , vol. 9, no. 3, p. e92097, 2014.[5] J. Ludescher, A. Gozolchiani, M. I. Bogachev, A. Bunde, S. Havlin, andH. J. Schellnhuber, “Very early warning of next el ni˜no,”
Proceedings of theNational Academy of Sciences , p. 201323058, 2014.[6] G. Tirabassi, J. Viebahn, V. Dakos, H. Dijkstra, C. Masoller, M. Rietk-erk, and S. Dekker, “Interaction network based early-warning indicators ofvegetation transitions,”
Ecological Complexity , vol. 19, pp. 148–157, 2014.[7] D. Raymer,
Aircraft Design: A Conceptual Approach . Reston: AmericanInstitute of Aeronautics and Astronautics, Inc., 2018.[8] H. E. Stanley,
An introduction to phase transitions and critical phenomena .Oxford: Oxford University Press, 1987.[9] C. Kuehn, “A mathematical framework for critical transitions: Bifurca-tions, fast-slow systems and stochastic dynamics,”
Physica D , vol. 240,no. 12, pp. 1020–1035, 2011.[10] V. Dakos, S. R. Carpenter, W. A. Brock, A. M. Ellison, V. Guttal, A. R.Ives, S. K´efi, V. Livina, D. A. Seekell, E. H. van Nes, and M. Scheffer,“Methods for detecting early warnings of critical transitions in time se-ries illustrated using simulated ecological data,”
PLoS ONE , vol. 7, no. 7,p. e41010, 2012.[11] R. Singh, J. D. Quinn, P. M. Reed, and K. Keller, “Skill (or lack thereof)of data-model fusion techniques to provide an early warning signal for anapproaching tipping point,”
PLoS One , vol. 13, no. 2, p. e0191768, 2018.[12] S. R. Carpenter, D. Ludwig, and W. A. Brock, “Management of eutroph-ication for lakes subject to potentially irreversible change,”
Ecological Ap-plications , vol. 9, no. 3, pp. 751–771, 1999.1913] C. Boettiger and A. Hastings, “Early warning signals and the prosecutor’sfallacy,”
Proceedings of the Royal Society B: Biological Sciences , vol. 279,no. 1748, pp. 4734–4739, 2012.[14] G. D. Peterson, S. R. Carpenter, and W. A. Brock, “Uncertainty and themanagement of multistate ecosystems,”
Ecology , vol. 84, no. 6, pp. 1403–1411, 2003.[15] H. Moradkhani, S. Sorooshian, H. V. Gupta, and P. R. Houser, “Dualstate-parameter estimation of hydrological models using ensemble kalmanfilter,”
Advances in Water Resources , vol. 28, no. 2, pp. 135–147, 2005.[16] J. Liu and M. West, “Combined parameter and state estimation insimulation-based filtering,” in
Sequential Monte Carlo Methods in Prac-tice (A. Doucet, N. de Freitas, and N. Gordon, eds.), pp. 197–223, NewYork: Springer, 2001.[17] S. J. Lade and T. Gross, “Early warning signals for critical transitions:A generalized modeling approach,”
PLoS Computational Biology , vol. 8,no. 2, p. e1002360, 2012.[18] D. Piovani, J. Gruji ¨A ‡ , and H. J. Jensen, “Linear stability theory asan early warning sign for transitions in high dimensional complex sys-tems,” Journal of Physics A: Mathematical and Theoretical , vol. 49, no. 29,p. 295102, 2016.[19] J. Honerkamp,
Stochastische dynamische systeme . Weinheim: Wiley VCH,1990.[20] N. G. van Kampen,
Stochastic processes in physics and chemistry . Ams-terdam: North Holland, 2007.[21] R. Steuer, J. Kurths, O. Fiehn, and W. Weckwerth, “Observing and in-terpreting correlations in metabolomic networks,”
Bioinformatics , vol. 19,no. 8, pp. 1019–1026, 2003.[22] J. R. Magnus and H. Neudecker,
Matrix differential calculus with applica-tions in statistics and econometrics . Weinheim: Wiley, 2019.[23] G. Strang,
Linear algebra and its applications . Cambridge, MA: AcademicPress, 1976.[24] L. N. Trefethen and D. Bau III,
Numerical linear algebra . Philadelphia:SIAM, 1997.[25] A. Brechtel, P. Gramlich, D. Ritterskamp, B. Drossel, and T. Gross, “Mas-ter stability functions reveal diffusion-driven pattern formation in net-works,”
Physical Review E , vol. 97, no. 3, p. 032307, 2018.[26] A. Brechtel, T. Gross, and B. Drossel, “Far-ranging generalist top predatorsenhance the stability of meta-foodwebs,”
Scientific Reports , 2019.2027] U. Brose, R. J. Williams, and N. D. Martinez, “Allometric scaling enhancesstability in complex food webs,”
Ecology Letters , vol. 9, no. 11, pp. 1228–1236, 2006.[28] A. Nyberg, T. Gross, and K. E. Bassler, “Mesoscopic structures and thelaplacian spectra of random geometric graphs,”
Journal of Complex Net-works , vol. 3, no. 4, pp. 543–551, 2015.[29] H. Nakao and A. S. Mikhailov, “Turing patterns in network-organizedactivator-inhibitor systems,”
Nature Physics , vol. 6, no. 7, pp. 544–550,2010.[30] M. Baurmann, T. Gross, and U. Feudel, “Instabilities in spatially extendedpredator-prey systems: Spatio-temporal patterns in the neighborhood ofturing-hopf bifurcations,”
Journal of Theoretical Biology , vol. 245, no. 2,pp. 220–229, 2007.[31] T. Gross and U. Feudel, “Generalized models as a universal approach tothe analysis of nonlinear dynamical systems,”
Physical Review E , vol. 73,no. 1, p. 016205, 2006.[32] T. Gross, L. Rudolf, S. A. Levin, and U. Dieckmann, “Generalized modelsreveal stabilizing factors in food webs,”
Science , vol. 325, no. 5941, pp. 747–750, 2009.[33] S. J. Plitzko, B. Drossel, and C. Guill, “Complexity-stability relations ingeneralized food-web models with realistic parameters,”
Journal of Theo-retical Biology , vol. 306, pp. 7–14, 2012.[34] M. L. Rosenzweig and R. H. MacArthur, “Graphical representation and sta-bility conditions of predator prey interaction,”
American Naturalist , vol. 97,no. 895, pp. 209–223, 1963.[35] H. Aufderheide, L. Rudolf, T. Gross, and K. D. Lafferty, “How to predictcommunity responses to perturbations in the face of imperfect knowledgeand network complexity,”
Proceedings of the Royal Society B: BiologicalSciences , vol. 280, no. 1773, p. 20132355, 2013.[36] M. J. Khatibipour, F. Kurtoˇglu, and T. C¸ akir, “Jacly: a jacobian-basedmethod for the inference of metabolic interactions from the covariance ofsteady-state metabolome data,”
PeerJ , vol. 6, p. e6034, 2018.[37] P. Kloeden, E. Platen, and I. Wright, “The approximation of multiplestochastic integrals,”
Stochastic Analysis and Applications , vol. 10, no. 4,pp. 431–441, 1992. 21
Generation of noisy timeseries
The stochastic system is modeled as an It¯o process. The calculation is done byusing the Euler-Maruyama method [37] x ( t + dt ) = x ( t ) + ˙ x ( t ) d t + a (cid:112) x ( t ) d W ( t ) , (33)where d W ( t ) is the increment of a Wiener process with a normal distributionaround the mean 0 and the standard deviation √ d t . The amplitude of the noiseis proportional to the square root of the population size if we suppose that thenoise is due to intrinsic fluctuations of the population dynamics due to stochasticbirth and death processes. While realistic for ecological systems, it is importantto note that this does not match the stochastic behaviour (additive noise) thatwas assumed in derivation of the relationship between the co-variance and jaco-bian matrices. The subsequent results demonstrate that this assumption in theformulation of the method is not a barrier to its application. For the time serieswe used the noise strength a = 0 .
01 and different step sizes d t . The number ofused time steps and step sizes is listed in table 1. B Jacobian-Covariance relationship
We summarize the derivation of the Jacobian-Covariance relationship followingthe presentation in [21]. The response of the system to small fluctuations of thevariables around the equilibrium values can be approximated bydd t X = JX , (34)We can model the system with noise using a Langevin-type equationd X i d t = (cid:88) j J ij X j + (cid:112) D i ξ i ( t ) , (35)where ξ i ( t ) is Gaussian white noise, with zero mean and unit variance and D i is the mean amplitude of fluctuations.The corresponding stationary Fokker-Planck equation for the state proba-bility distribution P ( X ) is − (cid:88) ij J ij ∂∂X i X j P + (cid:88) ij D ij ∂ P∂X i ∂X j = 0 . (36)Multiplying Eq. (36) by X k X l and integrating gives (cid:88) j J kj (cid:104) X l X j (cid:105) + (cid:88) j J lj (cid:104) X k X j (cid:105) + 2 D kl = 0 , (37)where (cid:104) X i X j (cid:105) is the co-variance of variables X i and X j . Equation (37) can bewritten in matrix form JΓ + ΓJ T = − D , (38)where Γ is the co-variance matrix with entries Γ ij = (cid:104) X i X j (cid:105) [20].22 Vectorization of matrix products
Following [22] we consider a product of three matrices M = XYZ . We find thevectorization of M by stacking its columns, i.e.vec( M ) = m m ... (39)where m i is the i -th column of M . We can obtain the i -th column of M byreplacing Z by its i -th column z i , which yields XY z i = (cid:88) j X y j Z j,i = (cid:88) j Z j,i X y j (40)where y j is the j -th column of Y . The sum on the right-hand-side is also theproduct of the factors ( Z ,i X , Z ,i X , . . . ) = z i T ⊗ X (41)and y y ... = vec( Y ) (42)Hence XY z i = ( z i T ⊗ X )vec( Y ) (43)Stacking these equations for the different values of i yields XY z XY z ... = ( z T ⊗ X )vec( Y )( z T ⊗ X )vec( Y )... (44)and hence the result vec( XYZ ) = ( Z T ⊗ X )vec( Y ) (45) D Pseudoinverse of overdetermined system
Consider a system of the form B v = w (46)If A has a row-dimension that is greater than the column dimension of A this isan overdetermined system, so for a given A and w , we cannot generally expectthat there is a v that solves the equation. Hence for any v there will be someresidue x : B v − w = x (47)23ur aim is now to find the v such that the residue x is minimized. Specificallywe seek to minimize the euclidean norm | x | = x T x (48)This expression has a unique minimum at which the gradient vanishes [24].Hence we can find the desired v by demanding0 = ∇ ( x T x ) (49)where ∇ = ( ∂/∂v , ∂/∂v , . . . ) T . Transforming this equation we find0 = 2( ∇ x T ) x (50)= (cid:0) ∇ ( v T B T − w T ) (cid:1) x (51)= B T x (52)= B T B v − B T w (53)We can write this condition as B T B v = B T w (54)on the left-hand-side the square matrix B T B appears. In contrast to the rect-angular matrix B this matrix can typically inverted. Hence we can multiply theinverse from the left to obtain the desired formula v = (cid:0) B T B (cid:1) − B T w (55) E Metafoodweb model
Our aim in this section is to formulate a model that can serve as a test casefor the Jacobian reconstruction method. For this purpose we want a large,complex, and nonlinear dynamical network model where the Jacobian matrix isnevertheless analytically accessible,The model consists of a spatial network of N habitat patches linked by av-enues of species dispersal. Each patch harbours a complex food web, consistingof P nodes that represent populations of different species, which are linked bypredator-prey interactions.The variables of the model X si denotes the population density of species s onpatch i . For simplicity we assume that all patches are identical. The populationdynamics is given by the general form˙ X si = G ( X si ) + (cid:15)F ( X si , T si ) Local gains − (cid:16)(cid:80) p c sp X si T pi F ( X pi , T pi ) (cid:17) − M ( X si ) Local losses − d s k i X si + d s (cid:80) j A ij X sj Dispersal (56)Let’s unpack this equation. The first of the gain terms, G represents primaryproduction of biomass growth, e.g. by photosynthesis. We assume that this24erm is zero for all predator species. By contrast, the second term with thefunctional response F , which describes growth by predation, is assumed to bezero for primary producers. This second term represents growth by predation,which depends on the density of the predator, and the total density of prey, T si that fits the predators diet. The total prey density for predator s in patch i canbe written as T si = (cid:88) p c ps X pi (57)where c ps is the relative contribution that p makes to the diet of s . For example c ps = 1 if species p is easy prey for species s , but c ps = 0 if s cannot prey on p .The first loss term in Eq. (56) captures losses by predation, where we assumedthat the biomass uptake by a predator is assigned to the predator’s prey speciesaccording to their diet. The function M ( X si ) describes losses due to natural(i.e. non-predatory) mortality.The final two terms describe the effect of diffusive emigration from andimmigration to the patch, with a species dependent diffusion rate d s . Thematrix A is the adjacency matrix of the geographic network, and k i the degreeof patch i . With this diffusive coupling there is always a solution where allpatches in the networks have the same densities in all populations X si ∗ = X s ∗ (58)In the following we call this steady state the homogeneous state.We note that at this stage the functions G , F and M are still unspecified.Using the approach of Brechtel et al. [25] the Jacobian matrix for this typeof model can be computed analytically. The result is a Jacobian J that stillcontains a set unknown, but ecologically interpretable parameters, that describeproperties of the unspecified functions in the model. In [25] we showed that theJacobian in the homogeneous state can be written as J = I ⊗ P − L ⊗ K (59)Here I is an N × N identity matrix, P is the Jacobian of a single isolated patch, L is the Laplacian matrix of the geographical network and K is a diagonalcoupling matrix K = d d . . . (60)Following the procedure in [25] one can show that the eigenvalues λ of J canbe computed as λ = Ev ( P − κ K ) (61)where κ is a laplacian eigenvalue of the geographic network. Solving this equa-tion for every κ yields the complete set of eigenvalues of J .The patch Jacobian P for this model was already derived in [31] in a non-spatial context. In summary these results allow for the very convenient compu-tation of a broad class of fairly realistic meta-foodweb models. Once we have25ound the desired bifurcations in a food web it is straightforward to specify theunspecified function in the model to lie at a specific point of the generalizedmodel bifurcation diagram. In this way meta-foodwebs that exhibit the desiredtransitions can be generated very efficiently.For most of the analysis we use the functions G ( X s ) = α s X sφ s (62) F ( X s , T s ) = α s (1 + K s ) X sψ s T s T s + K s (63) M ( X s ) = α s X sµ s (64)here α s is the species biomass turnover rate, which we calculate from its trophicposition using an allometric scaling law; φ s is the elasticity of primary pro-duction, µ s is the elasticity of mortality, ψ s is the elasticity of predtion withrespect to predator density, and K s is a halfsaturation constant which we set to K s = γ s / (1 − γ s ) where γ s is the desired elasticity of predation with respect tothe prey.For the simulations with the Rosenzweig-MacArthur model we used˙ X i = rX i (cid:18) − X i c (cid:19) − aX i X i ahX i + (cid:88) j d ( X j − X i ) (65)˙ X i = (cid:15) aX i X i ahX i − mX i − ( wX i ) + (cid:88) j d ( X j − X i ) (66)The full set of parameters used in all simulations is shown in the tablesbelow.Figure 3 5 7 (a), (b) 7 (c), (d) 6 (a), (b) 6 (c), (d) T max
200 4 200 2 000 80 000 20 000 50 000steps 2 · . · · · · · ( ∗ )d t − − − − − − Table 1: The parameters of the solver. The Length of the simulation T max and the number of timesteps. ( ∗ ) In the case of figure 6 (c), (d) only 2 · ,that means every 25th step, of the time steps were used for the eigenvaluesestimation. 26arameter InterpretationExponent φ Sensitivity of primary production to own population density γ Sensitivity of predation to total available prey density ψ Sensitivity of predation to predator density µ Exponent of closureScale α Biomass flow σ Fraction of biomass loss due to predation˜ σ Fraction of biomass loss due to respiration β Relative contribution to biomass loss due to predation by acertain predator δ Fraction of local growth by predation˜ δ Fraction of local growth by primary production χ Relative contribution of population as prey to a certainpredatorTable 2: Generalized parameters used to describe the meta-foodweb. i α i φ i µ i ψ i γ i σ i δ i d i φ γ φ, γ ) start ( φ, γ ) end Bifurcation(0.25, 0.07) (0.30, 0.02) Transcritical(0.30, 0.10) (0.35, 0.05) Transcritical(0.35, 0.13) (0.40, 0.08) Transcritical(0.40, 0.16) (0.45, 0.11) Transcritical(0.45, 0.18) (0.50, 0.13) Transcritical(0.50, 0.21) (0.55, 0.16) Transcritical(0.55, 0.24) (0.60, 0.19) Turing(0.60, 0.27) (0.65, 0.22) Turing(0.65, 0.30) (0.70, 0.25) Turing(0.70, 0.35) (0.75, 0.30) Turing(0.75, 0.38) (0.80, 0.33) Turing(0.80, 0.41) (0.85, 0.36) Turing(0.85, 0.46) (0.90, 0.41) Turing(0.90, 0.51) (0.95, 0.46) Hopf(0.95, 0.56) (1.00, 0.51) HopfTable 4: Start and end points of the parameter trajectories used in the errorstatistics plot (Fig. 4) together with the corresponding bifurcation.27 c (cid:15) m w a h d d i α i φ i µ i ψ i γ i σ i δ i d i γ γ γ φ i, j , , , , χ ij . . β ij / / / / φ is used asthe bifurcation parameter, all other parameters are kept constant. For the therelative gain χ ij and loss β ij only the nonzero entries are shown. In the usedexample A ij = χ ij . If φ is used as the bifurcation parameter we choose γ = 0 . γ is used as the bifurcation parameter we choose φ = 0 ..