An Exergy-Based Model for Population Dynamics: Adaptation, Mutualism, Commensalism and Selective Extinction

Following the critical analysis of the concept of “sustainability”, developed on the basis of exergy considerations in previous works, an analysis of possible species “behavior” is presented and discussed in this paper. Once more, we make use of one single axiom: that resource consumption (material and immaterial) can be quantified solely in terms of exergy flows . This assumption leads to a model of population dynamics that is applied here to describe the general behavior of interacting populations. The resulting equations are similar to the Lotka-Volterra ones, but more strongly coupled and intrinsically non-linear: as such, their solution space is topologically richer than those of classical prey-predator models. In this paper, we address an interesting specific problem in population dynamics: if a species assumes a commensalistic behavior, does it gain an evolutionary advantage? And, what is the difference, in terms of the access to the available exergy resources, between mutualism and commensalism? The model equations can be easily rearranged to accommodate both types of behavior, and thus only a brief discussion is devoted to this facet of the problem. The solution space is explored in the simplest case of two interacting populations: the model results in population curves in phase space that can satisfactorily explain the evolutionistic advantages and drawbacks of either behavior and, more importantly, identify the presence or absence of a “sustainable” solution in which both species survive.


List of Symbols
ji a Discharge from "j" used by "i" Greek symbols i c Minimum consumption for survival of "i", W i  Dimensionless exergy capture coefficient of "i" E  Exergy flow rate, W i  Inflow capture coefficient, W/person of "i" p i Global discharge coefficient of "i" i  Lyapunov exponents, 1/y r i Limit growth rate (births-deaths) of "i", 1/y i  Intrinsic mortality rate of "i", 1/y w i Exergy discharge coefficient of "i"

Introduction
The concept of sustainability is based on the qualitative or quantitative assessment of environmental interactions between a system and its environment at large, and implies the identification of a certain number of natural-and possibly finite-resources so that the dynamic interplay among species in a given ecological niche can be modeled to study the time evolution of the system.A first problem is thus that of "measuring" the resources, i.e., of finding a quantifier that is at the same time useful, convenient and general enough to represent material and immaterial fluxes.A second problem is posed by the empirical observation that most existing species display a multitude of "modes" of natural interrelation, like indifference, competition, cooperation, antagonism, adaptation, etc.A third problem is the "flexible attitude" exhibited by most species towards the substitutability of certain resources [1][2][3].
Since the very first attempts to construct behavioral and evolutionary models of living organisms, two substantially different approaches have been taken: a) An ad hoc approach that focuses on the particular dynamics of a single species or of a particular, generally small, ecological niche: models of such kind strongly depend on a set of very stringent initial assumptions, but thanks to their "specificity" have enjoyed a remarkable degree of success in reproducing experimental field data and-to a lesser extent-to predict future trends; b) The opposite approach is also possible: one tries to understand the population dynamics in a global sense, emphasizing and exploiting similarities in the behavior of different species in different environmental conditions and within different trophic chains.An original model of this second type has been proposed in previous papers by the present authors [4,5], and is based on the generally accepted principle that population dynamics depend substantially on the availability of primary resource and on the axiom that the consumption of such resources can be measured by an environmental indicator called extended exergy, EE in the following [6].
Extended Exergy is a physical quantity (measured in J, its fluxes in W) that explicitly measures the different forms in which natural resources are "available", has the logical attributes of a "cost", and thus is amenable to a resource accounting procedure -the method in fact goes under the name of Extended Exergy Accounting, EEA.Extended exergy is additive, explicitly includes irreversible losses, and its formulation covers the so-called "externalities": using a somewhat anthropological expression, it can be rightly said that the EE of a commodity is a measure of the biosphere "effort" to convert low-entropy into high-entropy resources while generating that commodity.(Note: The conceptual novelty of Extended Exergy is its inclusion of the equivalent Capital, Labor and Environmental Remediation costs into the exergy "embodied" in a product of whatever material-or energy conversion chain.In the context of this paper, Capital and Labor are of course absent, and the difference between the use of exergy and extended exergy consists in the inclusion in the latter of the amount of natural exergy resources required of the environment to "biodegrade" the environmental damage produced by a species.To provide a simple example, the restoration of the trees destroyed by the Southern Pine Beetle [3] requires a certain amount of cumulative exergy, in the form of solar irradiation and nutrients, that represents the "cost" for the biosphere to recover its previous state. Similarly, the oxidation of pollutants in a water basin requires solar irradiation, microbic action, etc., all of which can be expressed in terms of "exergy inflow").
Nature's intertwined trophic chains are far too complex a system to be amenable to a tractable comprehensive mathematical approach: for example, it is impossible to classify within general schemes voluntary or instinctive behaviors that, especially in superior species, can be treated only by recurring to species-specific, and therefore non-general, assumptions.Reasoning in an ethologic sense, though, it is legitimate to posit [7] the dependence of the dynamics of a certain species on the available resources (measured by extended exergy), formulate a proper mathematical model on its basis, and apply the model to predict the dynamics of a single species (local level) or of an entire ecological niche (global level): since a large amount of experimental data are available on the dynamics of certain populations, the results can subsequently be validated and verified.It was shown in two previous papers [4,5] that such an analysis of the sustainability concept leads to credible results.In the model adopted in those papers, "competition" was among the species and for the available resources.The underlying assumptions are quite simple: 1) The number N of individuals in a species at a given time depends on the global exergy resources a population can avail itself of.2) These resources may be quantified in terms of Extended Exergy (i.e., by their primary resource equivalent), and expressed as a flux of exergy that each species may tap.
The resulting equation for a single population (for the complete derivation, see [4,5]) couples the population size (its numerosity N) with the exergy rate of the global exergy resources exploited by the population in its ecological niche: where the dot denotes the time derivative, the constant r is the limit growth rate of the population (in the case of infinitely abundant flow of resources) and  is the intrinsic mortality rate of the species (calculated in the case in which the availability of resources drops to zero).It is useful to briefly discuss the behavior of the solutions of this equation in a very simple but realistic situation.Suppose that the incoming exergy rate is constant in time: this means the population ) (t N can avail itself of the same exergy inflow at every instant of time.What one expects is then a "logistic" behavior, in the sense that, after an initial transient, the solution ) (t N approaches a limit depending on the parameters of the model and on the available exergy flow.Indeed in this simple case the equation possesses an implicit solution given by: The behavior is indeed similar to that of a logistic-type equation under a carrying capacity constraint: here however the specific exergy consumption rate and the parameters of the model exactly define the carrying capacity.Note that in this case, under a constant influx condition (and therefore an infinite amount of cumulative exergy) at its disposal, the numerosity can sustain itself to the value * N for an infinite time.
Obviously this is no longer true when describes a finite amount of cumulative exergy (see [4,5]).Notice also that the carrying capacity * N is proportional to E  but also inversely proportional to c, a model constant representing the specific exergy consumption rate corresponding to "minimum survival" of the population; indeed if in a given interval of time the total energy consumption rate ) (t E  of the population is much higher than ) (t cN then the population increases in that interval, while in the opposite situation the population would decrease: this means that the curve ) (t N "follows" in some sense the curve ) (t E  .A last remark: in the previous simple situation when ) (t E  is constant we can have two qualitative different behavior of the solution ) (t N (see e.g.[4]): if at 0 t t  the initial population size is larger than * N , then the population will monotonically (and exponentially) decrease to the value * N , while if the population size at 0 t t  is smaller than * N , then the population will increase to the value * N .Some plots of the transient  In the case of Z populations the previous equation generalizes to: where again the r j , j=1…Z, are the limit growth rates of each population and the  j are the intrinsic mortality rate of each species.Now the "source term" represents the exergy rate allocated to population j and, as such, it will depend on the strength of the mutualism and competition between populations i and j.Once the characteristics of each species have been assigned (r j ,  j and a "scenario" is completely identified by specifying the evolution of the incoming total exergy rate: for each scenario, equations (3) describe the evolution of the Z populations.For a detailed analysis on the dynamics of competing populations governed by Equation (3), we refer readers to [4].The relative strength of the interactions between populations i and j, that we describe by means of a set of parameters characteristic of the species and of their ecological niche, determines the relative amount of the global incoming exergy allocated to each population.
In this paper we show that the system (3) is comprehensive and general enough to encompass the three following experimentally observed modes of species interaction [8,9]: a) Indifference: two species in the same environmental niche share the available incoming exergy flux but do not feed on each other's discharges; b) Commensalism: one of the species survives solely on the resources released by its "host"; c) Mutualism; the two species compete for the externally available resources, but feed on each other's discharges.

A Model of Two Interacting Populations Sharing a Common Renewable Resource
Consider two species that share a single renewable resource and who "feed" on their respective discharges (Figure 2).Notice that the relative strength of the interactions between populations i and j, that we describe by means of a set of parameters characteristic of the species and of their ecological niche, determines the relative amount of the global incoming exergy allocated to each population.Our model assumes that, at each instant of time, the rate of exergy accumulation or depletion (Ė 1 and Ė 2 here) is the force that drives the increase (or decrease) in the respective numerosities: it will be therefore regarded as "useful input".The exergy "balances" read: Notice that, since exergy is not a conserved quantity, in this paper the expression "exergy balance" is adopted to denote the closure of the exergy equation obtained by introducing the exergy destruction    We assume that the allocation of the exergy flows 1 E  and 2 E  increases with the numerosity of the respective population.So we set: ,  1 and  2 (both in J/individual/day) being the specific exergy consumptions of the two species (notice that  1 + 2 =1 and as requested by the imposed boundary conditions).In most real situations, the parameters i  could depend also on time: in this paper, for simplicity, we posit that they are constant.
Note that the parameter a nk is the amount of the discharge of species k "recycled" by species n: if a nk =a kn =0, the two populations are indifferent to each other, i.e., neither one exploits the discharges of the other as its own exergy input).Since the instantaneous rate of discharge and destruction of the exergy flow can be only a fraction of the incoming one, we set and where i w and i  are two dimensionless species-specific constants measuring the portion of incoming flow being discharged and destroyed by population i N : In the model presented here, these "discharges" are NOT equivalent to the detritus (which is a pure material flow rate): they include also any form of exergy released by species j (thermal, chemical, etc.) It is clear that this interpretation for the constants i w and i  implies the constraints Equations ( 4) and ( 5) then result in: we see that the above expressions are indeed, at every instant, non-negative, as requested by their physical interpretation.It is convenient to set : by substituting ( 6) and ( 7) into (3), we obtain: Where for simplicity we set . The strong coupling between the numerosity of each species and the time evolution of ALL the participating species is apparent.Notice that in this case, Z = 2, but the formulae can be easily extended to account for whatever number Z of interacting populations, see [1].We shall assume in the following that the exergy influx Ė in is constant and that the model coefficients , w, , r,  and c are known for all species considered.
The type and degree of interaction between the two species is represented by the terms k  and k  that respectively contain the " j " (that express the exergy voracity of species j) and the " ij ," that quantify the degree of mutualism.Let us examine, on the same line as proposed in [8], the three possible general behaviors.Notice that Jørgensen [8] uses "cooperation" and "parasitism" to indicate "mutualism" and "commensalism" respectively.

Indifference
Species 1 and 2 display a disjointed behavior: in spite of sharing the same ecological niche, they do not interact, in the sense that neither species feeds, directly or indirectly, on the waste of the other one, which means that 0 . Equations ( 8) and ( 9) simplify to: Even in this rather simple case, the strong non-linearity of the equations make them not amenable to an explicit solution.However the global behavior of the solutions, as we will see, strongly depends on the two quantities will be discussed at the end of the section.In the state plane ) , ( 21 N N , the system of Equations 8 and 9 possesses three fixed points, listed in Table 1: Table 1.Fixed points of the system of Equations (10-11). 1 The corresponding Lyapunov exponents are given by the following expressions: Table 2. Lyapunov exponents of the system of Equations (10)(11) at points a, b and c.
Point a is always unstable.The exponents 1

 and 2
 for points b and c are such that one point will be unstable and the other stable: the assumption are considered here because this semiplane is invariant under the flow).The three regions are divided by the so called isoclines 1 s and 2 s , defined as (see Figure 3): . Notice that the two lines have the same angular coefficient and, since here  In region I 1 N and 2 N are both increasing functions of time, in region III they are both decreasing functions of time, while in region II 1 N is increasing and 2 N decreasing.The behavior of the solution is characterized as follows (see Appendix 1 for the proof): if at some time 0 t the solution is in region I then there exists a finite time 0 t t  for which the solution will leave that region.If at some 1 t the solution is in region II then it will remain in this region and asymptotically approach point c for t∞.If at some time 3 t the solution is in region III then it will either remain there for all subsequent times, asymptotically approaching c, or it will enter at some time   : The arrows indicate the time direction of the evolution.

Commensalism
One of the two species displays a commensalistic behavior with respect to the other one, in the sense that rather than extracting its exergy feed rate from Ė in , it thrives solely on the "discharge flows" of the other species.Asymmetrically, the species "pestered" by the parasite does not feed on it-(in other words, any degree of mutualistic behavior is disregarded).Assuming species 2 to be the parasite, this implies: Under these circumstances, (10-11) become: Note that now the coupling between the species (or better of species 2 with species 1) is described only in terms of the two constants 1  and 1  : thus, the system (14) can be formally integrated, even if the solutions can be only written out implicitly: where the two values 1 N and 2 N are given by: : increasing the discharged and destroyed portion of the incoming exergy flows proportionally reduces the corresponding carrying capacity value for population 1.But population 2 feeds on this discharged flow: Equation 16shows that indeed the carrying capacity limit for 2 N is proportional to 1 w .Note that, in contrast to the case of competing populations, and assuming 0 1  w , the two populations will always both survive.In fact, species 1 meets no competition in tapping the exergy flow it needs and species 2 feeds on the "crumbs" of species 1.This result is quite obvious, but we stress that: a) it is derived on exact mathematical basis.b) the carrying capacity limits are directly linked to the parameters of the model.c) Equations 14 describe also how fast these limits are reached (Figure 6). Figure 6.An example of the evolution of the two populations in the commensalistic case.

Mutualism
This is the most general case, in which both 12 a and 21 a are strictly positive (>0).The model equations are 6 and 7. Again, it is convenient to examine the phase space ) , ( 21 N N , or better the invariant subspace given by   . In this case, the Lyapunov coefficients assume the following values (the fixed points a, b and c are the same as in Table 1): Table 3. Lyapunov exponents of the system of Equations (10)(11) at points a, b and c. 1 All three points are unstable.This means that the system tends to avoid situations in which one or both of the populations become extinct.This is an indication that in the present model, cooperation (mutualism) may lead to sustainable scenarios.Let us show the existence of another fixed point and look more closely at the dynamics.
If we consider the numerators of (8-9), repeated here for clarity in Equations 16, it is possible to verify that, for whichever value of the model parameters, the two curves described by ( 16) always intersect in a point (the fixed point) on the semi-plane   , which we shall denote as d.
It is useful to write down the solution of this system in terms of the parameter  given by the ratio of where the parameter  solves the following quadratic equation ( the two solutions   being: From (20) we note that   is always positive and   always negative.The negative solution must be discarded because 1 N and 2 N are both positive definite.So there exists another fixed point given by:   Furthermore this point is attractive.Indeed we can label different domains in the semi-plane according to the signs of 1 N  and 2 N  .Consider species 1 first.In order to we need: Equation 22 describes a curve dividing the semi-plane in two regions; region I, where 0 1  N  , and region II, where 0 1  N  . The two lines depicted in Figure 7 correspond to the two cases ) will be discussed at the end of this section.Repeating the same argument for species 2, we find two more curves dividing the semi-plane as in Figure 8, again corresponding to the two cases Note that the fixed point d is given by the intersection of a curve in Figure 7 with a curve in Figure 8.For example if we have , to find the fixed point d we can superimpose the two curves on the right in Figure 7 and Figure 8 to obtain the picture shown in Figure 9.
A plot of trajectories corresponding to initial values in the four different regions is given in Figure 10.It is apparent from Equation 20 that the value of   is independent of in E  : this is a strong constraint on the system because it means that the asymptotic value of the ratio and depends only on the parameters characterizing the species and their interaction (competitioncooperation).From a mathematically point of view this property of the solution is tautologic, but the question of whether a real system actually possesses this attribute is indeed far from obvious: since the ratio 2 1 N N and in E  are easily measurable in field experiments, the above constraint is suggesting a possible experimental "check" on the model.In a future paper we plan to examine some of the available experimental result in order to verify whether this statement is falsifiable.
Let us now briefly discuss the case . Looking at Figure 7 now the two points marked with black dots on the N 1 axis coincide, so that the curve dividing the semi-plane into two regions becomes a straight line, and the intersection point is now given by: where   is given by Equation 20.In this case the equation for 1 N decouples, since we have: This is the same equation as the first of (14).In that case species 1 (the host), having no competitors, fed on all of the available exergy flow, so that after discharges and destructions are taken into account, its net exergy income was (see Equation 4).On the contrary, in this case, 1 N captures only a fraction of the environmental exergy at its disposal, the complementary fraction being taken by species 2; but there is another "resource" flow coming from the fraction of the exergy of species 2 discharge "recycled" from species 1 and these two flows, because of . This is a rather special situation but it is interesting to note how the equation for 1 N decouples.In contrast to the second of (15), however, the equation for 2 N  depends on both populations, so we cannot give analytical results as in that case.It is clear however that the qualitative dynamics cannot differ much from the general case.

Conclusions
The first conclusion that can be derived from the arguments developed in this paper is that the axiom that is at the foundation of our reasoning, namely that resource consumption (material and immaterial) can be completely and unambiguously quantified in terms of exergy flows, is a convenient and fruitful working hypothesis.The model of population dynamics derived from such axiom results in strongly coupled and intrinsically non-linear evolution equations.In this paper, the time evolution patterns of two species sharing an ecological niche have been studied under three different scenarios: indifference, commensalism and mutualism.Our results indicate that: if the two species just coexist in the same niche without cooperating or preying on each other, their competition invariably leads (Section 2a) to the extinction of one of them, unless their genetic (r,), consumption (c, ) and efficiency (p) parameters are in a very precise and delicately tuned balance, in which case the two species may coexist sustainably; if one of the two species assumes a commensalistic type of behaviour (it feeds on the other's discharges, Section 2b), then both species reach their respective carrying capacity limits; if the two species display a mutualistic type of behaviour (each one feeds partially on the other's discharges, Section 2c), then there is always a stable sustainable point and both populations survive.
In all cases, the sustainability of the point in state-space of the surviving population(s) prescribes a certain well-identified numerosity.
In spite of the apparent success of our model, caution ought to be exercised in generalizing its results: nature is far more complex than any complicated and non-linear model makes it appear, and the results presented here are but a basis-and an incentive!-forfurther study.In particular, the following problems need to be carefully addressed if this type of approach is brought down from a descriptive to an applicative or predictive level: 1) Substitutability of exergy resources: our model explicitly lumps all of the "available" material and energy fluxes into the extended exergy inflow in E  .This poses no problem for the energy flows, but contains the implicit assumption that ALL material flows are perfectly substitutable.
In the simplest case, consider an animal species that has access to two or three different sources of vegetable food: obviously, each vegetable has a different nutritive value for the species, but also a different extended exergy content.Thus, our model would predict that the "optimal choice" for the animal would be to feed on the lower-embodied exergy food with the highest nutritive value.But this prediction disregards preferences, spatial distribution, physical accessibility, and similar factors.For more complex behaviours (for instance, related to the human species), the choice between, say, bronze and iron to build spears is even more multi-faceted than implicitly posited by the model.2) Actual possibility of species reaching a stable stationary state: this is of course an extremely important point, because if two species in a single ecological niche reach a stable state (point d in Figure 10), this is by definition a sustainable situation.For larger number of species and realistic conditions, the problem of the stability of such stationary states is of paramount importance, and has not been addressed here.3) Metadynamics: the model presented in this paper is spatially lumped and the niche in which the species thrive has "impermeable boundaries".Immigration and emigration are completely neglected, as well as spatially varying distributions of in E  .While the former could be modeled by inserting an "immigration term" in Equations ( 6) and (7), the spatially lumped nature of the model is difficult to modify.One solution might be that of "discretizing" the domain, and apply the equations separately to different areas, accounting for the variation of in E  in each subdomain.This would make our model somewhat comparable to the so-called "patch" [10] or "incidence function" models [11], but would also demand for very taxing numerical simulations, not considered here.4) Limiting factors ("exogenous factors" in [3]): while we maintain that (extended) exergy is a satisfactory indicator of species numerosity (in the sense that high exergy niches must logically and thermodynamically correspond to globally higher numerosities), it must be remarked that in real ecological niches there are phenomena occurring at some scales that affect the numerosity of species at a different scale [12].Viral and fungal infections to which species N is insensitive may affect one or more of its food sources; growth of high-foliage trees with which a vegetable species does not compete for nutrients may reduce its share of solar irradiation, etc.: at the present state of development it is not clear whether and how our model may accommodate such effects (like in "tritrophic" or "multi-trophic" models [3]).5) Though not addressed in this paper, the application of our model raises a question that could be of paramount importance for practical applications: is the efficiency of one single species in the niche compatible with the overall energy-conversion efficiency of the same niche, considered as a lumped system?In other words, if two or more species reach a sustainable point ("d" in Figure 10), do their respective numerosities satisfy both a local (for each species) and a global (for the entire niche) "optimal exergy use" criterion?This is a topic we have left for a future study. 1  N has a minimum for   t t .However this cannot be true because 1  N is an increasing function in region II.So the hypothesis that the solution crosses the line , but again this cannot be the case since 2  N is a decreasing function in region II.


Exergy flow rate discharged by "i", W i Exergy destruction coefficient of "i" i E   Exergy flow rate destroyed by "i", W i Coupling terms in Equation (7), W/person N Population numerosity  Parameter in Equation (15)

Figure 1 :
Figure 1: the gray curves represent the first case, the black ones the latter.

Figure 1 .
Figure 1.Qualitative transient behavior or the solutions of Equation (3).

.For
These two parameters contain the efficiency of each species to gain "net" exergy flows ( k k p  ), their "minimal survival" requests ( k c ) and the ratio of the growth rate to the mortality rate k k r in the respective limit cases of infinite and zero resources available.
c being asymptotically stable and point b unstable.We can divide the plane )

Figure 3 .
Figure 3.The three regions defined by the two isoclines.The fixed points are marked by a dot.
different initial conditions, is shown in Figure4.

Figure 4 .
Figure 4. Three different orbits for the system (10-11) in the case

Figure 7 .
Figure 7.The two regions identified by the curve (22) for 1 2 2 p    (left) and

Figure 9 . 1 N  and 2 N
Figure 9.The fixed point d as the intersection of the two curves in Figures 7-8.

Figure 10 .
Figure 10.Four different orbits for the general case.

2 1 N
N is also independent of in E  1s must be false.Suppose instead that it crosses the line2s at