# Biological meaning of Lotka-Volterra Jacobian matrices

We are searching data for your request:

Forums and discussions:
Manuals and reference books:
Data from registers:
Wait the end of the search in all databases.
Upon completion, a link will appear to access the found materials.

In a book about community ecology, I learned about Lotka-Volterra models and dynamic food web models (I am not an ecologist in the first place).

In one of the chapters, a Jacobian matrix of partial derivatives of the Lotka-Volterra model at equilibrium densities is given: egin{bmatrix}-&-&0&0+&0&-&0&+&0&-&0&+&0end{bmatrix} where + and - represent positive or negative values.

This matrix is later interpreted in terms of direction of trophic links in a food web, but the reasoning allowing to jump form one to the other is not detailed.
So what is the biological meaning of these values?

The Jacobian tells how the system changes along different state variables (which can be, for instance, the concentrations of the predator and the prey).

The Jacobian matrix by itself doesn't give you a lot of intuitive information. However, the eigenvalues of the Jacobian matrix at the equilibrium point tell you the nature of the steady state. For example if all the eigenvalues are negative then the system is stable at that point and the equilibrium point is called a stable node. If some of them are positive and some are negative then the system is metastable at that point and the point is called a saddle point. If all eigenvalues are positive then the point is unstable. If the eigenvalues have an imaginary part then the system exhibits oscillations. All these properties are local i.e. they just tell you the behaviour of the system near the equilibrium point.

I hope I am clear. This is not in a strict sense, a biological question. If you want to know more then read a book on nonlinear dynamics. Non-linear dynamics and Chaos by Steven Strogatz is a good place to start.

EDIT

When the system is locally linear (linearized at a given point) then the system dynamics at that point can be described as \$\$frac{dx}{dt}=J.x\$\$

Where \$x\$ is the state vector (a vector of the variables in the system) and \$J\$ is the Jacobian matrix. In this situation, the signs will denote the effect of a variable on the other variable. However, I am not sure (I feel it is highly unlikely) if this can be used to make any inferences about the food web.

The matrix above depicts the signs of the partial change of a function with respect to each species. Because there are 4 rows and columns, this means there are 4 species. Each row represents a partial change in each of the 4 differential equations representing the rates of change for each species with respect to each species, which is each column. To rewrite the Jacobian for species 1-4 ($$N_1, N_2, N_3, N_4)$$: $$mathbf{J} = left[egin{smallmatrix} frac{partial F_1}{partial N_1} & frac{partial F_1}{partial N_2} & frac{partial F_1}{partial N_3} &frac{partial F_1}{partial N_4} frac{partial F_2}{partial N_1} & frac{partial F_2}{partial N_2} & frac{partial F_2}{partial N_3} &frac{partial F_2}{partial N_4} frac{partial F_3}{partial N_1} & frac{partial F_3}{partial N_2} & frac{partial F_3}{partial N_3} &frac{partial F_3}{partial N_4} frac{partial F_4}{partial N_1} & frac{partial F_4}{partial N_2} & frac{partial F_4}{partial N_3} &frac{partial F_4}{partial N_4} end{smallmatrix} ight]$$ The signs above correspond the effect of species on each other's rates of change. For a first example, the first row and column is effect of species 1 on its own rate of change. It's negative, meaning that it's rate of change decreases with increasing density (crowding). A second example is the first row and second column-species 2 has a negative effect on the growth of species 1. In a food web context, this means that species 2 consumes species 1. A third example is the first row and third column-there is no effect of species 3 on species 1. Lastly, the second row and first column has a positive sign, meaning that species 1 has a positive effect on the growth rate of species two. In a food web context, that means that species 2 consumes species 1.

Notice that in food webs, for each positive effect, there is a corresponding negative effect and vice versa (i.e., species grow from eating another species, and species shrink from being eaten).

In sum, (1) species 1 is the only species who is affecting by crowding, (2) species 1 is consumed by species 2, (3) species 2 is consumed by species 3, (4) species 3 is consumed by species 4, and (5) no other species interact.

## Generalized Lotka–Volterra equation

The generalized Lotka–Volterra equations are a set of equations which are more general than either the competitive or predator–prey examples of Lotka–Volterra types. [1] [2] They can be used to model direct competition and trophic relationships between an arbitrary number of species. Their dynamics can be analysed analytically to some extent. This makes them useful as a theoretical tool for modeling food webs. However, they lack features of other ecological models such as predator preference and nonlinear functional responses, and they cannot be used to model mutualism without allowing indefinite population growth.

## Integrated Population Biology and Modeling, Part B

Paul Georgescu , . Hong Zhang , in Handbook of Statistics , 2019

### 3.3 An Unlikely Marriage: Two-Sex Reproduction in the Mutualistic Framework

Population models that separate the individuals by sex have a long history in the mathematical biology literature. All models, whether purely demographic or extended to larger ones as in the case of epidemic models, come with increased mathematical challenges. First, the number of equations doubles which makes the model far more difficult to analyze. Second, as shown below, the portion of the model that describes the two-sex reproduction/mating is usually not differentiable at the origin which makes standard stability analysis techniques unusable when trying to establish population extinction/persistence threshold conditions.

These models evolved into their own branch of mathematical biology, i.e., they are not usually presented as cases of mutualism for reasons that will be discussed shortly. However, in some cases, they still fit the mutualistic framework that is the subject of this chapter. Consider the following two-sex model in which F and M denote the female and male populations at time t (for a comprehensive survey of two-sex modeling see Iannelli et al. (2005) and references within):

The meaning of the parameters are: μf and μm are the mortality rates for females and males respectively, γf and γm are the probabilities of a newborn being a female or male respectively (hence γf + γm = 1), β is the per-mating birth rate and, finally, the central component of the model, M ( F , M ) is the number of matings per unit of time. It is this last term than can be interpreted as a mutualistic interaction between F and M. While demographers do not completely agree on what specific form M should take, there is, in general, agreement on a set of properties that this function should satisfy:

positivity: M ( F , M ) ≥ 0 whenever F ≥ 0 and M ≥ 0.

heterosexuality: M ( F , 0 ) = M ( 0 , M ) = 0 . This indicates that the number of matings is zero if individuals of one gender are absent. In the language of mutualistic interactions this would correspond to an obligatory mutualism.

degree one homogeneity: M ( α D , α M ) = α M ( F , M ) for every α > 0. This ensures the preservation of sex ratio if the number of singles in each gender changes by a factor α.

monotonicity: M ( F , M ) is increasing in F and M. This means that the number of matings increases with the number of individuals of each gender which, obviously, corresponds to the mutualistic assumption in the two-species models.

consistency: M ( F , M ) ≤ k F and M ( F , M ) ≤ k M for some constant k > 0 and for every F and M. This ensures that the number of matings is limited by the number of available individuals of a certain gender. Put it differently, this is a scarcity argument: if one gender is over-represented, the number of matings is controlled by the opposite (and scarce) gender. In the context of mutualistic models this correspond to a saturation effect in which the benefit from the other species approaches a certain limit if that species size increases to infinity.

Since all demographic models designed for long-term population modeling must contain some form of logistic assumptions on the natural mortality rate, we see that these models will not be suitable to the two-species mutualistic framework. Nevertheless, two-sex models can also be used for short-term modeling in which case the logistic effects on the mortality may be neglected. In such a case, the models do fit the mutualistic framework. Typical examples of the mating function analyzed in the literature are: the harmonic mean F M F + M , geometric mean F M or the generalized weighted mean a F b + ( 1 − a ) M b 1 / b .

One important additional reason why two-sex models are not typically treated as particular cases of mutualism is that, in many cases, the role of individuals of each gender does not appear only as a positive contribution to the other. Indeed, in the case of using the linear logistic mortality, the model ( 12 ) becomes

In this case, F and M, in addition to their positive contribution in the reproduction terms, they also increase the mortality due to logistic effect. As such, model ( 13 ) does not fit the mutualistic interaction framework discussed in this chapter.

## 3 Main results

Theorem 3.1 Assume that α < 1 and δ < 1 , then the following statements are true.

The equilibrium point O = ( 0 , 0 ) is locally asymptotically stable.

The equilibrium point Q = ( − 1 + α γ , 0 ) is unstable.

The equilibrium point R = ( 0 , − 1 + δ η ) is unstable.

Proof (i) The Jacobian matrix of the linearized system of (1.1) about the fixed point ( 0 , 0 ) is given by

Moreover, the eigenvalues of the Jacobian matrix J F ( 0 , 0 ) about ( 0 , 0 ) are λ 1 = α < 1 and λ 2 = δ < 1 . Hence, the equilibrium point ( 0 , 0 ) is locally asymptotically stable.

The Jacobian matrix of the linearized system of (1.1) about the fixed point ( − 1 + α γ , 0 ) is given by

The eigenvalues of the Jacobian matrix J F ( − 1 + α γ , 0 ) about ( − 1 + α γ , 0 ) are λ 1 = 1 α > 1 and λ 2 = γ δ − ϵ + α ϵ γ .

The Jacobian matrix of the linearized system of (1.1) about the fixed point ( 0 , − 1 + δ η ) is given by

The eigenvalues of the Jacobian matrix J F ( 0 , − 1 + δ η ) about ( 0 , − 1 + δ η ) are λ 1 = 1 δ > 1 and λ 2 = β − β δ + α η η . □

Theorem 3.2 The following statements are true.

If α > 1 , δ < 1 , and ϵ < γ − γ δ α − 1 , then the equilibrium point Q = ( − 1 + α γ , 0 ) is locally asymptotically stable.

If δ > 1 and α < 1 , then the equilibrium point R = ( 0 , − 1 + δ η ) is locally asymptotically stable.

Theorem 3.3 Assume that α > 1 , δ > 1 , and η > − β + β δ − 1 + α , then the unique equilibrium point P = ( β − β δ + ( − 1 + α ) η β ϵ + γ η , γ ( − 1 + δ ) + ( − 1 + α ) ϵ β ϵ + γ η ) is locally asymptotically stable if

Proof Assume that α > 1 , δ > 1 , and η > − β + β δ − 1 + α . Let L = β − β δ + ( − 1 + α ) η > 0 . Then a characteristic polynomial of the Jacobian matrix F J ( P ) about the unique equilibrium point P = ( L β ϵ + γ η , γ ( − 1 + δ ) + ( − 1 + α ) ϵ β ϵ + γ η ) is given by

Assume that Ω < ( β ( γ − γ δ + ϵ ) + α γ η ) 2 ( γ δ η + ϵ ( β + ( − 1 + α ) η ) ) . Then one has

Then, by Rouche’s theorem, S ( λ ) and S ( λ ) − T ( λ ) have the same number of zeroes in an open unit disk | λ | < 1 . Hence, the unique positive equilibrium point P is locally asymptotically stable. □

### 3.1 Global character

Theorem 3.4 Let I = [ a , b ] and J = [ c , d ] be real intervals, and let f : I × J → I and g : I × J → J be continuous functions. Consider the system (2.1) with initial conditions ( x 0 , y 0 ) ∈ I × J . Suppose that the following statements are true.

f ( x , y ) is non-decreasing in x and non-increasing in y.

g ( x , y ) is non-decreasing in both arguments.

If ( m 1 , M 1 , m 2 , M 2 ) ∈ I 2 × J 2 is a solution of the system

such that m 1 = M 1 and m 2 = M 2 , then there exists exactly one equilibrium point ( x ¯ , y ¯ ) of the system (2.1) such that lim n → ∞ ( x n , y n ) = ( x ¯ , y ¯ ) .

Proof According to the Brouwer fixed point theorem, the function F : I × J → I × J defined by F ( x , y ) = F ( f ( x , y ) , g ( x , y ) ) has a fixed point ( x ¯ , y ¯ ) , which is a fixed point of the system (2.1).

## Integrating the ODE using scipy.integrate

Now we will use the scipy.integrate module to integrate the ODEs. This module offers a method named odeint, which is very easy to use to integrate ODEs:

infodict is optional, and you can omit the full_output argument if you don’t want it. Type “info(odeint)” if you want more information about odeint inputs and outputs.

We can now use Matplotlib to plot the evolution of both populations:

The populations are indeed periodic, and their period is close to the value T_f1 that we computed.

## 2 Methods

### 2.1 Empirical food webs

Our local stability analysis was performed on Jacobian community matrix models (linearizations of systems of differential equations) of the same food-web structures as were used for the feedback analysis in Neutel and Thorne ( 2014 ). These were 23 observed food webs (Neutel & Thorne, 2014 , table S3): two Antarctic food webs with interaction strengths quantified from independent flux observations (Neutel & Thorne, 2014 , 2016b ) and 21 soil food webs with interaction strengths quantified from inferred fluxes based on time-averaged biomass observations (de Ruiter et al., 1995 de Ruiter, Neutel and Moore 2016 Neutel et al., 2007 , 2016 ). Our analysis was on the 23 trophic networks (+/− structure, consumer–resource, or predator–prey interactions) in these food webs, obtained by removing the detritus row and column from the matrices.

(1) (2)

### 2.3 Parameterization of Jacobian community matrices

The interaction strengths between the populations are the elements of the Jacobian community matrix, a linearization of the system around the nontrivial equilibrium (where each species has a positive population density). Thus, they are the partial derivatives of the population growth equations (dimension per time) evaluated at equilibrium (May, 1973 ). Using equation 2, the elements of the community matrix Α are effects of predator j on prey i, effects of prey h on consumer i, and intraspecific effects, (because at equilibrium, . Thus, the negative effects of predators on prey are the feeding rates of a predator on its prey divided by predator biomass, and the positive effects of prey on their predators are predator growth rates divided by prey biomass (Pimm, 1982 de Ruiter et al., 1995 Yodzis, 1989 ).

The interspecific interaction strengths of the Antarctic food webs were quantified from direct flux observations (Neutel & Thorne, 2014 , 2016b ), and those of the 21 soil food webs were quantified from inferred fluxes based on time-averaged biomass observations (de Ruiter et al., 1995 de Ruiter et al. 2016 Neutel et al., 2007 , 2016 ). The intraspecific interaction strengths (diagonal elements) were not obtained from directly observed fluxes, but the observations provided upper bounds. The rationale is as follows: For a given food web in equilibrium, total loss for each population equals total gain. Growth rates and predatory loss rates of each population were known hence, nonpredatory loss rates were also known, because the systems were in equilibrium. This total nonpredatory loss rate consists of intrinsic death and intraspecific competition: . The amount of intraspecific competition is hence contained within the energetic boundaries of the system . The upper bounds were used to quantify the diagonal elements of (un-normalized) community matrix Α: .

### 2.4 Scaled interaction strengths

Following Neutel & Thorne, 2014 interaction strengths were scaled by dividing each row in community matrix Α by the absolute value of its respective diagonal element, , which resulted in time-independent and dimensionless matrices Γ. This scaling procedure was introduced by Neutel and Thorne ( 2014 ) to translate the diagonal structure of the matrix into the off-diagonal structure, and obtain an eigenvalue which, for their observed food webs, was equivalent to a critical value of intraspecific competition for stability [specifically, it represents the proportion of total nonpredatory loss needed for stability see Neutel et al. ( 2002 )].

### 2.5 Determination of stability

The diagonal elements of the Jacobian and scaled matrices were then set at zero, obtaining matrices Α0 and Γ0. System stability was determined as the largest real part of the eigenvalues of these matrices, λd. By definition, with all diagonal elements set at zero, then λd ≥ 0 that is, the systems need some level of self-damping in order to be stable, and −λd is the amount of self-damping needed for stability. In the case of the normalized matrix Γ0, λd has a biological meaning and indicates a tipping point. It represents a critical level of intraspecific competition of the populations as a proportion of the maximum possible intraspecific competition (the upper bound of the diagonal elements) (Neutel & Thorne, 2014 ). In the case of the Jacobian matrix Α0, λd is related to the timescales of the systems and is not easily interpretable biologically given the different intraspecific interaction strengths of the populations. It is the opposite of a system's resilience. The inverse of λd of Α0 is the time with which the system moves away from the equilibrium, after a very small disturbance.

### 2.6 Synthetic parameterization of community matrices

For the synthetic parameterization of the community matrices, the empirical values of the nonzero off-diagonal elements were first replaced by values randomly drawn from uniform distributions (−10, 0) for effects of predators on prey and (0, 0.1) for effects of prey on predators (following Pimm & Lawton, 1978 see also Neutel & Thorne, 2014 ). This procedure was then repeated without the asymmetry in size ranges, drawing values from uniform distributions (−1, 0) and (0, 1) (following May, 1972 see also Neutel & Thorne, 2014 ).

### 2.7 Pairwise disturbance experiment

To show the effect of the patterning of interaction strengths on stability, we performed a disturbance experiment, following Yodzis ( 1981 ). For each empirical community matrix Α0, the pairs of nonzero off-diagonal elements were randomly permuted. This preserved the sign structure and the pairwise structure of the interaction strengths of the original matrix.

### 2.8 Feedback metric

The metric proposed by Neutel and Thorne ( 2014 ) expresses a ratio between three-link and two-link feedback: , where a2 and a3 are coefficients of the characteristic polynomial of the community matrix. In a matrix with zero-diagonal elements, the second coefficient, a2, represents the sum of all the two-link feedback loops (F2), resulting from the pairs of predator–prey interactions, which are by definition negative: In community matrix Α0, (the same holds for normalized matrix Γ0). The third coefficient, a3, in zero-diagonal matrices, represents the sum of all the three-link loops (F3), coming from the smallest omnivorous structures, each generating a positive and a counteracting negative feedback loop: , where i is the bottom prey, j is the intermediate predator, and k is the omnivore. The sum of the positive and counteracting negative loop in each three-link omnivorous structure is by definition positive, given the functional assumptions (Neutel & Thorne, 2014 ).

### 2.9 Pairwise metric

The metric proposed by Tang et al. ( 2014 ) quantifies the overall correlation between effects of predators on prey and vice versa: , where S is the number of “species,” E is the mean of the off-diagonal elements of the community matrix, V is their variance, and ρ is the overall pairwise correlation between the elements of the community matrix (αij, αji)ij (Tang et al., 2014 ).

## Results

### Independently measured fluxes vs. inferred fluxes

We quantified the two entire above-and below-ground food webs of the Antarctic dry and wet tundra ecosystems (Fig. 2a) including recycling through dead organic matter and found that the overall material flow in both systems approximated a balanced carbon cycle (Fig. 2b).

The material fluxes together with the observed biomass densities of the populations provided us directly with the element values of the Jacobian community matrices, including the structure of the diagonals. We translated the observed structure of the diagonals into the off-diagonal elements and evaluated the stability properties of the ecosystems without assuming any self-damping through intraspecific losses (competition) by the populations. The food web of the Antarctic wet tundra was substantially more vulnerable than that of the Antarctic dry tundra, λd = 0.26 and λd = 0.04 respectively. We compared the community matrices based on the independent material-flow observations with matrices of the Antarctic systems where interaction strengths were calculated from inferred material fluxes deriving them from the steady-state assumption (Fig. S2). System stability of the Antarctic webs was remarkably similar to that of their inferred-flux counterparts (λd = 0.24 for the inferred wet and λd = 0.08 for the inferred dry tundra). This indicated that the indirect method of inferring fluxes from steady-state biomass observations reflected the direct observations sufficiently to capture the essential stability properties. We tested the robustness of the results by sampling feeding rates and biomass from the observed variability in these values (Fig. S4).

### Effect of recycling on ecosystem stability

To investigate what determined the vulnerability of the Antarctic ecosystems, we first tested the effect of recycling of organic matter on the stability of the systems. We found that removing the recycling effects (interactions between populations and detritus) neither increased nor decreased the stability of the systems: they had no effect on system stability. We explored the generality of this result by drawing in 39 soil food webs with parameter values derived from steady-state biomass observations. Two categories of systems emerged: for 21 food webs recycling had no or a negligible effect on system stability, for 18 food webs recycling had a clear effect, either stabilising or destabilising the system (Fig. 3). It is worth noting that the latter were all communities from relatively young soils, low in organic matter, or with relatively simple network structure (Table S3). We then continued our analysis with the Antarctic matrices, their inferred-flux counterparts and the 21 soil food-web matrices that shared the property that only predator-prey feedbacks (feedbacks without interactions with detritus) affected system vulnerability.

Examination of all the feedback loops in the systems confirmed earlier findings that longer loops tended to be relatively weak (Neutel et al. 2002 , 2007 Emmerson et al. 2005 ). The coefficients in the characteristic polynomials of the matrices also tended to be relatively small (in absolute value) on higher levels. The signs of the coefficients showed that our systems, with no self-damping by the populations, were systems of predominantly positive feedback. If we increased the level of self-damping up to the tipping point at which there was just enough self-damping for stability, this was also the point where positive feedback ceased to be predominant and negative feedback prevailed. Thus, we concluded that the stability of our systems depended on the balance of positive and negative feedbacks, not on the absence of time lags caused by excessive negative feedback (see SI for details). We found that the vulnerability of the systems could roughly be predicted by the net positive feedback resulting from the 3-link omnivorous interactions and the negative feedback from the 2-link predator-prey interactions according to , where a2 and a3 are coefficients of the characteristic polynomial of the community matrix (Box 1 Fig. S3).

This can be understood as follows. For a system to be stable, the sum total of feedback of every length has to be negative. The net positive feedback of the 3-link loops can only be compensated for by the negative feedback resulting from combinations of pair-wise predator-prey feedback and self damping effects. The stronger the 2-link negative feedback, the less self-damping needed. One could think of the positive 3-link feedback as the elementary destabilising force in trophic networks, balanced by the stabilising feedback from 2-link loops and self-damping.

### Symmetry, asymmetry, complexity: Empirical vs. synthetic parameterisation

The absence of an effect of detritus feedbacks on system stability allowed us to make a direct comparison between the empirically-based parameterisations and synthetic (symmetric and asymmetric) parameterisations underlying the canonical ideas on the destabilising effects of diversity or link density (Gardner & Asby 1970 May 1972 ) and omnivory (Pimm & Lawton 1978 ). The ratio roughly predicted system stability irrespective of the type of parameterisation (Fig. 4a).

However, did not capture the difference between the symmetric webs. The feedback spectrum of the symmetric webs explained why this was the case (Fig. 4e). In the symmetric webs, the lack of structure in the interaction strengths created a certain homogeneity at system level, a symmetry in feedback strength between positive and negative feedback loops of the same length. This balance between positive and negative 3-link loops meant that it was just the 2-link loops that governed stability. But instead of contributing to stability, as in the asymmetric and empirical parameterisations, the 2-link predator-prey loops showed a negative relationship with system stability. The number of 2-link loops is equal to the number of interactions in the system and we found indeed that the number of links per species (degree of the network), represented by the classic complexity metric (May 1972 ), where n is the number of species and C is the connectance, was strongly correlated with the vulnerability of the symmetric webs (Fig. 4b).

In the asymmetric webs, the negative effects of predators on prey were much stronger than the positive effects of prey on predators. Thus, the strength of a feedback loop depended on the proportion of negative and positive effects in that loop, so that, on average, positive 3-link loops (containing two negative effects and one positive effect) were much stronger than negative 3-link loops (Fig. 4e). This made the asymmetric webs very vulnerable. In these webs, was effectively determined by the number of 3-link loops relative to 2-link loops, since values for positive and negative interactions were all randomly sampled from the same two intervals. The degree of omnivory (Pimm & Lawton 1978 ), here calculated as , where omn3 is the number of 3-link omnivorous structures and L is the total number of interactions, was strongly correlated with the vulnerability of the asymmetric webs (Fig. 4c).

### Affiliations

Department of Mathematics, University of Illinois, Urbana, IL, 61801, USA

Department of Plant Biology, University of Illinois, Urbana, IL, 61801, USA

Carl R. Woese Institute for Genomic Biology, University of Illinois, Urbana, IL, 61801, USA

You can also search for this author in PubMed Google Scholar

You can also search for this author in PubMed Google Scholar

### Contributions

S.B. and J.O.D. contributed equally to this work.

## Genome-scale Molecular Analysis𠅊 PANOMICS Platform Combined With Data-Driven Modeling

In the following, I will introduce techniques for genome-scale molecular analysis as well as mathematical and statistical hybrid modeling as fundamental requirements to fulfill this proposed cycle of understanding the dynamics of an organism. The comprehensive analysis of organisms, microbes, plants, fungi, animals, and human begins nowadays with genome sequencing using “next-generation sequencing” technologies [12]. On the basis of a full genome sequence and following the information flow from genome to transcriptome to proteome to metabolome, subsequent genome-scale molecular techniques such as transcriptomics (RNAseq), proteomics, and metabolomics𠅊 PANOMICS strategy𠅊re employed to reveal a complete picture of molecular dynamics of the respective organism (Figure 1). In parallel, comprehensive genome-scale metabolic reconstructions are derived from the genome sequence by comparing orthologous genes and their functions by database search, e.g., using UNIPROT with annotated protein sequences and functions of up to 160 million predicted proteins [15, 23]. A comprehensive discussion of all these initial steps up to the PANOMICS platform to investigate an organism is described in a recent publication entitled the “unpredictability of metabolism from genome sequences” [12]. This obviously contrasting title comes from the fact that a genome-scale metabolic reconstruction does not reflect any dynamics or plasticity of metabolism. In other words, the phenotypic plasticity of an organism that is the result of the interaction with its environment cannot be derived from the genome sequence alone, at least not at this stage of biological and biochemical knowledge. A hope is, of course, to apply system-theoretical ideas and, thus, systems biology to this problem of functional genome annotation to finally give a causal𠄿unctional interpretation of a newly sequenced genome up to its molecular and phenotypical dynamics in relation to its environment. It is further obvious that such a prediction of function and phenotype cannot be derived by looking at separated parts of the organism but rather by looking to the relation or network dynamics of all biochemical and physiological components. Thus, the organism has to be understood as a whole. This is not a contradiction to sophisticated biochemical analysis of isolated parts of the system, and, thus, often a misinterpretation of system theory and consequently systems biology. The more knowledge about single parts of the system there is, the better the description and understanding of the whole will be.

Figure 1. A workflow which combines multiomics analysis [36] and causal modeling [22] of molecular co-regulation (also called correlation or association) networks. Because the linkage of genome-scale metabolic reconstruction (RECON) and the data covariance matrix (COV) is central in this approach we call it COVRECON. The calculation of the differential Jacobian is demonstrated for the first time by Sun and Weckwerth [22].

Another very prominent example of functional annotation of genomes derives from genome-wide association studies (GWAS). Here, allelic polymorphisms in the genome of one species and its corresponding ecotypes or phenotypes are screened systematically to reveal their correlation with phenotypic and adaptive traits. In most cases, single-nucleotide polymorphisms (SNPs) are molecular neutral mutations however, there are phenomena such as linkage disequilibrium and accumulation of SNPs in a specific genomic region pointing to potential adaptation processes in the genome [24�]. Initially, these potentially functional SNPs are revealed by correlation analysis with phenotypic traits, and in a tedious and laborious way, these spurious correlations need to be tested for their causal relevance. A very elegant study has demonstrated the full workflow 2008 in maize starting with GWAS and finally identifying the causal relationship in some maize seeds [29]. Other studies have combined GWAS and metabolomics to identify causal genes involved in the human metabotype [30]. We have recently applied in situ eco-metabolomics in combination with SNP enrichment and metabolic modeling to reveal potential biochemical adapation processes of Arabidopsis thaliana to the natural habitate and micro-environment [31]

The final and conclusive step to link genome-scale molecular analysis with genome sequence information is model building [12]. These models can be pure structural, kinetic, or mixed forms, which allow one to analyze the stability of the system by systematically perturbing kinetic parameters [32, 33].

However, all these approaches rely on “speculated” forward conjectures about network structure and kinetic parameters. The problem is that all these data are not available, especially the detailed kinetics of protein interaction, enzymatic reactions, and biochemical regulation. We have recently introduced a novel idea by linking correlation networks of molecular components in an organism with the underlying biochemical regulation [34]. The dynamic molecular interaction or correlation networks are derived by applying the comprehensive PANOMICS platform introduced above and using statistics to reveal correlations between the components [34�] (Figure 2). Here, the terminus correlation is equivalent to covariance, co-regulation or association often used in literature. Any kind of pairwise correlation analysis or multivariate statistical analysis is using the same basic assumption that there is a relation between the components that can be resolved by statistics. In a second step, we use genome data for a metabolic reconstruction as introduced above [12, 17]. The linkage of the molecular correlation networks and the genome-scale metabolic reconstruction somehow fills the static information of the genome with life. Eventually, we have demonstrated in 2012 that one can use the molecular covariance networks to calculate the biochemical perturbation of the system using the inverse stochastic Lyapunov matrix equation [22] (Figure 1). In several follow-up studies, we have demonstrated that the �lculated” predictions—in contrast to “speculated” predictions—reveal biochemical key points of regulation [16, 39�]. Also, other groups have used the proposed concept to predict biochemical perturbation [42, 43]. The aim of the following paragraph is to demonstrate that the concept of the stochastic Lyapunov matrix equation is systematically derived from very basic system-theoretical principles, thereby presenting a fundamental system-theoretical equation.

Figure 2. Derivation of metabolite correlation networks from large-scale metabolomic analysis. These correlation networks show differential connectivity in dependence of the analyzed genotype. From this fundamental multivariate structure several research fields are envisioned: Firstly, topology analysis of the correlation network resulted in a power-law-like pattern of the probability of metabolite connectivity [34, 35]. The power-law pattern is observed in real world networks with only a few hub nodes showing high connectivity in contrast to most others. This pattern is also highly conserved in biochemical networks, thus reflecting biochemical network structure and regulation. Secondly, based on the correlation pattern we developed a stochastic model of metabolism which led to the stochastic Lyapunov matrix equation (for more details see text).

## Old good Gaussian fit¶

To illustrate how to use ABC within PyMC3 we are going to start with a very simple example estimating the mean and standard deviation of Gaussian data.

Clearly under normal circumstances using a Gaussian likelihood will do the job very well. But that would defeat the purpose of this example, the notebook would end here and everything would be very boring. So, instead of that we are going to define a simulator. A very straightforward simulator for normal data is a pseudo random number generator, in real life our simulator will be most likely something fancier.

Defining an ABC model in PyMC3 is in general, very similar to defining other PyMC3 models. The two important differences are: we need to define a Simulator distribution and we need to use sample_smc with kernel="ABC" . The Simulator works as a generic interface to pass the synthetic data generating function (normal_sim in this example), its parameters, the observed data and optionally a distance function and a summary statistics. In the following code we are using the default distance, gaussian_kernel , and the sort summary_statistic. As the name suggests sort sorts the data before computing the distance.

Finally, SMC-ABC offers the option to store the simulated data. This can he handy as simulators can be expensive to evaluate and we may want to use the simulated data for example for posterior predictive checks.