A closed form for Jacobian reconstruction from timeseries and its application as an early warning signal in network dynamics

The Jacobian matrix of a dynamical system describes its response to perturbations. Conversely one can estimate the Jacobian matrix by carefully monitoring how the system responds to environmental noise. Here we present a closed form analytical solution for the calculation of a system's Jacobian from a timeseries. Being able to access a system's Jacobian enables us to perform a broad range of mathematical analyses by which deeper insights into the system can be gained. Here we consider in particular the computation of the leading Jacobian eigenvalue as an early warning signal for critical transition. To illustrate this approach we apply it to ecological meta-foodweb models, which are strongly nonlinear dynamical multi-layer networks. Our analysis shows that accurate results can be obtained, although the data demand of the method is still high.


Introduction
As humans we are dependent on the functioning of complex systems on many scales, ranging from our own body with its interlinked metabolic, signalling and microbial networks, via technical and organizational systems such as power grids and political systems, to the planetary-scale supply chains and the climate systems. All of these are complex nonlinear many variable systems, and as such at a risk of undergoing sudden, qualitative and potentially irreversible transitions [1]. Over the past decades there has been a steadily growing interest in methods that provide early warning of such transitions [2,3]. In particular there is a growing awareness that such methods are needed for systems that are spatial or network based and hence inhherently high-dimensional [4,5,6].
The traditional approach to anticipating qualitative transitions in complex systems is mechanistic modelling. Many well-understood systems can be modelled so precisely that the model can accurately predict the threshold parameter values where transitions occur. For technical applications such as aircraft flight or stability of structures the model-based stability analysis is well established and in many cases part of a legally mandated licensing process [7]. However, model-based approaches tend to produce poor results in systems that are less well known as seemingly minor details of the model can sometimes drastically affect transition points.
It has long been known that critical transitions are generally preceded by critical slowing down [8,9]. This phenomenon leads to a distinctive increase in the auto-correlation and cross-correlations of timeseries before a transition. The advantage of correlation-based warning signs is that detailed understanding of the system is not needed. Its disadvantage is that high frequency time series data are necessary for robust warning, which imposes a strong, and for many applications prohibitive, constraint.
The model-based and correlation-based approaches to predicting critical transitions can be seen as two extreme strategies. The former is entirely based on structural knowledge of the system and uses real world data at most to fit model parameters. The latter eliminates the need for structural information, at the cost of a high time series data demand [10,11].
In many potential applications there are aspects of the system that are well understood because different variables are related via physical laws or subject to logical constraints [11]. If such bits of 'structural' knowledge are available in a system it is desirable to exploit them for the construction of early warning signals. However, the same systems may also contain aspects that are considerably less well understood and hence make purely model based predictions unreliable. This defines the need for a middle way, where available structural information on a system is used to reduce the data demand, while time series data is used to close gaps in the understanding in areas where structural information is not available [2].
A common approach to finding a middle way in the construction of early warning signs is to use data assimilation approaches to continuously improve a dynamical model [11]. A classical application of these approaches is the eutrophication of shallow lakes [12,13], for which good early warning signs can be constructed using techniques such as Bayesian learning and Kalman filters [14,15,16]. However, a broad investigation of such approaches shows that they may result in warning signs that are "faint and late" [11], while Boettiger and Hastings 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 general to encompass the whole class of conceivable models into which a given system could potentially fall. Using the so-called generalized modelling approach, the Jacobian matrices were computed that govern the system's response to pertur-bations and from which the bifurcation points can be computed. Because the underlying models contain unknown functional relationships the Jacobian matrices that are thus obtained still contain unknown parameters. The authors used time series data to eliminate these remaining uncertainties. They showed that in this way accurate early warning signals can be constructed that only require very limited timeseries data.
The present paper takes the idea of reconstructing a systems Jacobian from data one step further by avoiding assumptions on the functional form of the interaction terms. In a system that is subject to some noise, the cross correlations in timeseries encode very similar information to the Jacobian matrix of its deterministic backbone. Considered in isolation the cross correlations are not sufficient to reconstruct the full Jacobian, however a complete reconstruction is possible if we have some additional structural information. In systems that can be described as complex networks the network structure imposes constraints on which variables can interact directly, which in turn implies that some entries of the Jacobian must be zero. In a sufficiently sparse network, and particularly in multi-layer networks, knowing these zeroes provides sufficient information to reconstruct the remainder of the Jacobian from time series correlations.
For illustration we apply Jacobian reconstruction approach to an ecological meta-foodweb model, formulated as a dynamical multi-layer network. By comparing with the known ground truth of the model, i.e., its exact Jacobian, we show that the Jacobian can be reconstructed faithfully and demonstrate its value as an early warning signal. We find that despite leveraging the structural information, the amount of timeseries data required for accurate results is at present still prohibitively high for the ecological application. However we discuss several avenues of future research that may reduce the data requirements to a point where the method becomes widely applicable.

Mathematical Background
In every dynamical system that is in the vicinity of some form of long-term behaviour, the response of the system to small perturbations in the variables can be captured by some matrix. In the simplest case the system is a system of ordinary differential equations (ODEs) in which a vector of variables x evolves in time in a way that is dependent on a set of parameters p. The simplest form of long-term behavior is rest in a stationary state x * (p), such that The response to sufficiently small perturbations in the the steady state is then described by the Jacobian matrix J, whose elements are computed as The steady state under consideration is stable when all eigenvalues of the Jacobian have negative real parts. A smooth change in the parameters p will generally cause the the eigenvalues to change smoothly. When such a change causes one or more eigenvalues to acquire positive real parts, then a bifurcation occurs in which the dynamics change qualitatively. Because the Jacobian is a real matrix its eigenvalues are real or form complex conjugate pairs. Hence there are two fundamental types of bifurcations. In bifurcations of fold-type a single real eigenvalue acquires a positive real part as it passes through zero. Several different forms of this bifurcation are commonly encountered in dynamical systems (fold bifurcation, pitchfork bifurcation, transcritical bifurcation), but generally speaking these are associated with a change in the number of steady states in the system or the exchange of stability properties.
In a Hopf bifurcation a complex conjugate eigenvalue pair acquires a positive real part by crossing the imaginary axis of the complex plane. This bifurcation is generally associated with the onset of at least transient oscillations as the system departs from the steady state.
Both types of bifurcations can occur in several different forms, some of which cause only relatively mild non-catastrophic transitions (say, replacing stationary behaviour 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 computing normal 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 on a spatial network of discrete nodes, the fundamental bifurcations can come in two flavors: A bifurcation may affect all points in space simultaneously in the same way, or it may affect points in space differently leading to the formation of spatial patterns. For clarity the bifurcations of the latter type are called Turing bifurcations (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 cause for concern. In this spirit our aim in the following is to reconstruct the Jacobian of a dynamical system from data in order to determine its leading eigenvalue. If the real part of this eigenvalue approaches zero, we interpret this as a warning signal for an impending bifurcation.

Jacobian Reconstruction
Our aim in this section is to formulate a method that can reconstruct the Jacobian matrix of a dynamical system from time series. We do this by expanding on 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 the covariance matrix Γ of timeseries, and to the fluctuation matrix D of the noise, via the equation In the following, we show that we can use this relationship to compute J. Consider that a Jacobian matrix of linear dimension N contains N 2 independent entries. By contrast, Eq. (4) equates two symmetric matrices, and hence imposes only N (N + 1)/2 constraints on the elements of J. Therefore in any application with multiple variables (N > 1) the system is underdetermined such that we cannot recover the complete Jacobian purely from the time series. We can still recover the full Jacobian if we have additional information that we can leverage. Fortunately, in many applications some structural information is easily accessible. In particular in large spatially complex systems or reasonably sparse networks we know that certain variables cannot interact and hence the corresponding elements of the Jacobian must be zero. This yields a set of structural constraints of the form J ij = 0 for certain pairs (i, j). If we can identify at least N (N − 1)/2 such zero entries, then the timeseries contain enough information to reconstruct all remaining Jacobian entries. Below in Sec. 6 we demonstrate that this is quite generally the case.
Given a timeseries of the system's N variables and a set of G additional structural constraints, with G ≥ N (N − 1)/2 we now solve the Eq. 4 to estimate the nonzero entries of the Jacobian. To understand how this equation is solved let us first consider the two-dimensional example which implies the independent conditions 2J 11 Γ 11 + 2J 12 Γ 12 = −2D 11 (5) with the second condition applying to both off-diagonal terms. The left-hand side of these equations is a linear system. Hence we can write the conditions in the form where B is a matrix, j is a column vector that contains the entries of the Jacobian, i.e. j = (J 11 , J 21 , J 12 , J 22 ) T and d is the corresponding vector for D. For the two dimensional example this reads The form of this equation suggests that we can solve it for j by multiplying B −1 from the left. However, we have to take care because B is not invertible because the two center rows are identical, which is a consequence of the missing information. We can fix this problem by using the additional constraints. Imposing structural constraints on the two-variable system makes this example almost pointless, but, for the purpose of illustration, let us assume that we know that variable 1 cannot depend on variable 2 and hence J 12 = 0. We can represent this constraint in the same matrix equation as the system by adding it as an additional line, If this is the only constraint then we can now drop the third row of the matrix and the third entry of the vector on the right-hand side and solve for j by matrix inversion. In practice there are typically additional constraints and hence the system will be overdetermined. In this case we use least squares optimization to find an approximate solution. As we will see below the form of Eq. (10) is very convenient for finding the least squares solution. In particular it allows us to obtain a analytical closed form solution for j.
Let us now generalize from the two-dimensional example to systems with an arbitrary number of variables. For this purpose we define the vectorization operator [22] vec(X) = (X 11 , X 21 , . . . X N 1 , X 12, . . .) T We can now write Eq. (8) as i.e. the vectorization of a matrix is the columns of the matrix stacked on top of each other. To find the general form of B we start from Eq. (4) and vectorize both sides, which yields Because vectorization is a linear operator we can pull the -2 out of the vectorization on the right hand side and apply the vectorization separately to the two terms on the left hand side, hence It is known [22] that for matrices X, Y, Z the following identity holds (cf. appendix C): where ⊗ is the Kronecker product of matrices defined by We now substitute Y = J, Z = Γ and X = I, where I is the identity matrix of appropriate size. This yields which brings the first term from Eq. (14) into the desired form. If we try the same for the second term we find 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 matrix and the vectorization of its transpose are not identical but closely related.
Consider that for our two-dimensional the two vectorizations are related by where is a permutation matrix. In the general case we can still write where C is now a permutation matrix of size N 2 × N 2 , one can construct this matrix from blocks of size N × N as Here, C nm is a N × N -matrix defined by where δ is the Kronecker delta.
Using the commutation matrix we can write Eq. (18) as Substituting this relationship and Eq. (17) into Eq. (14) we get the desired form We now know how to construct the matrix B such that we can write the system in the form 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 purpose we can gather the respective elements of j in an ordered list U and then define U as an |U| × N matrix with Each row of this matrix represents one of the constraints that we wish to impose. For example if in the fourth row the only nonzero entry is in the 8th column this means our forth condition is that the 8th element of j must be zero. In analogy to the small example we can now impose the additional conditions on the system by stacking them below the matrix B such that Eq. (26) becomes where 0 is a column vector containing |U| zeroes. The () that appear in this equation should be read as a block-wise notation for matrices/vectors, where rows are stacked on top of each other in analogy to Eq. (10).
To simplify the notation we introducê which allows us to write the whole set of conditions once again in the form In practise this will now be an overdetermined system, such that no exact solution exists. However, finding an approximation that minimizes the squares of the deviations in each row is a well-known problem. The known solution [23,24] (cf. appendix D) for this problem is where the expression (B TB ) −1BT , the pseudoinverse ofB,appears.
The equation provides a closed form solution by which the Jacobian elements can be computed from the covariance matrix of the timeseries, a known or estimated fluctuation matrix, and a set of additional structural constraints on the Jacobian.
Typically the least squares fit will not set the Jacobian elements governed by the structural constraints exactly to zero. One may therefore enforce these known zeros by setting the respective elements of the Jacobian explicitly to zero after the computation of Eq. (31) has finished. In numerical experiments, described below, we found that this improved the accuracy of Jacobian eigenvalues estimated by this method.
To summarize, the Jacobian J of an N -dimensional system close to a steady state can be reconstructed as follows: 1. Compute the co-variance matrix from the time series data, Γ ij = X i X j .
2. Construct the diagonal fluctuations matrix D, and compute d = vec(D).
In some systems these fluctuations can be measured directly, otherwise a reasonable approximation may be derived based on assumptions on the underlying noise process [20].

Compute the matrix
where I is the N × N identity matrix.
5. Define j = vec(J) and identify at least N (N −1)/2 elements of j that must be 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, and row 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 = 1 if the nth structural condition reads j m = 0.
and recover the Jacobian J from its vectorization j.
8. Set the elements of J governed by structural constraints explicitly to zero.

Application to a meta-foodweb model
In the following we explore if the Jacobian matrix can be reconstructed sufficiently accurately to warn of impending critical transitions. For this purpose, we constructed a test system of realistic complexity for which the Jacobian matrix is nevertheless known analytically, and we created noise timeseries. We used a meta-foodweb model already studied in [25,26] (see appendix E for details). The model consists of a spatial network of P habitat patches linked by avenues of species dispersal (Fig. 1). Each patch harbours a complex foodweb, consisting of S nodes that represent populations of different species, which are linked by predator-prey interactions. The dynamics of the system are given by a set of differential equations that govern the changes in variables due to diffusive dispersal between patches and biological processes occurring within a 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 availability of different food sources (prey switching). The dynamics of different species occurs on different time scales according to biological scaling relationships which relate a species position in the foodweb to its expected biomass turnover rate [27]. Similar scaling relationships also govern the rate at which different species disperse across the spatial network [25].
We consider two different versions of the patch topology. The smaller of the two consists 6 patches which are connected in such a way that they form the smallest completely asymmetric network (Fig. 1D). The larger one contains 15 patches and was generated as a random geometric graph (Fig. 1E). It has a comparatively large diameter and high clustering and contains several symmetries that are characteristic for this type of spatial networks [28]. For the food webs we used a predator-prey system consisting of two species (Fig. 1C), as well as a 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 approach to meta-foodwebs [25] showed that network models behave analogously to dynamical systems in continuous space. Hence the theory of pattern formation in partial-differential equation systems can transferred almost exactly to these dynamical networks. This means that our test systems can exhibit instabilities that are best described as pattern-forming bifurcations, specifically Turing and Wave instabilities. Pattern forming instabilities in predator-prey systems in continuous space were studied in detail in [30] which made it easy to locate these bifurcations also in our multi-layer network predator-prey system (Fig. 2).
Using a generalized modelling approach [31,32,33] we analytically computed the Jacobian that describes the Jacobian matrices of the two example food webs on arbitrary spatial topologies. We then picked specific realizations of models in which relevant bifurcations occurred. To produce noisy timeseries we simulated these models with added noise using the Euler-Maruyama method (see appendix A). Figure 1: Schematic representation of the model system. We consider a multilayer network where copies of an ecological food web exist in different geographical patches (A). We consider an intra-guild predation food web with an additional predator (B) and a predator-prey system (C). For the spacial network, we use the smallest completely asymmetric graph (D) and a larger random topology, generated as a random geometric graph(E).

Results
As a first test we generated simulated noisy timeseries containing 2 · 10 5 data points from the two-species system on the six-patch topology. Parameter values for these simulations were chosen to lie on four transects through the parameters space that crossed bifurcation points. For each of these time series we then reconstructed the Jacobian matrix and computed the leading eigenvalue. Before the bifurcation point the estimated eigenvalues are in very good agreement with the known ground truth provided by the analytic eigenvalues (Fig. 3).
As we follow the transects the real part of the leading eigenvalue crosses zero in a bifurcation. When the bifurcation occurs the reconstructed eigenvalue departs from the analytical value. This behavior is expected as the analytical solution continues to show the eigenvalues around the, now unstable, steady . The bifurcation diagram was computed using the master stability function approach from [25] (see appendix). It corresponds directly to Fig. 1 from [30] which studies a predator-prey system in continuous space. Lines (a-d) indicate the transects used for the corresponding simulations in Fig. 3.
state, whereas the reconstruction algorithm computes the leading eigenvalue associated with the new dynamics, which has now departed from the previous steady state. We note that in transect d) the reconstructed eigenvalue is almost exactly zero in a wide region after the bifurcation. This happens because the system approaches a stable limit cycle for which the leading Lyapunov exponent is zero. The recovery of this zero by the algorithm provides some (unexpected) evidence suggesting that the method reveals some salient information even in non-stationary states.
Careful examination of the transects (Fig. 3) suggests that the accuracy of reconstruction improves as we approach the bifurcation point. To explore this further we generated 15 transects in the vicinity of bifurcation points and generated three sets of timeseries along every one of the transects (see appendix for details). The results confirm that the accuracy of the estimate improves as the bifurcation point is approached (Fig. 4).
We now consider the case where parameters are slowly changing over a long simulation run. Our aim is to estimate the leading eigenvalue of the Jacobian over time as this slow change in the system is taking place. For this purpose we apply the proposed method to reconstruct the Jacobian in sliding time win- dow of length τ . The results (Fig. 5) show that the estimates based on the sliding window are more noisy, undergoing visible fluctuations around the true value. Nevertheless the trend of the eigenvalue approaching zero is still clearly captured.
To test the applicability to larger networks, we apply Jacobian reconstruction to a model system with 4 species on 15 patches. Even in this larger system the accuracy is still quite good but the estimated eigenvalue is systematically slightly less than the true value. For the case of fixed parameter values the difference once again disappears as we approach the bifurcation. But for the case of sliding parameter values a small difference remains. We suspect that this may be the combined effect of the non-autonomous nature of this systems and the noise leading to bifurcation delay. This delay effect can be expected to be more pronounced in the larger food web due to the presence of higherlevel predator whose dynamics happen on correspondingly longer timescales. If this is the case then the reconstructed Jacobian eigenvalue may actually offer a better estimate of the relevant transition point than the analytical solution for the system without noise. So far we studied systems that were designed using the generalized modelling approach. We complement this by the analysis of a well established ecological model, the Rosenzweig-MacArthur predator-prey model [34] with quadratic mortality and diffusion on the 6-patch network (see appendix for details). The results of Jacobian reconstruction (Fig. 7) show that the leading eigenvalue can be recovered with reasonable accuracy, again the accuracy of the estimate improves significantly as the system approaches the bifurcation point.

Summary and discussion
In this paper we expanded on previous work by Honerkamp, van Kampen, Steuer and others to derive a closed form expression for the reconstruction of Jacobian matrices from time series data.
For illustration we applied the mathematical formula to the an ecological meta-foodweb model. This example illustrated that a relatively robust reconstruction of the leading eigenvalue of the Jacobian is possible even in a strongly nonlinear multi-layer network with dynamics on multiple timescales. However, the example also revealed that the required amount of data is still prohibitive Figure 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 parameters, governing the nonlinearity of primary production, φ (panel a,c) and the nonlinearity of the functional response of attack rate to prey density, γ (b,d). The leading eigenvalue of the Jacobian was reconstructed from timeseries for fixed parameters (a,b) and parameter transsects (c,d). The steady state under consideration is unstable subsequent to a Hopf bifurcation (region shaded grey). For the transects estimates are shown in the center of the sampling window (window width is indicated by blue area). See appendix for parameters and details.
for the ecological application.
There is reason to believe that future research and particular a deeper mathematical understanding of the Jacobian reconstruction can significantly reduce the required amount of data. Particularly interesting in this respect is the ob- Figure 7: Eigenvalue estimations in the Rosenzweig-MacArthur with quadratic mortality and diffusion on the six-patch network. The steady state under consideration loses stability due to a Fold (panels a,c) and Hopf (b,d) bifurcations. The leading eigenvalue was reconstructed from simulation runs with fixed parameters (a,b) and slowly changing parameters (c,d). The accuracy of the estimates improves as the system approaches the bifurcation. (Estimated values are shown in the middle of the sampling window, indicated in blue. See appendix for details) served increase in accuracy close to the bifurcation point. A promising goal for future exploration would be to understand how far this region of heightened accuracy extends. If a real world application under consideration is already close to bifurcation one might find that much less data is required to reach the desired accuracy.
Alternatively, it might be possible to reduce the data demand by optimizing the 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 Jacobian. Using the initial estimate one could then identify the relevant timescales and important entries in the covariance matrix and optimize the sampling effort accordingly.
A third alternative may be to use additional information that may be available. The advantage of our jacobian-based approach is that it can take advantage of additional knowledge on the systems that may be available in ecological applications, such as closure exponents or the predator-dependence of the predation rate. The previous example [17] suggests that the use of this information may reduce the data demand considerably.
Meanwhile the method proposed here may be useful in fields where data is more readily available, such as studies of metabolism, power grids, or economic data. In the study of metabolism Jacobian reconstruction is already frequently used [36], for this application the present work yields an analytic closed form solution to a problem that is so far solved by machine learning methods. For power grids, reconstructing Jacobians may be particularly interesting because it could yield deeper insights into the functioning of the system in addition to providing an early warning signal. For economics it may be particularly interesting that the Jacobian can be seen as a representation of causality in the system. Closed form jacobian reconstruction thus offers a way to infer causality from correlation and is guaranteed to be exact in the large data limit. An important caveat is that the method requires some additional information to avoid underdeterminedness. However, as illustrated here, it is already sufficient if we can identify N (N − 1)/2 variable pairs that do not interact directly among the N variables of the system. In networks language, this is equivalent to saying that the mean degree z of the network of interactions must obey a condition that should be easy to satisfy in many applications and becomes easier to meet in larger systems. In summary, we find that Jacobian reconstruction is a promising approach to the analysis of complex systems near critical transitions, although the data requirements presently still limit its applicability. We expect that the closed form solution derived in the present paper inspires future mathematical work to alleviate these requirements.

A Generation of noisy timeseries
The stochastic system is modeled as an Itō process. The calculation is done by using the Euler-Maruyama method [37] x(t + dt) = x(t) +ẋ(t) dt + a x(t) dW (t) , where dW (t) is the increment of a Wiener process with a normal distribution around the mean 0 and the standard deviation √ dt. The amplitude of the noise is proportional to the square root of the population size if we suppose that the noise is due to intrinsic fluctuations of the population dynamics due to stochastic birth and death processes. While realistic for ecological systems, it is important to note that this does not match the stochastic behaviour (additive noise) that was assumed in derivation of the relationship between the co-variance and jacobian matrices. The subsequent results demonstrate that this assumption in the formulation of the method is not a barrier to its application. For the time series we used the noise strength a = 0.01 and different step sizes dt. The number of used time steps and step sizes is listed in table 1.

B Jacobian-Covariance relationship
We summarize the derivation of the Jacobian-Covariance relationship following the presentation in [21]. The response of the system to small fluctuations of the variables around the equilibrium values can be approximated by We can model the system with noise using a Langevin-type equation 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 probability distribution P (X) is Multiplying Eq. (36) by X k X l and integrating gives where X i X j is the co-variance of variables X i and X j . Equation (37) can be written in matrix form where Γ is the co-variance matrix with entries Γ ij = X i X j [20].

C Vectorization of matrix products
Following [22] we consider a product of three matrices M = XYZ. We find the vectorization of M by stacking its columns, i.e.
where m i is the i-th column of M. We can obtain the i-th column of M by replacing Z by its i-th column z i , which yields where y j is the j-th column of Y. The sum on the right-hand-side is also the product of the factors and    Hence Stacking these equations for the different values of i yields and hence the result D Pseudoinverse of overdetermined system Consider a system of the form Bv = w If A has a row-dimension that is greater than the column dimension of A this is an overdetermined system, so for a given A and w, we cannot generally expect that there is a v that solves the equation. Hence for any v there will be some residue x: Our aim is now to find the v such that the residue x is minimized. Specifically we seek to minimize the euclidean norm This expression has a unique minimum at which the gradient vanishes [24]. Hence we can find the desired v by demanding where ∇ = (∂/∂v 1 , ∂/∂v 2 , . . .) T . Transforming this equation we find We can write this condition as on the left-hand-side the square matrix B T B appears. In contrast to the rectangular matrix B this matrix can typically inverted. Hence we can multiply the inverse from the left to obtain the desired formula

E Metafoodweb model
Our aim in this section is to formulate a model that can serve as a test case for the Jacobian reconstruction method. For this purpose we want a large, complex, and nonlinear dynamical network model where the Jacobian matrix is nevertheless analytically accessible, The model consists of a spatial network of N habitat patches linked by avenues of species dispersal. Each patch harbours a complex food web, consisting of P nodes that represent populations of different species, which are linked by predator-prey interactions.
The variables of the model X si denotes the population density of species s on patch i. For simplicity we assume that all patches are identical. The population dynamics is given by the general forṁ Let's unpack this equation. The first of the gain terms, G represents primary production of biomass growth, e.g. by photosynthesis. We assume that this term is zero for all predator species. By contrast, the second term with the functional response F , which describes growth by predation, is assumed to be zero 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 can be written as 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 assumed that the biomass uptake by a predator is assigned to the predator's prey species according 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 and immigration to the patch, with a species dependent diffusion rate d s . The matrix A is the adjacency matrix of the geographic network, and k i the degree of patch i. With this diffusive coupling there is always a solution where all patches in the networks have the same densities in all populations 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 type of model can be computed analytically. The result is a Jacobian J that still contains a set unknown, but ecologically interpretable parameters, that describe properties of the unspecified functions in the model. In [25] we showed that the Jacobian in the homogeneous state can be written as 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 diagonal coupling matrix Following the procedure in [25] one can show that the eigenvalues λ of J can be computed as where κ is a laplacian eigenvalue of the geographic network. Solving this equation for every κ yields the complete set of eigenvalues of J.
The patch Jacobian P for this model was already derived in [31] in a nonspatial context. In summary these results allow for the very convenient computation of a broad class of fairly realistic meta-foodweb models. Once we have found the desired bifurcations in a food web it is straightforward to specify the unspecified function in the model to lie at a specific point of the generalized model bifurcation diagram. In this way meta-foodwebs that exhibit the desired transitions 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 M (X s ) = α s X s µs (64) here α s is the species biomass turnover rate, which we calculate from its trophic position using an allometric scaling law; φ s is the elasticity of primary production, µ s is the elasticity of mortality, ψ s is the elasticity of predtion with respect 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 to the prey. For the simulations with the Rosenzweig-MacArthur model we useḋ The full set of parameters used in all simulations is shown in the tables below. 80 000 20 000 50 000 steps 2 · 10 5 4.2 · 10 6 2 · 10 5 2 · 10 6 2 · 10 6 5 · 10 7 ( * ) dt 10 −3 10 −3 10 −2 10 −2 10 −2 10 −3 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 · 10 6 , that means every 25th step, of the time steps were used for the eigenvalues estimation.
Parameter Interpretation Exponent φ Sensitivity of primary production to own population density γ Sensitivity of predation to total available prey density ψ Sensitivity of predation to predator density µ Exponent of closure Scale α Biomass flow σ Fraction of biomass loss due to predatioñ σ Fraction of biomass loss due to respiration β Relative contribution to biomass loss due to predation by a certain predator δ Fraction of local growth by predatioñ δ Fraction of local growth by primary production χ Relative contribution of population as prey to a certain predator Table 2: Generalized parameters used to describe the meta-foodweb. Table 3: Parameters used for the two species predator-prey system.  Table 4: Start and end points of the parameter trajectories used in the error statistics plot (Fig. 4) together with the corresponding bifurcation.    Table 6: Parameters used for the 4 species food web system. φ 4 is used as the bifurcation parameter, all other parameters are kept constant. For the the relative gain χ ij and loss β ij only the nonzero entries are shown. In the used example A ij = χ ij . If φ is used as the bifurcation parameter we choose γ = 0.95. If instead γ is used as the bifurcation parameter we choose φ = 0.93.