Analysis of trophic networks: an optimisation approach
Jean-Guy Caputo, Valerie Girardin, Arnaud Knippel, Hieu Nguyen, Nathalie Niquil, Quentin Nogues
aa r X i v : . [ n li n . AO ] F e b Analysis of trophic networks: an optimisationapproach
Jean-Guy Caputo ∗ , Valerie Girardin † , Arnaud Knippel ‡ , HieuNguyen § , Nathalie Niquil, ¶ and Quentin Nogu`es, ‖ Laboratoire de Math´ematiques, INSA de Rouen Normandie,76801 Saint-Etienne du Rouvray, France. UNICAEN, CNRS, Laboratoire de Math´ematiques NicolasOresme, 14000 Caen, France UNICAEN, Laboratoire Biologie des ORganismes et Ecosyst`emesAquatiques, FRE 2030 BOREA (MNHN, UPMC, UCBN, CNRS,IRD-207) CS 14032, 14000 Caen, FranceFebruary 4, 2021
Abstract
We introduce a methodology to study the possible matter flows of anecosystem defined by observational biomass data and realistic biologicalconstraints. The flows belong to a polyhedron in a multi dimensionalspace making statistical exploration difficult in practice; instead, we pro-pose to solve a convex optimization problem. Five criteria correspondingto ecological network indices have been selected to be used as convex goalfunctions. Numerical results show that the method is fast and can beused for large systems. Minimum flow solutions are analyzed using flowdecomposition in paths and circuits. Their consistency is also tested byintroducing a system of differential equations for the biomasses and exam-ining the stability of the biomass fixed point. The method is illustratedand explained throughout the text on an ecosystem toy model. It is alsoapplied to realistic food models. keywords : Convex optimization Ecosystem Trophic networkAMS classification: 92C42, 49N30 ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] ¶ [email protected] ‖ [email protected] cknowledgements This work is part of the ECUME project, co-financed bythe European Union with the European regional development fund (ERDF) andby the Normandie Regional Council.
Functional ecology is based on seminal works from the XXth century centeredon the object ecosystem. From the first definition of an ecosystem by [40] to theconstruction of its main concepts by [26], [34] or [28], among others, ecosystemshave been described as entities gathering living organisms and their habitat,and described as dynamic entities, based on exchanges of organic matter. Fromthose works were derived a system analysis of these exchanges based on emergentproperties; see [35], [36], [41], [12].The description of ecosystems is often based on networks of interactions, ofdifferent types. For terrestrial ecosystems, recent developments concern differenttypes of interactions, sometimes gathered into a common model called multiplex[10]. In marine ecology, the most studied interactions are trophic, i.e. theinteractions between predators and preys; they form a network called a foodweb . Food webs in marine ecosystem are highly complex, compared to theterrestrial ones [1] and have been described by numerous models. These modelshave been widely used to describe the impact of human activities on marineecosystems [19]. They are also important tools for the sustainable managementof marine and coastal environments [24].The trophic modeling of food webs has been mainly based on weighted net-works; see for example the Ecopath-Ecosim-Ecospace models [6]. There, eachlink corresponds to a transfer of organic matter between two trophic compart-ments, collecting individuals of similar feeding behaviors and metabolisms, andwith the same predators. Some fluxes can be estimated using laboratory exper-iments that are often associated to field studies, however many of them remainunknown. To take into account these unknown flux values within food webs, aclass of models was developed called Linear Inverse Modeling (LIM) [32]. LIMassumes a steady state for the biomass of all compartments – a mass balancedsystem. This yields a set of linear equations (equalities) describing the steadystate or mass balance. Then, constraints are added from field measurements ofmass transfers like local estimations of primary production, respiration or dietcontents. Additional constraints come from experiments or the study of otherecosystems. All these constraints constitute a set of linear equalities and in-equalities defining a bounded multidimensional polyhedron, called a polytope,within which lie all realistic solutions to the problem. Such solutions are termedflows in graph theory.In the literature on ecosystems, the polytope is explored using a randomwalk method, called Monte Carlo Markov Chain – see [21],[22], [44] and [45] –or Monte Carlo Linear Inverse Modeling (MC-LIM) – see [21] and [44]. Linear2 i ff f f Figure 1: Schematic drawing of the polytope with the optima for different goalfunctions.Inverse Monte Carlo Markov Chain (LIM-MCMC) models are mass balanced.This stochastic approach is an indirect way to consider the variability of theliving; see [32] [45]. As such, they provide a wide range of possible results,and not a single value like other approaches. However, for large systems, thisexploration of the polytope can be very long, with a very large number ofsimulation runs.Gathering indices from various domains, from information theory to input-output analysis in econometry, several Ecological Network Analysis indices havebeen introduced for describing the organization of the flows and the functioningof the ecosystem [25], or as criterions of ecosystem maturity [42]; One suchsimple index is the sum of the flow components squared [46]; see [18], [22], [38],and [43] for many others. These indices assume that the flows are given, butLIM results can be used to compute approximated values.In this article, we propose to use these indices as goal functions and combinethem with the constraints to set up an optimization problem. This procedurehas a low computing time compared to the LIM MC-MC method and directlyyields a unique flow solution within the polytope, if the goal function is con-vex. Note that [38] combined the LIM-MCMC exploration of the polytope tocomparison of different indices to select a unique flow vector.Fig. 1 presents a schematic picture of the polytope in flow space f togetherwith three minima corresponding to three convex goal functions. To comparethe optima obtained through optimization of goal functions, say f , f and f ,a first step is to examine the main flows from an ecological point of view and seeif they appear reasonable. In a more quantitative approach, one can decomposethese flows into paths and circuits and again check these using ecology know-how. We also suggest, using simple rules, to introduce a dynamical systemsatisfied by the biomasses and whose coefficients depend on the flow solution f , f or f . This dynamical system has a fixed point – the given biomasses,whose stability can be determined. If the fixed point is stable, then the model3s consistent, say for example f in Fig. 1. Then the optimum f yields anacceptable solution to describe the ecosystem.Using a six species toy model inspired by a realistic ecosystem, we proceedto illustrate this methodology. We do not pretend that the model is realisticbut we focus on the analysis and present it in as much detail as possible. Thisdetailed presentation is easy to follow on the six species system and naturallyextends to an ecosystem of any size. We write explicitly the constraints definingthe polytope in Section 2. In Section 3, we discuss flow decomposition intopaths and circuits, a general result from the theory of polyhedra. The fulloptimization problem is presented in Section 4. First we present the five convexgoal functions, three of which are independent of the constraints and two dependon the constraints. The results of the optimization problem are analyzed usingthe flow decomposition of Section 3. From these flow solutions, we write theformal dynamical system for the biomasses and examine the stability of the fixedpoint in Section 5. We show that the fixed point is always marginally stable,in the absence of detritus, in Section 6. We show that the detritus controlsthe stability and give a sufficient condition for the fixed point to be stable.Conclusions and application to real data are presented in Section 7. A realistic model of a marine ecosystem with nineteen species was introducedand analyzed by the authors in [33]. To focus on the method of analysis, wesimplified this model and reduced it to an ecosystem of six species. This method-ology can be extended to ecosystems or arbitrary size; in Section 6, we give someresults for the realistic ecosystem studied in [33].The graph of this simplified 6-species ecosystem is presented in Fig. 2. Theordered types of living organisms with biomasses – circles in Fig. 2 – are:Phytoplankton ≡ PHY1 , Zooplankton ≡ ZOO3 , Bivalve ≡ BIV4 , Fish benthic feeders ≡ FBF5 , Bacteria ≡ BAC6 . The other vertices – rectangular boxes in Fig. 2 – are:Detritus ≡ DET2 , Photosynthesis of phytoplankton ≡ FIX7 , Respiration ≡ RES9 , Fishing, Trawl, Dredge,... ≡ LOS10 , Import to system (ability of a species to move geographicallyin order to feed) ≡ IMP8 . The arrows between the compartments (edges between vertices) represent mat-ter flows that are inferred by ecologists. Following graph theory, terminologysuch oriented edges will be denoted arcs.We denote by S the set of all vertices including the detritus DET2, S = { P HY , DET , ZOO , BIV , F BF , BAC } , ET2ZOO3 BIV4 BAC6FBF5IMP8LOS10 RES9PHY1 FIX7
Figure 2: The 6-species model ecosystem described in Section 2.5ith S ′ = S − { DET } the set of all species vertices, and E the set of all arcs– denoted by ij when going from vertice i to vertice j .The central mathematical object of this article is the flow. Definition 2.1
Let the successors and predecessors of a vertice i be N + ( i ) = { vertices j, ij ∈ E } and N − ( i ) = { vertices j, ji ∈ E } . A flow is a vector, of dimension the number of arcs with non negative com-ponents, satisfying Kirchoff law on all vertices of S , X j ∈ N + ( i ) f i,j − X j ∈ N − ( i ) f j,i = 0 , i ∈ S, (1)In many situations, all biomasses B i of the species i can be measured withaccuracy. On the other hand, the flow components f i,j are much more difficultto evaluate. Therefore, we will adopt here the standard viewpoint that thebiomasses are given and the flows between nodes are unknown. The flow components f i,j between the compartments satisfy biological con-straints. For example, a fish cannot eat more that a certain percentage ofits biomass. When defining the constraints, we gather all available information,if possible from studies of the local ecosystem, if not, from ecosystems similar toours. In absence of information, the constraints are derived from experimentsor from empirical equations. Definition 2.2
For all species i ∈ S , the production P i is P i = P j ∈ S ′ f j,i − f i, res − f i, det , where f i, res is the respiration flow and f i, det is the excretion flow – assuming itgoes to the detritus. Among non species compartments (vertices), the detritus plays a singularrole as the only one for which flows go in and out of. It is then natural to assumefor it a Kirchoff law where in-going equal out-going flows. For the species verticesin S ′ , the following constraints are imposed by biological observations throughnonnegative coefficients c . These coefficients come from field measurements ofmass transfers like local estimations of primary production, respiration or dietcontents. They can also be estimated from experiments or the study of otherecosystems. Positivity of the flow components f i,j ≥ , for all ij ∈ E. irchoff law at the vertices Equations (1) express the conservation of massat each species vertex.
Primary production constraint
The production P of the entry in theecosystem Phy1 is bounded, c − pro ≤ P ≤ c +pro . This constraint comes from a local study estimating the carbon incorpo-rated, with enrichments in 13Cu a stable isotope of carbon, compared withstudies based on an estimation of the activity of photosystems using pulseamplitude modulation [30].
Respiratory constraints
The respiration flow f i, res of each species i is bounded, c − res ,i X j ∈ N − ( i ) f j,i ≤ f i, res ≤ c +res ,i X j ∈ N − ( i ) f j,i , for all i ∈ S ′ . Excretion constraints
The excretion of all species i but the phytoplanktonis bounded, c − det ,i X j ∈ N − ( i ) f j,i ≤ f i, det ≤ c +det ,i X j ∈ N − ( i ) f j,i , for all i ∈ S ′ , i = PHY1 . The phytoplankton excretion is bounded too, c − det , phy P phy ≤ f phy , det ≤ c +det , phy P phy . Food conversion efficiencies constraints
The production of a species i isconstrained by bounds depending on the entering flow of species ic − eff ,i X j ∈ N − ( i ) f j,i ≤ P i ≤ c +eff ,i X j ∈ N − ( i ) f j,i , for all i ∈ S ′ . Production related to biomass constraints
The production of a species i is constrained by bounds depending on its biomass B i , c − bio ,i B i ≤ P i ≤ c +bio ,i B i , for all i ∈ S ′ , Diet constraints
The entry flow f j,i for species i is constrained by boundsdepending on the sum of the entry flow of this species, c − diet ,i X j ∈ N − ( i ) f j,i ≤ f j,i ≤ c +diet ,i X j ∈ N − ( i ) f j,i , for all i ∈ S ′ . Diet information is derived from stomach content analyses; see e.g. [9].7ifferent empirical equations are gathered to define relationships betweenproduction and biomass or consumption and biomass, or respiration and con-sumption; see [3], [11]. These equations use the shape of the caudal fin, theindividual weight, temperature, growth, etc. The individual mass and totalbiomass per km values are estimated from local field studies [8] or field studiesfrom a similar ecosystem if not available; see, e.g., for zooplankton [37].All above constraints can be summarized in the following set of equations,where f i is the total flow entering a species i . f i,j ≥ , ij ∈ E, (2) f i = P j ∈ N − ( i ) f j,i , i ∈ S, (3) P j ∈ N + ( i ) f i,j − P j ∈ N − ( i ) f j,i = 0 , i ∈ S, (4) f − f , res − f , det − c +pro ≤ , (5) − f + f , res + f , det + c − pro ≤ , (6) − c +res ,i f i + f i, res ≤ , i ∈ S ′ , (7) c − res ,i f i − f i, res ≤ , i ∈ S ′ , (8) − c +det ,i f i + f i, det ≤ , i ∈ S ′ , i = PHY1 (9) c − det ,i f i − f i, det ≤ , i ∈ S ′ , i = PHY1 (10) − c +det , phy P phy + f phy , det ≤ , (11) − f i + f i, res + f i, det + c − eff ,i f i ≤ , i ∈ S ′ , (12) f i − f i, res − f i, det − c +bio ,i B i ≤ , i ∈ S ′ , (13) − f i + f i, res + f i, det + c − bio ,i B i ≤ , i ∈ S ′ , (14) f j,i − c +diet ,i,j f i ≤ , j ∈ S, i ∈ S ′ , (15) − f j,i + c − diet ,i,j f i ≤ , j ∈ S, i ∈ S ′ . (16)For the 6-species ecosystem defined in Fig. 2, these constraints lead tothe bounds on the flows shown in Table 1. More precisely, for example, thebounds on f , are obtained by minimizing or maximizing f , together withthe constraints (2) to (16). Since the constraints (2) to (16) are all linear, the defined domain is polyhedral.A polyhedron in R n is an intersection of a finite number of half-spaces, in otherwords P = { x ∈ R n | Ax ≤ b } , where A is an n ′ × n matrix with n ′ > n and b ∈ R n . A bounded polyhedron is called a polytope.The decomposition P = Q + C for all polyhedrons P , with Q a polytope and C a cone, is classical; see Nemhauser and Wolsey [31]. Since Q and C are convex8omponents bounds f , [10.300 ; 51.600] f , [0.000 ; 136.181] f , [35.818 ; 172.000] f , [0.000 ; 96.498] f , [5.963 ; 95.828] f , [0.000 ; 129.668] f , [0.000 ; 27.892] f , [9.768 ; 92.615] f , [0.000 ; 271.796] Components bounds f , [0.000 ; 82.036] f , [0.000 ; 16.407] f , [0.000 ; 39.002] f , [0.000 ; 35.445] f , [0.000 ; 49.222] f , [0.000 ; 78.113] f , [11.722 ; 111.138] f , [0.000 ; 35.445] f , [18.235 ; 197.447] f , [0.000 ; 78.260]Components bounds f , [0.638 ; 35.445] f , [0.000 ; 6.380] f , [3.190 ; 31.901] f , [0.000 ; 6.380] Components bounds f , [0.000 ; 112.581] f , [0.000 ; 180.966] f , [0.000 ; 192.996] f , [119.263 ; 319.428] f , [0.000 ; 35.445]Table 1: Bounds on the flow components for the 6-species ecosystem of Fig. 2.sets, any point in P can be expressed as the sum of a convex combination of theextreme points of Q (the vertices) and a combination of the extreme rays of C with non negative coefficients : x = X i ∈ I α i q i + X j ∈ J β j r j , (17)where I is the index set of the vertices of Q , J the index set of the extreme raysof C and the coefficients α i satisfy P i ∈ I α i = 1 and α i , β j are non negative.For webs of flows, the vertices of Q are indicator vectors of elementary pathsand the extreme rays are – up to a constant coefficient – indicator vectors ofelementary circuits of the web (paths that begin and ends at the same vertex).Therefore any numerical solution of such an ecological system can be interpretedin terms of paths and circuits. A linked important notion is the flow value. Definition 3.1
For a given path of circuit, the flow value is the smallest fluxfor all arcs of the path or circuit.
Usually, the decomposition (17) is not unique. However, the largest flow valuesof the network will make some paths or circuits necessary in any decomposi-tion, and hence important for interpreting the obtained ecosystem solution. Analgorithm to find all necessary paths or circuits is the following:Repeat • Find the path or circuit with largest value α ;9 Take out α from the flow components of the arcs of this path or circuit;until stopping criterion is met.Such an analysis assumes that the numerical values of the flows are knownand satisfy the biological constraints. Unfortunately, such exact numerical val-ues are generally unavailable, and their a priori approximated values do not sat-isfy the constraints especially the conservation ones. The next section presentsan optimization approach addressing this issue. A flow solution is computedfrom the knowledge of only the biomasses and some approximated values ofthe flows, or intervals of approximated values. Then, examples of paths andcircuits extracted using the above method of flow decomposition are given forthe 6-species network. We consider five goal functions, the most classical least squares, Ecological Net-work Analysis indices, and functions adapted from information theory. All areconvex, yielding a unique optimum corresponding to a unique functioning stateof the system. An important characteristics is whether the goal function in-cludes information from the constraints or not. This separates the functions intwo classes.In one class, no information is used from the constraints. The first functioncorresponds to the quadratic energy, the classical least squares method, and thesolution will be the minimum of F ( f ) = X ij ∈ E f i,j . (18)The second function is minus Shannon entropy from information theory [5],introduced for ecological systems in 1955 in [28], F ( f ) = X ij ∈ E p i,j ln( p i,j ) = X ij ∈ E f i,j f .. ln (cid:18) f i,j f .. (cid:19) , (19)where the sum of all flows is f .. = P ij ∈ E f i,j , and the proportion of flows fromvertex i to vertex j is p i,j = f i,j f .. . (20)Classically, entropy is a concave function and is to be maximized. Here, for prac-tical purposes, we only deal with convex goal functions and therefore changedthe sign. With this sign, F is convex – see Appendix 2.10inally, the third function is minus the system redundancy (overhead) in-troduced in [42] as F ( f ) = X ij ∈ E p i,j ln (cid:18) p i,j p i. p i,j p .j (cid:19) = X ij ∈ E f i,j f .. ln f i,j f i. f .j ! , (21)where the marginal proportions are p i. = X j ∈ S : ij ∈ E p i,j and p .j = X i ∈ S : ij ∈ E p i,j ;note that the sum is on all i or j such that ij is an arc of E . As for the entropy,with this choice of sign, F is convex – see Appendix 2.In information theory terms, − F is the Shannon entropy of the system and − F is a symmetrized conditional entropy; see [5] and also [43] for details onsuch entropic indexes. Note that, although it is a well-known Ecological NetworkAnalysis index, we do not consider the ascendency of [15] because it is neitherconvex nor concave – see Appendix 2.In the other class, the goal functions incorporate information from the con-straints. The most classical quantity in this aim in information theory is theKullback-Leibler divergence introduced in [23], that measures a ”distance” be-tween two distributions, F ( f ) = K ( f | f ∗ ) = X ij ∈ E p i,j ln p i,j p ∗ ij ! , (22)where the proportions are given in (20). The divergence is not a mathematicaldistance because it is not symmetric in p and p ∗ . Still, it is nonnegative andnull only if p = p ∗ , and minimizing K determines the projection in terms ofdivergence of the reference f ∗ (or p ∗ ) on the set of solutions f to the constraints;see [5] and [7]. A natural way, that makes sense in ecology, to include theinformation from the constraints is to set all f ∗ ij as the middle of the constraintintervals [ f ij min , f ij max ]. Note that F is the Kullback-Leibler divergence where f ∗ is the uniform distribution on the flows. This uniform distribution is usuallynot a flow; this confirms the importance of using graph theory to describe suchsystems.Finally, a simple generalization of the quadratic function F is F ( f ) = X ij ∈ E ( f i,j − f ∗ i,j ) . (23) Combining the biological constraints (2) to (16) with the goal functions (18) to(23) yields the well posed convex optimization problem in a positive polytope11f R n , min f F, (24)with constraints (2) − (16) , where F can be any one of the five goal functions given in (18) to (23).A well suited method to solve (24) is the Sequential Quadratic Program-ming (SQP) which uses an Augmented Lagrangian Solver; see [2] containedin the R library NlcOptim [4]. We used the R software infrastructure for theoptimization. The algorithm is presented in Appendix 1.The flows corresponding to the optimum for the five different goal functionsand the 6-species ecosystem are presented in Table 2. Using the flow decom-position estimation introduced in the previous section, we computed the mainpaths and the circuits for each flow solution. These are shown respectively inblack and green in Fig. 3.For instance for the graph minimizing F , the main paths/circuits are:FIX7 −→ PHY1 −→ LOS10 , FIX7 −→ PHY1 −→ BIV4 −→ RES9 , (25)BIV4 −→ DET2 −→ BIV4 . Fig. 3 shows the graph with the largest obtained flows. The flows obtainedfor F and F appear as the most expected results as the largest flow componentsare located at the low trophic levels, corresponding to the classical pyramidalview of energy flows proposed by [26] and [17].Note that the common path (25) appears for the five optimization functions,showing the well-known importance of phytoplankton in coastal ecosystems.This path is the main input of energy in the system. Pathways with largestflows also display the role of the bivalves in the system: they are importantconsumers of phytoplankton and producers of detritus. Their role as recyclersshows in solutions for all goal functions. In such ecosystems, bivalves are notthe only recyclers, bacteria also recycle through the bacterial loop. This showsin solutions yielded by F and F . One odd path is the importation path tofish FBF in the solution of F ; its importance displays a state of the food webwhere the connection to other systems is crucial for the high trophic level FBFcompartment. This solution corresponds to an extreme case, where the systemis unable to sustain the FBF compartment and thus relies more on importation.12omponent Solution1 Solution2 Solution3 Solution4 Solution5 f , f , f , f , f , f , f , f , f , f , f , f , f , f , f , f , f , f , f , f , f , f , f , f , f , f , f , f , F to F . 13 ET2ZOO3 BIV4 BAC6FBF5IMP8LOS10 RES9PHY1 FIX7119 3925 1033 29 1217 1815 DET2ZOO3 BIV4 BAC6FBF5IMP8LOS10 RES9PHY1 FIX7121 3646 27 502326 2323 25DET2ZOO3 BIV4 BAC6FBF5IMP8LOS10 RES9PHY1 FIX71313641 27 33 2016182317DET2ZOO3 BIV4 BAC6FBF5IMP8LOS10 RES9PHY1 FIX71815645 4271 2027274434 29 DET2ZOO3 BIV4 BAC6FBF5IMP8LOS10 RES9PHY1 FIX72418451 74140 55585261 52
Figure 3: The largest flows for the 6-species system and the goal functions: top F (left) and F (right), middle F , bottom F (left) and F (right). Red linesare paths and green lines are circuits. 14 Consistency of the solution: stability of theecosystem
Solving the optimization problem (24) for each goal function F i yields a seriesof flows f between the species (including detritus) vertices. From these flows,using a set of biological rules to be given below, we obtain a formal dynamicalsystem where the variables are the biomasses B i . This system involves couplingcoefficients α that will differ for the different goal functions.All these systems have, by construction, the same fixed point correspondingto the given biomasses. Then, by examining the eigenvalue with maximal realpart of the Jacobians of these different dynamical systems at the fixed point, wecan evaluate the structural stability of the ecosystem for this set of biomasses. Each flow solution f of the optimization problem (24) can be used to derive asystem of differential equations satisfied by the biomasses. The following ruleswill be used to derive these systems of differential equations.If living organisms of species i eat organisms of (either living or detritus)species j , the flow f i,j is given by a law of mass action or Lotka-Volterra coupling– see [29], f i,j = α i,j B i B j . (26)Other flows are assumed to be proportional to the biomass of the start vertex,with a coefficient of proportionality depending also on the finish vertex, f i,j = α i,j B i . (27)Thus, the dynamics of the system is given by the flow f = ( f i,j ), that is tosay by α = ( α i,j ).For the 6-species ecosystem of Fig. 2, the formal rules (26) and (27) lead tothe system of differential equations˙ B = f , − f , − f , − f , − f , − f , − f , , (28)˙ B = f , + f , + f , + f , + f , − f , − f , − f , , (29)˙ B = f , + f , − f , − f , − f , − f , − f , , (30)˙ B = f , + f , + f , − f , − f , − f , − f , , (31)˙ B = f , + f , + f , − f , − f , − f , , (32)˙ B = f , + f , − f , − f , − f , . (33)By construction, the fixed point of the dynamical system is obtained at the15ost function F F F F F max( R e ( λ )) 0.093 0.021 0.104 0.149 -0.193Table 3: Maximum values of the real parts of the eigenvalues of the Jacobianassociated to the five goal functions for the 6-species ecosystem.known biomasses of the 6-species ecosystem, B = 3 . , B = 19 , B = 1 . , B = 19 . , B = 3 . , B = 0 . . (34)At the equilibrium point ˙ B i = 0 for all i , and we get, using (28)-(33), α , B − α , B − α , B − α , B B − α , B B − α , B B − α , B = 0 ,α , B + α , B + α , B + α , B + α , B − α , B B − α , B B − α , B B = 0 ,α , B B + α , B B − α , B − α , B − α , B B − α , B B − α , B = 0 ,α , B B + α , B B + α , B B − α , B − α , B − α , B B − α , B = 0 ,α , B + α , B B + α , B B − α , B − α , B − α , B = 0 ,α , B B + α , B B − α , B − α , B − α , B = 0 . A standard way to compute the stability of the fixed point is to evaluate theeigenvalues of the Jacobian matrix of the system. If the real part of one or moreeigenvalues is positive, then the fixed point is unstable.For the 6-species system, the Jacobian is J = (35) − α , B − α , B − α , B α , − α , B − α , B − α , B α , − α , B α , − α , B α , α , − α , B α , B α , B − α , B − α , B α , B α , B α , B − α , B
00 0 α , B α , B α , B α , B Using the optimal solutions obtained in Table 2 for the different goal func-tions, the biomass values (34), and rules (26) and (27), we can compute themaximal real part of the eigenvalues of J , shown in Table 3.Table 3 shows that only F yields a stable biomass fixed point. Further,Fig. 4 shows the spectra in the complex plane of the Jacobians for the five goal16 Re( ) -30-20-10010 I m () -5 -4 -3 -2 -1 0 1 Re( ) -40-30-20-10010203040 I m () Figure 4: The spectra in the complex plane of the Jacobians for the goal func-tions: top F (black × ), F (red o), F (blue *); bottom F (black × ), F (redo). 17unctions. Again, the systems are unstable for F to F . For all four cases, theeigenvalues with maximum real part are close and have a large imaginary partindicating strong oscillations. These goal functions give rise to similar dynamicalbehaviors. On the contrary, F gives rise to a weakly stable fixed point with nooscillations. There are then at least two main regions in the parameter space,giving different behaviors. From line 2 of the Jacobian J in (35), we see that the detritus acts in a differentway than the regular species. To understand this effect, we will detail theequations governing the dynamical systems with and without detritus. First, we assume that the detritus species is absent so that we have a pure Lotka-Volterra coupling, of which the theory is well established. Since the consideredecosystem involve no coupling B i − B i , a Lyapunov function can be found toshow that the equilibrium point is (globally) stable. We will follow [27], see also[14].Consider the Lotka-Volterra model for n species without detritus dB i dt = B i β i + n X j =1 α ji B j , i = 1 , . . . , n, where α ij = − α ji for all i, j for the couplings (26) and β i = P j α ji for thecouplings (27).The nontrivial equilibrium ( ˜ B , . . . , ˜ B n ) is the solution of the system β i + n X j =1 a ji B j = 0 , i = 1 , . . . , n. This equilibrium is feasible, that is ˜ B i > , for all i. Moreover, A is a skew-symmetric matrix so that A + A T = 0, where A T denotes the transpose of A .These two conditions induce that the Lotka-Volterra model is globally stable,see [27].Thus, the detritus is what will determine the stability. Its in and out flowsdepend on the goal function.Let us now consider the system with detritus. The Lotka-Volterra model18ith n species including detritus is dB i dt = B i (cid:16) β i + n X j =1 α ji B j (cid:17) , i = 2 ,dB dt = n X i =1 ,i =2 B i ( α i − α i B ) , where B is the biomass of the detritus, α ij = − α ji for all i, j = 2, α i ≥ α i ≥ i = 2.The Jacobian matrix of this system is J = J + R, (36)where J ij = α ji B i for all i, j , R i = α i , and R = − P i α i B i , and R ij = 0for i = 2. One can see that J = DA, (37)where D is the diagonal matrix of biomasses, with d jj = B j , and A is a skew-symmetric matrix.For example, for the six species system, we have J = − α , B − α , B − α , B − α , B − α , B − α , B α , B α , B − α , B − α , B α , B α , B α , B − α , B
00 0 α , B α , B α , B α , B = DA, where D = B B B B B
00 0 0 0 0 B and A = − α , − α , − α , − α , − α , − α , α , α , − α , − α , α , α , α , − α ,
00 0 α , α , α , α , . Finally, the matrix R is R = α − α B − α B − α B α α α α . .2 Sufficient condition for stability Evaluating how the detritus will change the stability of the system is a difficultproblem. Nevertheless, perturbation theory can help to understand how eigen-values of J get displaced to the ones of J = J + R when the norm of R is muchsmaller than the norm of J , say || R || ≪ || J || ; see [39]. The standard complexinner product on C n will be denoted by ( ., . ) and z is the conjugate of z . Proposition 6.1
With the notation of the previous section, a sufficient condi-tion for the stability of the system with detritus is R e ( v i , Rv i ) ≤ , i ∈ S, (38) where the v i are eigenvectors associated to the eigenvalues λ i of J . proof We will prove that (38) implies that the real parts of the eigenvalues of J are all negative.An approximation λ iJ of an eigenvalue of J is given by λ iJ = λ i + ( w i ,Rv i )( w,v i ) , (39)see [39], where the vector w satisfies J T w = λ i w. (40)In our special context, the eigenvalues of J = DA are zero or pure imaginary(see Lemma 8.1 in Appendix 3), so that without detritus, the fixed point ofthe biomasses is marginally stable. The presence of the detritus will shift thisstability. Let us prove that w = D − v i satisfies (40).Indeed, since D is invertible, we deduce from J v i = DAv i = λ i v i that Av i = λ i D − v i . Since J T = ( DA ) T = − AD , we compute J T w = − ADD − v i = − Av i = − λ i D − v i = − λ i w. (41)Thanks to (37), Lemma 8.1 applies to show that − λ i = λ i . Therefore, J T w = λ i w, and w = D − v i is an eigenvector of J T , satisfying (40).Thus, according to (39), λ iJ = λ i + ( D − v i , Rv i )( D − v i , v i ) . The real parts of λ i are all null, so R e ( λ iJ ) ≤ i ≡ R e (cid:16) ( D − v i ,Rv i )( D − v i ,v i ) (cid:17) ≤ . F F F F F max( R e ( λ )) 0.075 0.127 0.242 0.126 0.233Table 4: Maximum values of the real parts of the eigenvalues of the Jacobianassociated to the five goal functions for the 19-species ecosystem.Since ( D − v i , v i ) = n X j =1 | v ij | B j > i, we only need to consider the sign of the real part of( D − v i , Rv i ) = 1 B ( v i , Rv i ) . Then ∆ i ≤ R e ( v i , Rv i ) ≤
0, and (38) is indeed a sufficientcondition for the stability of the ecosystem.
Our methodology can easily be applied to general realistic ecosystems with alarger number of species and flow components. To show this, let us present theanalysis of the 19-species system defined and studied in [33].Fig. 5 shows the largest flow components similarly to Fig. 3. As for the 6-species system, ZOO3 and BIV4 are the principal phytoplankton-eating animalsin the 19 species system.We built for each solution of the optmisation problem and different goalfunction a dynamical system and estimate the stability of the fixed point. Themaximum of the real parts of the eigenvalues of the Jacobians are reported inTable 4, similar to Table 3. All the flows correspond to an unstable fixed pointshowing that the constraints may need to be refined.We also examined the ecosystems studied in [43]. For all the cases, thebiomass fixed point was found to be stable. For many ecosystems, as for examplethe Crystal River Creek in [43], the detritus component has a large biomass.This induces stability, as was shown above.Above methods may also be applied to human environmental networks, orin economics, for example to economic resource trade flow networks; see [20],[16], and the references therein.To conclude, we introduced a methodology to study the possible flows of anecosystem defined by observational biomass data and realistic biological con-straints. We formalized the constraints and described precisely the polytopecontaining the solutions. We presented a convex optimization problem based21
ET2ZOO4 BIV9 BAC3FBF14 RES21PHY1 FIX20 NEM5 SUP6IPR11IDF7 11920 1138 261214 11 162319 11 DET2ZOO4 BIV9 BAC3FBF14 RES21PHY1 FIX20 NEM5 SUP6IPR11IDF7 18449 36 31 29353324 244125
57 30
DET2ZOO4 BIV9 BAC3FBF14 RES21PHY1 FIX20 NEM5 SUP6IPR11IDF7 20152 36 32 283434 2944
54 36
DET2ZOO4 BIV9 BAC3FBF14 RES21PHY1 FIX20 NEM5 SUP6IPR11IDF7 19943 40 324132 3344
69 36
DET2ZOO4 BIV9 BAC3FBF14 RES21PHY1 FIX20 NEM5 SUP6IPR11IDF7 24758 57 10059 5993
143 6458
Figure 5: The largest flows for the goal functions: top F (left) and F (right), middle F , bottom F (left) and F (right).22n ecological network indices used as goal functions. The method is fast, canbe used for large systems and provides a solution within the polytope accordingto the indices.The minimal flow for each goal function can be analyzed using two differentcomplementary tools. First, the flow is decomposed into principal paths andcircuits. These enable ecologists to discriminate between the different solutions.Second, the consistency of the flow is examined by introducing a dynamicalsystem and studying the stability of the biomasses fixed point. References
Elements of information theory
Appendix 1: SQL local algorithm of Section 4.2
The SQL local algorithm solves a nonlinear minimization problem under con-straints, say: min f F ( f ) ,h i ( f ) = 0 , i = 1 , . . . , m,g j ( f ) ≤ , j = 1 , . . . , n. Let L ( f, α, β ) = F ( f ) − α T h ( f ) + β T g ( f ) denote the Lagrangian functionof this problem where α and β are Lagrange multipliers of dimensions m and n respectively. Let W denote the Hessian matrix of L , defined by W k ≡ W ( f k , α k , β k ) = ▽ ff L ( f k , α k , β k ) , and A ( f ) the Jacobian matrix of the constraints, A ( f ) T = [ ▽ h ( f ) , . . . , ▽ h m ( f ) , ▽ g ( f ) , . . . , ▽ g n ( f )] . Since the cost function F is convex, W k is a positive definite matrix and A ( f ) is a full rank matrix. Then, at an iterate f k , a basic sequential quadraticprogramming algorithm defines an appropriate search direction p k as a solutionto the quadratic programming subproblemmin p p T W k p + ▽ f Tk p, ▽ h i ( f k ) T p + h i ( f k ) = 0 , i = 1 , . . . , m, ▽ g j ( f k ) T p + g j ( f k ) = 0 , j = 1 , . . . , n. This is solved by the following algorithm:26oal function Computing time function evaluations F F F F F F F
7s 55684 F
12s 60905 F
9s 78447 F
1s 1742Table 5: Computing times and function evaluations: top bottom f , α , β );2. For k = 0, 1, 2, 3, . . .Evaluate F f k , ▽ F f k , W k = W ( f k , α k , β k ) , h k , g k , ▽ h k and ▽ g k ;Solve the quadratic subproblem for obtaining p k , α k , β k ; f k +1 ←− f k + p k ; α k +1 −→ u k ; β k +1 −→ v k ; If the condition of convergence is satisfied; STOP with the approximate solution;3.
End(For).
The computing time on an Intel I7 processor and number of function calls aregiven in Table 5. The code will be made available on github.
Appendix 2: Convexity of goal functions
In this appendix, we show that all the goal functions considered in Section 4.1are convex.Both F and F are sums of squares so they are obviously convex. Further,if a function is twice continuously differentiable, then it is convex if and only ifits Hessian is positive semidefinite. We therefore calculate the Hessians H i of F i for i =2, 3, 4. First H = H = p i,j ... p k,l ... ... ... ... ... ... p r,s , F andthe Kullback Leibler divergence F are convex.Let us rewrite the redundancy F as F = X ( i,j ) ∈ E p i,j ln (cid:18) p i,j p i. (cid:19) + X ( i,j ) ∈ E p i,j ln (cid:18) p i,j p .j (cid:19) . The goal is to show that each term T ij ≡ p i,j ln( p i,j p i. ) is convex. The Hessian of T ij is H ij = ( p i. − p i,j ) p i,j p i. − p i. − p i,j p i. ... − p i. − p i,j p i. − p i. − p i,j p i. p i,j p i. ... p i,j p i. ... ... ... ... − p i. − p i,j p i. p i,j p i. ... p i,j p i. Let us use the Sylvester’s criterion to show that T ij is convex; see [13]. Thiscriterion says that a Hermitian matrix is positive semidefinite if and only if allof its leading principal minors are positive.The Hessian H ij is a Hermitian matrix. Let M kij denote the principal minorsof this Hessian. We have : M ij = ( p i. − p i,j ) p i,j p i. ≥ M ij = det " ( p i. − p i,j ) p i,j p i. − p i. − p i,j p i. − p i. − p i,j p i. p i,j p i. = 0 . Moreover, M kij = 0 for all k ≥
3, because then M kij has at least two equal lines(or columns). Therefore, the Sylvester’s criterion is satisfied and each term T ij is convex.Similarly, each term T ij ≡ p i,j ln( p i,j p .j ) is also convex, and hence F is convexas a sum of convex terms.Let us now show, using a counter-example that the ascendency, defined by[15] as A ( f ) = X ij ∈ E f i,j f .. ln (cid:18) f i,j f .. f i. f .j (cid:19) is not convex. In this aim, consider the graph with three flows f , f and f shown on Fig. 6.We compute A = f ln (cid:18) f ( f + f ) ( f + f ) (cid:19) − f ln ( f + f ) − f ln ( f + f ) . The gradient of A is given by A f = − ln ( f + f ) − , A f = − ln ( f + f ) − ,A f = ln (cid:18) f f + ( f + f ) f + f f (cid:19) − . Figure 6: Example of 3-species model ecosystemThere is only one extremum, f = f = 0 and f = e − ≈ . . The Hessian of A is H A = − f + f − f + f − f + f − f + f − f + f − f + f − f − f f f +( f + f ) f + f f f . For f = 1 , f = 1 , f = 0 .
4, we get H = − .
714 0 − . − . − . − . − .
714 1 . , which is clearly not positive semidefinite. Appendix 3: Lemma 8.1
A classical result is that the eigenvalues of a skew-symmetric real matrix arepure imaginary or zero. Indeed, let A be a skew-symmetric matrix and B = iA ,then B ∗ = − iA T = iA = B and therefore B is Hermitian. Since B has all realeigenvalues λ , ..., λ n , all the eigenvalues of A are of the form − iλ , . . . , − iλ n and thus all pure imaginary.The following modified version deserves to be proven. Lemma 8.1
Let J be a matrix such that J = DA , where D is diagonal withnon zero elements and A is skew-symmetric. Then the eigenvalues of J areeither imaginary numbers or zero. proof Let λ be an eigenvalue of J . Let v be an associated eigenvector, such that29 Av = λv . Since D is invertible, we have Av = λD − v. (42)First, the product of both sides of (42) with v T gives v T Av = λv T D − v. (43)Since λv T D − v is a scalar, taking transpose of both sides of (43) yields ( Av ) T v = − v T Av , and hence λv T D − v = ( Av ) T v = − v T Av. (44)Second, the complex conjugate of (42) is Av = λD − v. The product of bothsides with v gives v T Av = λv T D − v. (45)Finally, since v T D − v = ( D − v, v ) = ( D − v, v ) >
0, identifying the twoexpressions of v T Av in (44) and (45) yields λ = − λ, so that λλ