## Abstract

We carry out mathematical analyses, à la Helmholtz’s and Boltzmann’s 1884 studies of monocyclic Newtonian dynamics, for the Lotka–Volterra (LV) equation exhibiting predator–prey oscillations. In doing so, a novel ‘thermodynamic theory’ of ecology is introduced. An important feature, absent in the classical mechanics, of ecological systems is a natural stochastic population dynamic formulation of which the deterministic equation (e.g. the LV equation studied) is the infinite population limit. Invariant density for the stochastic dynamics plays a central role in the deterministic LV dynamics. We show how the conservation law along a single trajectory extends to incorporate both variations in a model parameter *α* and in initial conditions: Helmholtz’s theorem establishes a broadly valid conservation law in a class of ecological dynamics. We analyse the relationships among mean ecological activeness *θ*, quantities characterizing dynamic ranges of populations *α*, and the ecological force *F*_{α}. The analyses identify an entire orbit as a stationary ecology, and establish the notion of an ‘equation of ecological states’. Studies of the stochastic dynamics with finite populations show the LV equation as the robust, fast cyclic underlying behaviour. The mathematical narrative provides a novel way of capturing long-term dynamical behaviours with an emergent *conservative ecology*.

## 1. Introduction

In the mathematical investigations of ecological systems, conservative dynamics are often considered non-robust, and thus unrealistic as a faithful description of reality [1,2]. Through recent studies of stochastic, nonlinear population dynamics, however, a new perspective has emerged [3–5]: the stationary behaviour of a complex stochastic population dynamics almost always exhibits a divergence-free cyclic motion in its phase space, even when the corresponding system of differential equations has only stable, node-type fixed points [6]. In particular, it has been shown that an underlying volume-preserving conservative dynamics is one of the essential keys to understanding the long-time complexity of such stochastic systems [7–9].

The aim of this work, following a proposal in [9], is to carry out a comprehensive stochastic dynamic and thermodynamic analysis of an ecological system with sustained oscillations. In the classical studies on statistical mechanics, developed by Helmholtz, Boltzmann, Gibbs and others, the dynamical foundation is a Hamiltonian system [10–12]. The theory in [9,13] generalized such an approach that requires no *a priori* identification of variables as position and momentum; it also suggested a possible *thermodynamic structure* which is purely mathematical in nature, independent of Newtonian particles. In the context of population dynamics, we shall show that the mathematical analysis yields a *conservative ecology*.

Among ecological models, the Lotka–Volterra (LV) equation for a predator–prey system has played an important pedagogical role [1,2,14], even though it is certainly not a realistic model for any engineering applications. We choose this population system, in this work, because of its mathematical tractability, and its stochastic counterpart in terms of a birth and death process [15,16]. It can be rigorously shown that a smooth solution to the LV differential equation is the law of large numbers for the stochastic process [17]. In biochemistry, the birth–death process for discrete, stochastic reactions corresponding to the mass-action kinetics has been called a Delbrück–Gillespie process [3].

The LV equation in non-dimensionalized form is [1]
*x*(*t*) and *y*(*t*) represent the populations of a prey and its predator, each normalized with respect to its time-averaged mean populations. The *xy* term in *f*(*x*,*y*) stands for the rate of consumption of the prey by the predator, and the *αyx* term in *g*(*x*,*y*;*α*) stands for the rate of prey-dependent predator growth.

It is easy to check that the solutions to (1.1) in phase space are level curves of a scalar function [1]
*Γ*_{H=h} to denote the solution curve *H*(*x*,*y*)=*h*, and *Γ*_{H=h}. Figure 1 shows the contours of *H*(*x*,*y*) with *α*=1 and *H*(*x*,*y*)=2.61 with different *α*s.

Let *τ* be the period of the cyclic dynamics. Then, it is easy to show that [1]
*Γ*_{H=h}, using the Lebesgue measure in the *xy*-plane. The appropriate measure for computing the area will be further discussed in §2. The parameter *α* represents the relative temporal variations, or dynamic ranges, in the two populations: the larger the *α*, the greater the temporal variations and range in the predator population, and the smaller in the prey population.

The paper is organized as follows. In §2, an extended conservation law is recognized for the LV system. Then, the relationship among three quantities—the ‘energy’ function *H*(*x*,*y*), the area *Γ*_{H=h} and the parameter *α*—is developed. According to the Helmholtz theorem, the conjugate variables of *α* are found as time averages of certain functions of population *x*(*t*) and *y*(*t*). Analysis on those novel ‘state variables’ demonstrates the tendency of change in mean ecological quantities such as population range or ecological activeness when the parameter *α* or energy *H* varies. In §3, we show that the area *Γ*_{H=h} is related to the concept of entropy. In §4, the conservative dynamics is shown to be an integral part of the stochastic population dynamics, which necessarily has the same invariant density as the deterministic conservative dynamics. In the large population limit, a separation of time scale between the fast cyclic motion on *Γ*_{H=h} and the slow stochastic crossing of *Γ*_{H=h} is observed in the stochastic dynamical system. The paper concludes with a discussion in §5.

## 2. The Helmholtz theorem

Equation (1.1) is not a Hamiltonian system, nor is it divergence free,
*G*(*x*,*y*)=*xy*. One can, in fact, understand this scalar factor as a ‘local change of measure’, or time *G*^{−1}(*x*,*y*) is an invariant density of the Liouville equation for the deterministic dynamics (1.1), and, more importantly, the invariant density of the Fokker–Planck equation for the corresponding stochastic dynamics. As will be demonstrated in §2b,c, the statistical average of quantities according to the invariant measure *G*^{−1}(*x*,*y*) d*x* d*y* can be calculated through the time average of those quantities along the system’s instantaneous dynamics. Knowledge about the system’s long-term distribution is not needed during the calculation. These facts make *G*^{−1}(*x*,*y*) the natural measure for computing the area

Any function of *H*(*x*,*y*), *ρ*(*H*) is conserved under the dynamics, as is guaranteed by the orthogonality between the vector field of (1.1) and gradient ∇*ρ* [9],

### (a) Extending the conservation law

The nonlinear dynamics in (1.1), therefore, introduces a ‘conservative relation’ between the populations of predator and prey according to (1.2). If we call the value *H*(*x*,*y*) an ‘energy’, then the phase portrait in figure 1*a* suggests that the entire phase space of the dynamical system is organized according to the value of *H*. The deep insight contained in the work of Helmholtz & Boltzmann [10,12] is that such an energy-based organization can be further extended for different values of *α*: therefore, the energy-based organization is no longer limited to a *single* orbit, nor a *single* dynamical system, but rather encompasses the entire class of parametric dynamical systems. In the classical physics of Newtonian mechanical energy conservation, this yields the mechanical basis of the fundamental thermodynamic relation as a form of the first law, which extends the notion of energy conservation far beyond mechanical systems [18,19].

More specifically, we see that the area, *α* and initial energy value *h*=*H*(*x*(0),*y*(0),*α*). Therefore, there must exist a bivariate function *x*,*y*) which is continuously varying with time; rather it reflects the geometry of an entire orbit. Then, equation (2.4) implies that any such ecological state has an ‘*h*-energy’, *if* one recognizes a geometric, state variable

Equation (2.4) can be written in a differential form
*h*-energy for an ecological system with fixed *α* via the factor *α*-force’ corresponding to the parameter *α*. In classical thermodynamics, the latter is known as an ‘adiabatic’ process.

The Helmholtz theorem expresses the two partial derivatives in (2.5) in terms of the dynamics in equation (1.1).

### (b) Projected invariant measure

For canonical Hamiltonian systems, the Lebesgue measure is an invariant measure in the whole phase space. On the level set *Γ*_{H=h}, the projection of the Lebesgue measure, called the Liouville measure, also defines an invariant measure on the submanifold. If the dynamics on the invariant submanifold *Γ*_{H=h} is ergodic, the average with respect to the Liouville measure is equal to the time average along the trajectory starting from any initial condition (*x*_{0},*y*_{0}) satisfying *H*(*x*_{0},*y*_{0})=*h*.

As we shall show below, the invariant measure for the LV system (1.1) in the whole phase space is *Γ*_{H=h} can be found by changing (*x*,*y*) to intrinsic coordinates (*h*,ℓ):
*Γ*_{H=h}; and *Γ*_{H=h}.

It is worth noting that d*μ*=d*t* on the level set *Γ*_{H=h}. Because the dynamics on *Γ*_{H=h} is ergodic, the average of any function *ψ*(*x*,*y*) under the projected invariant measure on *Γ*_{H=h} is equal to its time average over a period

### (c) Functional relation between ln A , *α* and *h*

Under the invariant measure *G*^{−1}(*x*,*y*), the area *Γ*_{H=h} is
*τ*(*h*,*α*) is the time period for the cyclic motion. Furthermore,
^{t} is the time average, or phase space average according to the invariant measure. We can also find the derivative of the area *Γ*_{H=h} with respect to the parameter of the system *α* as

In this setting, the Helmholtz theorem reads
*θ*(*h*,*α*) here is the mean *α*-force is then defined as

It is important to note that the definition of *F*_{α} given in the right-hand side of (2.20) is completely independent of the notion of *F*_{α}(*h*,*α*) is a function of both *h* and *α*, however. Therefore, the value of *α*-work *F*_{α}(*h*,*α*)*dα* depends on how *h* is constrained: there are iso-*h* processes, iso-*θ* processes, etc. [20].

### (d) Equation of state

The notion of an equation of state first appeared in classical thermodynamics [18,19]. From a modern dynamical systems standpoint, a fixed point as a function usually is continuously dependent upon the parameters in a mathematical model, except at bifurcation points. Let *α* be a parameter, then the function *equation of state* for the long-time ‘equilibrium’ behaviour of the dynamical system.

If a system has a globally asymptotically attractive limit set that is not a simple fixed point, then every geometric characteristic of the invariant manifold, say *α*. In this case,

The situation for a conservative dynamical system with centre manifolds is quite different. In this case, the long-time behaviour of the dynamical system, the foremost, is dependent upon its initial data. An equation of state therefore is a functional relation among (i) geometric characteristics of a centre manifold *α*, and (iii) a new quantity, or quantities, that identifies a specific centre manifold, *h*. This is the fundamental insight of the Helmholtz theorem.

In ecological terms, the area under the invariant measure, *α*, on the other hand, is the proportion of the predator’s over the prey’s population ranges of time variations:
*θ* can be viewed as a measure of the mean ecological ‘activeness’,
*x* and *y*, respectively, to the fixed populations in equilibrium (1,1). For population dynamic variable *u*, equation (2.22) suggests a norm *θ*=〈*α*∥*x*−1∥〉^{t}=〈∥*y*−1∥〉^{t}, and an averaged norm of *per capita* growth rates in the two species,
*α*. In other words, when |*F*_{α}| is greater, more *h*-energy change is needed to vary *α*. It is also worth noting that |*F*_{α}| is positively related to the prey’s average population range. In fact, we can define another ‘distance’ of the prey’s population *x* to 1 as: *F*_{α}=−〈∥*x*−1∥_{F}〉^{t}−1. Note that for very small *u*: ∥*u*∥≈*u*^{2}≈2∥*u*∥_{F}.

Figure 2 shows various forms of ‘the equation of state’, e.g. relationships among the triplets (*α*,|*F*_{α}|=−*F*_{α},*h*), (|*F*_{α}|,*θ*,*h*) and (*α*,*θ*,*h*) in the first row; among the triplet (*α*,|*F*_{α}|,*θ*) in the second row; and among the triplet *α*,|*F*_{α}|,*θ*) is just like that among (*V*,*P*,*T*) in the ideal gas model: the mean ecological activeness *θ* increases nearly linearly with the ecological force |*F*_{α}| for constant *α*, and with the proportion *α* of the predator’s population range over the prey’s, for constant |*F*_{α}|; ecological force |*F*_{α}| and the proportion *α* of population ranges are inversely related under the constant mean activeness *θ*. And when *θ*=0, *α*(*F*_{α}+1)=0. Other features can be observed by looking into the details of each column.

Figure 2*a*,*d*,*g* demonstrates that, as the proportion *α* of the population ranges increases, the ecological force |*F*_{α}| is alleviated (for a given *h*-energy or ecological activeness *θ*). This is due to the positive relationship between the ecological force |*F*_{α}| and the prey’s population range (as shown in equation (2.24)). Because *α* is the proportion of the predator’s population range over the prey’s, |*F*_{α}| and *α* would be inversely related when any resource, *h*-energy or activity, *θ*, remains constant. This fact means that on an iso-*h* or iso-*θ* curve, when the proportion *α* is large, relatively less *h*-energy change is needed to reduce it. Figure 2*a*,*d*,*g* also demonstrates an inverse relationship between *α* and the total population range *θ*, which reflects the fact that, as the proportion of the predator’s population range over the prey’s increases, the total population range of the two species would actually decrease.

Figure 2*b*,*e*,*h* demonstrates that the ecological force |*F*_{α}| and the total population range *θ* (with a given *h*-energy or *α*). This observation means that it would also take more *h*-energy to change the proportion *α* of the predator’s population range over the prey’s, if mean ecological activeness rises, and that a greater population range would be explored with more ecological activeness *θ*.

Figure 2*c*,*f*,*i*, which shows the relation between *θ* and *α*, is interesting: under constant *h*-energy, as the proportion *α* of the population ranges increases, the ecological activeness *θ* decreases, in accordance with the drop in the total population range *F*_{α} is to remain constant, the ecological activeness *θ* actually increases with *α*. This means that, under constant resource (*h*-energy), the proportion *α* of the predator’s population range over the prey’s restricts the mean ecological activeness. But if we fix the ecological force or total population range (supplying more *h*-energy), an increase in the predator’s population range over the prey’s can increase ecological activeness.

## 3. Liouville description in phase space

The nonlinear dynamics described by equation (1.1) has a linear, first-order partial differential equation (PDE) representation
*ρ*(*H*(*x*,*y*)) is a stationary solution to equation (3.2), it is not a stationary invariant density to (3.1).

This is due to the fact that vector field (*f*,*g*) is not divergence free, but rather as in (2.1) the scalar factor *G*(*x*,*y*)=*xy*. Then, it is easy to verify that *G*^{−1}(*x*,*y*)*ρ*(*H*(*x*,*y*)) is a stationary solution to (3.1),

### (a) Entropy dynamics in phase space

It is widely known that a volume-preserving, divergence-free conservative dynamics has a conserved entropy *G*(*x*,*y*), the Shannon entropy should be replaced by the relative entropy with respect to *G*^{−1}(*x*,*y*) (see the electronic supplementary material, appendix B, for a detailed calculation):
*canonical conservative* with respect to *G*^{−1}(*x*,*y*) in [8]. In classical statistical physics, the term

We can, in fact, show a stronger result, with arbitrary differentiable *Ψ*(⋅) and *ρ*(⋅) over an arbitrary domain *u*(*x*,*y*,*t*) is observed that determines whether a system is invariant—not the initial data *u*(*x*,*y*,0) [9].

### (b) Relation between A , Shannon entropy and relative entropy

Because a ‘state’ is defined as an entire orbit, it is natural to change the coordinates from (*x*,*y*) to (*h*,*s*) according to the solution curve to (1.1), where we use *s* to denote time, 0≤*s*≤*τ*(*h*,*α*). We have
*Ω*_{B}(*h*) is known as Boltzmann’s entropy in classical statistical mechanics. We see that the *ρ*=*Ψ*=1, and *u*(*x*,*y*,0)=*G*^{−1}(*x*,*y*). Then, *ρ*(*h*)=e^{−h/θ}.

The dynamics (1.1) is not ergodic in the *xy*-plane; it does not have a unique invariant measure, as indicated by the arbitrary *ρ*(*H*) in equation (3.3). However, the function *G*(*x*,*y*), as indicated in equations (2.3) and (3.8), is the unique invariant measure on each ergodic invariant submanifold *Γ*_{H=h}. It is non-uniform with respect to the Lebesgue measure. On the ergodic invariant manifold *Γ*_{H=h}: *u*(*t*)/d*t*=*r*(*t*)*u*(*t*). The regular time average of the *per capita* growth rate is
*G*(*x*(*t*),*y*(*t*)) weighted *per capita* growth rate or ‘kinetic energy’

## 4. Stochastic description of finite populations

Here, we show that the conservative dynamics in (1.1) is an emergent caricature of a robust stochastic population dynamics. This material can be found in many texts [17]. But for completeness, we shall give a brief summary.

Assume that the populations of the prey and the predator, *M*(*t*) and *N*(*t*), reside in a spatial region of size *Ω*. The discrete stochastic population dynamics follows a two-dimensional, continuous time birth–death process with transition probability rate

For a very large *Ω*, the population *densities* at time *t* can be approximated by continuous random variables as *X*(*t*)=*Ω*^{−1}*M*(*t*) and *Y* (*t*)=*Ω*^{−1}*N*(*t*). Then, equation (4.3) becomes a PDE by setting *x*=*m*/*Ω*, *y*=*n*/*Ω* and *u*(*x*,*y*,*t*)=*p*_{m,n}(*t*)/*Ω*:
*ϵ*=*Ω*^{−1}, we can perform the Kramers–Moyal expansion to obtain
**F**(*x*,*y*)=(*f*(*x*,*y*),*g*(*x*,*y*;*α*))^{T} and symmetric diffusion matrix

Equation (4.4) should be interpreted as a Fokker–Plank equation for the probability density function *X*(*t*),*Y* (*t*)) following the Itō integral [15–17]:

### (a) Potential–current decomposition

It can be verified that the stationary solution to (4.4) is actually *G*^{−1}(*x*,*y*)=(*xy*)^{−1}, which is consistent with the discrete case (cf. equation 4.1), and also a stationary solution to the Liouville equation (3.1).

As suggested in [7,9], the right-hand side of equation (4.4) has a natural decomposition:
*ϵ*≠0) is unstable when *x*,*y*>0. The system behaves as an unstable focus as shown in figure 3. The eigenvalues at the fixed point (1+*ϵ*,1−*ϵ*) are

On the other hand, the potential–current decomposition reveals that the system (1.1) will be structurally stable in terms of the stochastic model: any perturbation of the model system will yield corresponding conserved dynamics close to (1.1). The conservative ecology is a robust emergent phenomenon.

Equations such as (4.5) and (4.7) are not mathematically well defined until a precise meaning of integration
*X*(*t*) whose corresponding probability density function *f*_{X(t)}(*x*,*t*) follow different linear PDEs. The fundamental solution to any PDE, however, provides a Markov transition probability; there is no ambiguity at the PDE level. On the other hand, the only interpretation of (4.8) that provides a Markovian stochastic process that is non-anticipating is that of Itō’s [23]. The differences in the interpretations of (4.8) become significant only in the modelling context, when one’s intuition expects that *E*[*X*(*t*)]=0 even for interpretations other than Itō’s.

### (b) The slowly fluctuating *H*_{t}=*H*(*X*(*t*),*Y* (*t*))

With the (*X*(*t*),*Y* (*t*)) defined in (4.4) and (4.5), let us now consider the stochastic functional
*ϵ*, this suggests a separation of time scales between the cyclic motion on *Γ*_{H=h} and a slow, stochastic level crossing *H*_{t}. The method of averaging is applicable here [24,25]:
*ψ*(*x*,*y*)〉^{ΓH=h}=〈*ψ*(*x*,*y*)〉^{t} denotes the average of *ψ*(*x*,*y*) on the level set *Γ*_{H=h}. Then, using the Itō integral, the distribution of *H*_{t} follows a Fokker–Planck equation:
*H* under different *α*s are shown in figure 4. The steady-state distribution *p*^{ss}(*H*) does not depend on the volume size *Ω*=*ϵ*^{−1}.

When *H* is big enough, *p*^{ss}(*H*) increases with *H* without bound, since *b*(*H*) is a positive increasing function. Hence, *p*^{ss}(*H*) is not normalizable on the entire *A*(*H*) approaches zero when *H* approaches *α*+1. Consequently, the absorbing effect at *H*=*α*+1 makes *p*^{ss}(*α*+1) another possible local maximum.

## 5. Discussion

It is usually an obligatory step in understanding an ODE *x** as an implicit function of the parameters (*α*,*β*) [1]. One of the important phenomena in this regard is the Thom–Zeeman catastrophe [1,26]. From this broad perspective, the analysis developed by Helmholtz and Boltzmann in 1884 is an analysis of the geometry of a ‘non-constant but steady solution’, as a function of its parameter(s) and initial conditions. In the context of LV equation (1.1), the geometry is characterized by the area encircled by a periodic solution, *Γ*_{H=h}, where *h* is specified by the initial data: *H*, equation (5.1) can, and should be, interpreted as an extended *H* conservation law, beyond the dynamics along a single trajectory that includes variations both in *α* and in *h*. The partial derivatives in (5.1) can be shown as time averages of ecological activeness *α*, conserved quantity *H* and encompassed area

For the monocyclic LV system, the dynamics are relatively simple. Hence, the state variables have monotonic relationships, the same as those observed in ideal gas models. When the system’s dynamics become more complex (e.g. have more than one attractor, Hopf bifurcation), relations among the state variables will reflect that complexity (e.g. develop a cusp, exhibiting a phase transition in accordance [26]).

When the populations of predator and prey are finite, the stochastic predator–prey dynamics is unstable. This fact is reflected in the non-normalizable steady-state distribution *G*^{−1}(*x*,*y*) on

Indeed, all ecological population dynamics can be represented by birth–death stochastic processes [17]. Except for systems with detailed balance, which rarely holds true, almost all such dynamics have underlying cyclic, stationary conservative dynamics. This work shows that a hidden conservative ecological dynamics can be revealed through mathematical analyses. To recognize such a conservative ecology, however, several novel quantities need to be defined, developed and to become a part of ecological vocabulary. This is the intellectual legacy of Helmholtz’s and Boltzmann’s mechanical theory of heat [28].

## Data accessibility

There are no supporting data in this paper.

## Authors' contributions

Y.-A. M. and H. Q. contributed equally to this work.

## Competing interests

We declare we have no competing interests.

## Funding

There is no funder to report for this work.

## Acknowledgements

We thank R.E. O’Malley, Jr. and L.F. Thompson for carefully reading the manuscript and many useful comments.

- Received July 2, 2015.
- Accepted October 2, 2015.

- © 2015 The Author(s)