Standard Lagrangians in Population Modeling and their Ecological Implications

1. Introduction

In modern theoretical physics, all fundamental equations describing inorganic matter are derived by using the Lagrangian formalism [1,2,3], which requires prior knowledge of functions called Lagrangians [4,6]. A number of different methods have been proposed [7,8,9,10,11,12] to obtain the Lagrangians for the most basic equations of modern classical [5] and quantum [13] physics.

There have been attempts to formulate mathematical models of some biological systems, thereby establishing mathematical oriented theoretical biology [14]. Different areas of mathematics have become increasingly important in biology in recent decades, specifically, statistics in experimental design, pattern recognition in bioinformatics, and mathematical modeling in evolution, ecology, and epidemiology [15]. However, as pointed out in [15] some of these attempts can be classified as ’uses’ but others must be considered as ’abuses’ because most biological systems are complex many-body dynamical systems [16].

Among a large variety of biological systems, the population dynamics plays a special role since it is the key to understanding the relative importance of competition for resources and predation in complex communities, and for preserving biodiversity [16,17]. Population dynamics models that describe interacting species are typically expressed by ordinary differential equations (ODEs), which are first-order, coupled, damped, and nonlinear [16]. Despite the presence of damping and nonlinearities in such models, no clear demonstration of the onset of chaos has yet been shown [18]. However, some studies suggested that insect population dynamics can undergo transitions between stable and chaotic phases for models near a transition point between order and chaos [19].

In theoretical biology, Kerner [20] was the first who applied the Lagrangian formalism to biology and obtained Lagrangians for several selected biological systems described by first-order ordinary differential equations (ODEs). Later, Paine [21] investigated the existence and construction of Lagrangians for similar sets of ODEs following the original work of Helmholtz [6]. First specific applications of the Lagrangian formalism to population dynamics were done by Trubatch and Franco [22], who obtained in an ad hoc manner Lagrangians for the Lotka-Volterra [23,24], Verhulst [25], Gompertz [26] and Host-Parasite [27] population dynamics models. These ad hoc found Lagrangians for the models were formally derived by Nucci and Tamizhmani [28], who used the method based on Jacobi Last Multiplier [8].

The previously found Lagrangians [22,28] must be considered as the generating functions, whose specific forms allow deriving the equations of motion for the population dynamics models by substituting these Lagrangians into the Euler-Lagrange (E-L) equation. The obtained Lagrangians do not represent the difference between the kinetic and potential energies [22], therefore, they must be classified as non-natural [6] or non-standard Lagrangians [9,10,11] to make distinctions between them and the standard Lagrangians (SLs), in which the kinetic and potential energy-like terms are easily identified [1,2,3,4,5].

The main goal of this paper is to develop the Lagrangian formalism for six population dynamics models, which include the above five models as well as the SIR model [29]. The formalism is based on standard Lagrangians and a method to derive these Lagrangians is developed and presented. The emphasis on the SLs makes our approach and obtained results different from those previously found [22,28]. The SLs derived in this paper allow us to identify the kinetic and potential energy-like terms in the population dynamics models, and discuss novel roles of these terms in the models. New and interesting ecological and environmental implications of the derived SLs on the considered population dynamics models are also presented and discussed. The developed method to derive SLs can be used to find SLs for other physical and natural science systems described by second-order ODEs.

The paper is organized as follows. Section 2 presents a brief overview of the Lagrangian formalism and standard Lagrangians; in Section 3, the models of the population dynamics are described and the obtained results are reported and discussed; in Section 4, new ecological and environmental implications are presented and discussed, and Section 5 concludes the paper.

2. Lagrangian formalism

2.1. Action and its Lagrangian function

The functional $\mathcal{S}\left[x\right(t\left)\right]$ is called action and is defined by an integral over a scalar function L that depends on differentiable function $x\left(t\right)$ that describes the evolution of a property of any dynamical system (given by x) in time (represented by t), and on the time derivative $\dot{x}\left(t\right)=dx/dt$. The function $L(\dot{x},x)$ is called the Lagrangian function or simply Lagrangian and, in general, it may also depend explicitly on t and written as $L(\dot{x},x,t)$.

According to the principle of least action, or Hamilton’s principle [2,3,5,6], the action $\mathcal{S}\left[x\right(t\left)\right]$ must obey the following requirement $\delta \mathcal{S}=0$, which guarantees that the action is stationary, or has either a minimum, maximum, or saddle point. The necessary condition that $\delta \mathcal{S}=0$, is known as the Euler–Lagrange (E–L) equation, whose operator $\widehat{EL}$ is defined as

$$\widehat{EL}\equiv \frac{d}{dt}\left(\frac{\partial }{\partial \dot{x}}\right)-\frac{\partial }{\partial x},$$

(1)

and it acts on the Lagrangian, $\widehat{EL}\left[L(\dot{x},x,t)\right]=0$ giving a second-order ODE that becomes an equation of motion for a dynamical system whose $L(\dot{x},x,t)$ is known. The process of deriving the equation of motion is called the Lagrangian formalism [2,3,5], and it has been extensively used in theoretical physics to derive its fundamental equations [2,3,13].

For dynamical systems whose total energy is conserved, the existence of Lagrangians is guaranteed by the Helmholtz conditions [4], which can also be used to obtain Lagrangians. The procedure for finding Lagrangians is called the inverse (or Helmholtz) problem of the calculus of variations [5] and it shows that there are three separate classes of Lagrangians, namely, standard [2,3,5], non-standard [6,8,9,10,11] and null [7,30,31,32] Lagrangians. Both standard and non-standard Lagrangians give the same equations of motion after they are substituted into the E-L equation, or $\widehat{EL}\left[L(\dot{x},x,t)\right]=0$. However, null Lagrangians satisfy the E-L equation identically and therefore they do not give any equation of motion.

In this paper, we establish the Lagrangian formalism for ODEs that describe the time evolution of different models of the population dynamics and derive their standard Lagrangians. We also discuss new biological insights of these Lagrangians into the population dynamics.

2.2. Standard Lagrangians

The Lagrangian formalism requires prior knowledge of a Lagrangian. In general, there are no first principle methods to obtain Lagrangians, which are typically presented without explaining their origin. In physics, most dynamical equations were established first, and only then their Lagrangians were found, often by guessing. Once the Lagrangians are known, the process of finding the resulting dynamical equations is straightforward and it requires substitution of these Lagrangians into the E-L equation. There has been some progress in deriving standard and non-standard Lagrangians for physical systems described by different ODEs (e.g., [8,9,10,11,12]). However, the ODEs considered in this paper are more general than those previously studied, therefore, our paper presents a new approach to solve this problem, and to develop a method that can be used to find standard Lagrangians for the considered population dynamics models. The developed method can also be used to obtain SLs for other physical and natural science systems described by second-order ODEs.

As originally shown by Lagrange [1], Lagrangians for one-dimensional dynamical systems represent the difference between the kinetic and potential energy of these systems, and they can be written as $L(\dot{x},x)={\dot{x}}^2/2-V\left(x\right)$, with $V\left(x\right)$ being the potential energy [2,3]. It is common to call such Lagrangians standard, and this paper is mainly devoted to SLs, which are derived for the population dynamics models by modifying and extending the previously developed method [33].

The Lagrangian formalism based on standard Lagrangians have been well-established in most fields of modern physics [2,3,5,13]. There have also been several attempts to establish the Lagrangian formalism in biology and ecology, specifically in population dynamics [20,21,22,34]. However, all obtained Lagrangians were non-standard. Hence, in this paper, we derive standard Lagrangians and use them to gain new biological insights into the considered population models.

2.3. Method to derive standard Lagrangians

The main objective of this paper is to solve the inverse (Helmholtz) problem of calculus of variations [3,4,5] and derive the standard Lagrangian for a given second-order ODE. Let us consider the following general ODE

$$\ddot{x}+\alpha \left(x\right){\dot{x}}^2+\beta \left(x\right)\dot{x}+\gamma \left(x\right)x=C_0,$$

(2)

where $\alpha \left(x\right)$, $\beta \left(x\right)$ and $\gamma \left(x\right)$ are at least twice differentiable functions of the dependent variable only, and $C_0$ is a constant driving force. The presence of both linear and quadratic damping terms, the term $\gamma \left(x\right)x$, which can make the ODE nonlinear, as well as the driving force makes this ODE general enough to represent all population dynamics models considered in this paper (see Section 3. In a special case of $\alpha \left(x\right)=\beta \left(x\right)=C_0=0$, the equation represents a harmonic oscillator [2,3].

Based on definition of null Lagrangians [7,30,31,32], the term $\beta \left(x\right)\dot{x}$ is by itself a null Lagrangian that identically satisfies the E-L equation [21,22,23], regardless of the form of the function $\beta \left(x\right)$. This implies that derivation of the equation of motion with this term present is not possible [4,5,35,36] or is only restricted to very special Lagrangians [37,38,39]. Therefore, we follow [2,3] and account for this linear damping term by writing equation (2) as

$$\ddot{x}+\alpha \left(x\right){\dot{x}}^2+\gamma \left(x\right)x=F(x,\dot{x}),$$

(3)

where the force-like term is given by

$$F(x,\dot{x})=C_0-\beta \left(x\right)\dot{x}.$$

(4)

Note that in case, $F(x,\dot{x})=0$, the above equation reduces to that considered in [32]. In the following, the originally developed method to find SLs [32] is now generalized to the above form of the nonhomogeneous ODEs.

Let us consider the following integral transform

$$x\left(t\right)=x_1\left(t\right)e^{I_{\varphi }\left(x_1\right)},$$

(5)

where

$$I_{\varphi }\left(x_1\right)={\int }^{x_1}\varphi \left(x_2\right)d{x}_2,$$

(6)

and calculate

$$\dot{x}=(1+x_1\varphi ){\dot{x}}_1{e}^{I_{\varphi }\left(x_1\right)},$$

(7)

$${\dot{x}}^2=(1+x_1\varphi {)}^2{\dot{x}}_1^2{e}^{2I_{\varphi }\left(x_1\right)},$$

(8)

We now substitute these transformed variables into equation (3) and obtain

$$\begin{array}{cc}\hfill {\ddot{x}}_1& +\frac{1}{1+x_1\varphi }[2\varphi +x_1{\varphi }^2+x_1(\frac{d\varphi }{d{x}_1})+b(x_1{e}^{I_{\varphi }\left(x_1\right)}\left)\right(1+x_1\varphi \left){}^2{e}^{I_{\varphi }\left(x_1\right)}\right]{\dot{x}}_1^2\hfill \\ \hfill & +\frac{1}{1+x_1\varphi }\gamma \left(x_1{e}^{I_{\varphi }\left(x_1\right)}\left)x_1=\frac{e^{-I_{\varphi }\left(x_1\right)}}{1+x_1\varphi }F\right(\right(1+x_1\varphi \left){\dot{x}}_1{e}^{I_{\varphi }\left(x_1\right)},x_1{e}^{I_{\varphi }\left(x_1\right)}\right),\hfill \end{array}$$

In order to remove the term with ${\dot{x}}_1^2$, the function $\varphi$ must obey the following nonlinear ODE

$$2\varphi +x_1{\varphi }^2+x_1\left(\frac{d\varphi }{d{x}_1}\right)+b(x_1,e^{I_{\varphi }\left(x_1\right)}){(1+x_1\varphi )}^2{e}^{I_{\varphi }\left(x_1\right)}=0,$$

(9)

whose solutions $\varphi \left(x_1\right)$ were previously found [33]. The fact that the solutions to this nonlinear equation are found simply means that there is a function $\varphi \left(x_1\right)$ that allows removing the nonlinear damping term ${\dot{x_1}}^2\left(t\right)$ from the resulting equation of motion for the transformed variable $x_1\left(t\right)$. This is an important result since the standard Lagrangian for a linear ODE of this form is already known [33]. Thus, the known Lagrangian for the transformed equation of motion can be now transformed back to the original variable $x\left(t\right)$ as it is shown below.

The inverse integral transformation can be performed, and the transformed standard Lagrangian for the variable $x_1\left(t\right)$ and the Euler-Lagrange equation can now be expressed in terms of the original variable $x\left(t\right)$. This gives

$$L(\dot{x},x)=\frac{1}{2}{\dot{x}}^2{e}^{2I_{\alpha }\left(x\right)}-\underset{}{\overset{x}{\int }}\tilde{x}\gamma \left(\tilde{x}\right)e^{2I_{\alpha }\left(\tilde{x}\right)}d\tilde{x},$$

(10)

where

$$I_{\alpha }\left(x\right)=\underset{}{\overset{x}{\int }}\alpha \left(\tilde{x}\right)d\tilde{x},$$

(11)

and $\widehat{EL}\left[L(\dot{x},x)\right]=F(\dot{x},x)e^{2I_{\alpha }\left(x\right)}$ or more explicitly

$$\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{x}}\right)-\frac{\partial L}{\partial x}=F(\dot{x},x)e^{2I_{\alpha }\left(x\right)}.$$

(12)

It must be noted that the presence of the term $F(\dot{x},x)e^{2I_{\alpha }\left(x\right)}$ is justified by the fact that this term does not arise from any potential [2].

3. Applications to the population dynamics models

3.1. Selected models

The models of the population dynamics considered in this paper are listed in Table 1. Our selection process was guided by the previous work of Trubatch and Franco [22], and Nucci and Tamizhmani [28]. Both papers considered the well-known population models that involve two interacting species described by coupled nonlinear ODEs, namely, the Lotka-Volterra, Gompertz, Verhulst and Host-Parasite models as shown in Table 1. The authors of these papers determined Lagrangians corresponding to the ODEs representing mathematically the models, either ad hoc [22] or using the method of Jacobi Last Multiplier [28]; the Lagrangians of the same form were obtained, and they were treated as the generating functions for the ODEs representing the models.

As it is well-known, Lagrangians can be of different forms and yet they would give the same equation of motion [5,9,10]. In most cases, the forms of these Lagrangians do not resemble the SLs in which the kinetic and potential energy-like terms can be identified [8,9,10,11,12]. However, the main objective of this paper is to derive the SLs and for some selected models of the population dynamics and compare the obtained SLs to those previously found [22,28]. Because of the specific physical meaning of the SLs derived here, we are also able to address the role and meaning of these SLs in the population dynamics.

In selecting models of the population dynamics, we adopted the four models used in the previous studies [22,28]. In addition, we selected the SIR model (see Table 1).

The first four models of the population dynamics presented in Table 1 describe two interacting species (preys and predators) of the respective populations $w_1\left(t\right)$ and $w_2\left(t\right)$ that evolve in time t, which is denoted by the time derivatives ${\dot{w}}_1\left(t\right)$ and ${\dot{w}}_2\left(t\right)$. The coefficients a, A, b, B$f_1$, $f_2$, $m_1$ and $m_2$ are real and constant parameters that describe the interaction of the two species. The Lotka-Volterra, Verhulst and Gompertz models are symmetric, which means that the dependent variables can be replaced if, and only if, the constants are replaced, $a\to A$, $b\to B$, $f_1\to f_2$ and $m_1\to m_2$. However, the Host-Parasite model is asymmetric in the dependent variables.

The SIR model presented in Table 1 describes the spread of a disease in a population and the dependent variables $w_1\left(t\right)$ and $w_2\left(t\right)$ represent susceptible and infectious populations, with a and b being the recovery and infection rates, respectively. Similar to the Host-Parasite model, the SIR model is also asymmetric but the origin and nature of this asymmetry in both models is significantly different.

3.2. Standard Lagrangians for selected models

Our method to derive standard Lagrangians for the models presented in Table 1 requires that the systems of coupled nonlinear first-order ODEs are cast into one second-order ODE for a selected variable. Since the considered models of the population dynamics are either symmetric or asymmetric, we decided to derive the equations of motion for both variables, so we may see the effects caused by symmetry and asymmetry on the equations of motion for both variables.

All derived second-order ODEs can be expressed in the same form as the equation (3), which can be written as

$${\ddot{w}}_i+{\alpha }_i\left(w_i\right){\dot{w}}_i^2+{\gamma }_i\left(w_i\right)w_i=F_i({\dot{w}}_i,w_i),$$

(13)

where $i=1$ and 2. Since $w_i\left(t\right)$ represents the population of species, its derivative with respect of time ${\dot{w}}_i\left(t\right)$ describes the rate with which the population changes, and ${\ddot{w}}_i\left(t\right)$ its acceleration. Despite the presence of the damping-like term ${\alpha }_i\left(w_i\right){\dot{w}}_i^2$, the LHS of the above equation is conservative [9,33] and it describes oscillations of the population of species with respect to its equilibrium. These oscillations are modified by the force-like term on the RHS of the equation. Let us now describe this term.

Typically, the presence of any term with ${\dot{w}}_i\left(t\right)$ corresponds to friction forces in classical mechanics. In the approach presented in this paper, all friction-like terms that explicitly depend on ${\dot{w}}_i\left(t\right)$ are collected on the RHS of the equation as $F_i({\dot{w}}_i,w_i)$, which becomes the force-like term. Since $F_i({\dot{w}}_i,w_i)$ arises directly from the friction-like terms, its origin is not potential, and therefore this force-like term may appear on the RHS of the E-L equation (see equation 12) as it is shown in [2].

In our derivations of the standard Lagrangians for the models of the population dynamics presented in Table 1, we follow the procedure described in Section 2.3.

3.2.1. Lotka-Volterra Model

The Lotka-Volterra model was developed by Lotka [23] and Volterra [24], and this model describes the interaction of two populations (predator-prey) based on the assumption that the prey increases exponentially in time without the predator, and the predator decreases exponentially without the prey. The model is symmetric and it is represented mathematically by a system of coupled nonlinear first-order ODEs given in Table 1.

We cast the first-order ODEs are cast into the second-order ODEs of the form given by equation (3), and obtain

$${\ddot{w}}_1-\frac{1}{w_1}{\dot{w}}_1^2+(B{w}_1+A)a{w}_1=F_1({\dot{w}}_1,w_1),$$

(14a)

and

$${\ddot{w}}_2-\frac{1}{w_2}{\dot{w}}_2^2+(b{w}_2+a)A{w}_2=F_2({\dot{w}}_2,w_2),$$

(14b)

where the forcing terms are: $F_1({\dot{w}}_1,w_1)=(B{w}_1+A){\dot{w}}_1$ and $F_2({\dot{w}}_2,w_2)=(b{w}_2+a){\dot{w}}_2$. The coefficients in these equations represent: a - reproduction rate of prey, b - mortality rate of predator per prey, A - mortality rate of predator, and B - reproduction rate of predator per prey.

Our method gives the following standard Lagrangians for these equations

$$L_1({\dot{w}}_1,w_1)=\frac{1}{2}\left(\frac{{\dot{w}}_1}{w_1}\right){}^2-a(B{w}_1+Aln\left|w_1\right|),$$

(15a)

and

$$L_2({\dot{w}}_2,w_2)=\frac{1}{2}\left(\frac{{\dot{w}}_2}{w_2}\right){}^2-A(b{w}_2+aln\left|w_2\right|).$$

(15b)

The kinetic and potential energy-like terms are easy to recognize in these Lagrangians. Substituting the derived SLs and $F({\dot{w}}_i,w_i)$ into the Euler-Lagrange equations

$$\frac{d}{dt}\left(\frac{\partial L}{\partial {\dot{w}}_i}\right)-\frac{\partial L}{\partial w_i}=F({\dot{w}}_i,w_i)e^{2I_{\alpha }\left(w_i\right)},$$

(16)

where $i=1$ and 2, we obtain equations (14a) and (14b). This verifies that the presented method to derive the SLs is valid.

3.2.2. Verhulst Model

This logistical (or Verhulst) equation was first introduced by Verhulst [25]. The Verhulst model describes the organisms’ growth dynamics in a habitat of finite resources, which means the population is limited by a carrying capacity. This model is valuable for the optimization of culture media by developing strategies and the selection of cell lines. In this paper, the Verhulst model describes the population of interacting species by considering self-interacting terms that prevent the exponential increase or decrease in the size of the populations observed in the Lotka-Volterra model [22]. The system of coupled nonlinear ODEs given in Table 1 shows that the model is symmetric.

The second-order equations for the dynamical variables of this model are:

$${\ddot{w}}_1-(1+b)\frac{1}{w_1}{\dot{w}}_1^2+\left[\right(f_2-b)B{w}_1^2+(A{f}_2-2Ab-a)w_1+A(a-Ab)]w_1=F_1({\dot{w}}_1,w_1),$$

(17a)

and

$${\ddot{w}}_2-(1+B)\frac{1}{w_2}{\dot{w}}_2^2+\left[\right(f_1-B)b{w}_2^2+(a{f}_1-2aB-A)w_2+a(A-aB)]w_2=F_2({\dot{w}}_2,w_2),$$

(17b)

where the forcing terms are: $F_1({\dot{w}}_1,w_1)=-{\dot{w}}_1[(2b-1)B{w}_1-f_2{w}_1^2+(2Ab-a)+(f_2-b)B]$ and $F_2({\dot{w}}_2,w_2)=-{\dot{w}}_2[(2B-1)b{w}_2-f_1{w}_2^2+(2aB-A)+(f_1-B)b]$.

Then, the resulting standard Lagrangians are given as

$$L_1({\dot{w}}_1,w_1)=\frac{1}{2}[{\left(\frac{{\dot{w}}_1}{w_1}\right)}^2-\frac{(f_2-b)B}{(1-b)}{w}_1^2-\frac{2(A{f}_2-2Ab-a)}{(1-2b)}{w}_1+\frac{A(a-Ab)}{b}]w_1^{-2b}$$

(18a)

and

$$L_2({\dot{w}}_2,w_2)=\frac{1}{2}[{\left(\frac{{\dot{w}}_2}{w_2}\right)}^2-\frac{(f_1-B)b}{(1-B)}{w}_2^2-\frac{2(a{f}_1-2aB-A)}{(1-2B)}{w}_2+\frac{a(A-aB)}{B}]w_2^{-2B}$$

(18b)

The kinetic and potential energy-like terms are easy to identify, and the functions $F_1({\dot{w}}_1,w_1)$ and $F_2({\dot{w}}_2,w_2)$ are given above. Substitution of these Lagrangians into the E-L equations (see equation 16) validates the method.

3.2.3. Gompertz Model

Gompertz [26] proposed a model to describe the relationship between increasing death rate and age. The model is also useful in describing the rapid growth of a certain population of organisms such as the growth of tumors [27], as well as, modeling the amount of medicine in the bloodstream. Here, we follow [22,28] and consider the Gompertz model for the population dynamics. This model generalizes the Lotka-Volterra model by including self-interaction terms that prevent an unbounded increase of any isolated population [22]; the self-interacting terms in the Gompertz model are different than those in the Verhulst model. The mathematical representation of this model given by the coupled and nonlinear ODEs in Table 1 shows that the model is symmetric.

The equation describing time evolution of each model variable are given as

$${\ddot{w}}_1-\frac{1}{w_1}{\dot{w}}_1^2+[Alog(\frac{w_1}{m_1}\left)\right]w_1^2=F_1({\dot{w}}_1,w_1),$$

(19a)

and

$${\ddot{w}}_2-\frac{1}{w_2}{\dot{w}}_2^2+[alog(\frac{w_2}{m_2}\left)\right]w_2^2=F_2({\dot{w}}_2,w_2),$$

(19b)

where the forcing terms are: $F_1({\dot{w}}_1,w_1)=[A{m}_1+b{w}_1+g_1({\dot{w}}_1,w_1)]{\dot{w}}_1-g_1({\dot{w}}_1,w_1)A{w}_1$ and $F_2({\dot{w}}_2,w_2)=[a{m}_2+B{w}_2+g_2({\dot{w}}_2,w_2)]{\dot{w}}_2-g_2({\dot{w}}_2,w_2)a{w}_2$, with

$$g_1({\dot{w}}_1,w_1)=alog\left[\frac{1}{m_{2}B}\right(\frac{{\dot{w}}_1}{w_1}-Alog\left(\frac{w_1}{m_1}\right)\left)\right],$$

(20a)

and

$$g_2({\dot{w}}_2,w_2)=Alog\left[\frac{1}{m_{1}b}\right(\frac{{\dot{w}}_2}{w_2}-alog\left(\frac{w_2}{m_2}\right)\left)\right].$$

(20b)

Then, the following standard Lagrangians for this model are obtained

$$L_1({\dot{w}}_1,w_1)=\frac{1}{2}(\frac{{\dot{w}}_1}{w_1}{)}^2-A[log\left(\frac{w_1}{m_1}\right)-1]w_1,$$

(21a)

and

$$L_2({\dot{w}}_2,w_2)=\frac{1}{2}(\frac{{\dot{w}}_2}{w_2}{)}^2-a[log\left(\frac{w_2}{m_2}\right)-1]w_2.$$

(21b)

In both Lagrangians the kinetic and potential energy-like terms are seen, and the forcing functions $F_1({\dot{w}}_1,w_1)$ and $F_2({\dot{w}}_2,w_2)$ are given above. If we substitute these Lagrangians into equation (16), the second-order ODEs for the variables $w_1$ and $w_2$ are obtained.

3.2.4. Host-Parasite Model

This model describes the interaction between a host and its parasite. The model takes into account the nonlinear effects of the host population size on the growth rate of the parasite population [22]. The system of coupled nonlinear ODEs (see Table 1) is asymmetric in the dependent variables $w_1$ and $w_2$. The time evolution equations for these variables are:

$${\ddot{w}}_1-\frac{1}{w_1}(1+\frac{B}{b{w}_1}){\dot{w}}_1^2+aA{w}_1=F_1({\dot{w}}_1,w_1),$$

(22a)

and

$${\ddot{w}}_2-\frac{2}{w_2}{\dot{w}}_2^2+A(b{w}_2-a)w_2=F_2({\dot{w}}_2,w_2),$$

(22b)

where the forcing terms are: $F_1({\dot{w}}_1,w_1)=B\frac{a^2}{b}+(A-\frac{2aB}{b{w}_1}){\dot{w}}_1$ and $F_2({\dot{w}}_2,w_2)=(b{w}_2-a-A){\dot{w}}_2$.

The standard Lagrangians for the variables $w_1$ and $w_2$ are

$$L_1({\dot{w}}_1,w_1)=\frac{1}{2}(\frac{{\dot{w}}_1}{w_1}{)}^2{e}^{2B/b{w}_1}+aAEi(\frac{2B}{b{w}_1}),$$

(23a)

and

$$L_2({\dot{w}}_2,w_2)=\frac{1}{2}(\frac{{\dot{w}}_2}{w_2}{)}^2\frac{1}{w_2^2}-A[\frac{1}{2}\frac{a}{w_2}-b]\frac{1}{w_2},$$

(23b)

where the exponential integral $Ei(2B/b{w}_1)$ is a special function defined as

$$Ei\left(z\right)={\int }_{\infty }^z\frac{e^{\tilde{z}}}{\tilde{z}}d\tilde{z},$$

(24)

with $z=2B/b{w}_1$. It must be noted that $Ei\left(z\right)$ is not an elementary function. Now, the standard Lagrangian for the variable $w_2$ is given by

It is seen that there are significant differences between the Lagrangian for $w_2$ and that for $w_1$ in both the kinetic and potential energy-like terms. The differences are especially prominent in the potential energy-like terms, whose explicit dependence on the exponential integral $Ei(2B/b{w}_1)$ is a new phenomenon. The differences are caused by the asymmetry between the dependent variables in the original equations (see Table 1), which makes this model different than the fully symmetric Lotka-Volterra, Verhulst and Gompertz models, whose standard Lagrangians are also fully symmetric.

3.2.5. SIR Model

Kermack and McKendrick [29] derived the system of the first ODEs (see Table 1) describing the spread of a disease in a population. It is one of the simplest models, dividing the population into three distinct sub-populations: a susceptible population denoted by $w_1\left(t\right)$, the infectious population represented by $w_2\left(t\right)$, and a recovered population, we denote as $w_3\left(t\right)$. It is seen that the dependent variable $w_3\left(t\right)$ does not appear explicitly in the set of ODEs given in Table 1 because it is related to $w_1\left(t\right)$ and $w_2\left(t\right)$ through the following population conservation law: $d/dt(w_1+w_2+w_3)=0$, which means that the sum of the three populations must remain constant in time. Moreover, $a>0$ is the recovery rate and $b>0$ is the rate of infection, which means that the terms $-b{w}_1{w}_2$ and $-a{w}_2$ represent newly infected and recovered individuals, respectively.

The time evolution equations for the variables $w_1\left(t\right)$ and $w_2\left(t\right)$ can be written as

$${\ddot{w}}_1-\frac{1}{w_1}{\dot{w}}_1^2=F_1({\dot{w}}_1,w_1),$$

(25a)

and

$${\ddot{w}}_2-\frac{1}{w_2}{\dot{w}}_2^2+ab{w}_2^2=F_2({\dot{w}}_2,w_2),$$

(25b)

where the forcing terms are given by $F_1({\dot{w}}_1,w_1)=(b{w}_1-a){\dot{w}}_1$ and $F_2({\dot{w}}_2,w_2)=-b{w}_2{\dot{w}}_2$.

Our method gives the following standard Lagrangians

$$L_1({\dot{w}}_1,w_1)=\frac{1}{2}(\frac{{\dot{w}}_1}{w_1}{)}^2,$$

(26a)

and

$$L_2({\dot{w}}_2,w_2)=\frac{1}{2}(\frac{{\dot{w}}_2}{w_2}{)}^2-ab{w}_2.$$

(26b)

The fact that the SIR model is asymmetric is shown by the lack of the potential energy-like term in $L_1({\dot{w}}_1,w_1)$ and its presence in $L_2({\dot{w}}_2,w_2)$. However, the kinetic energy-like terms are the same for the SLs for both variables, and they are also similar to such terms in the SLs obtained for the other population dynamics models.

3.3. Discussion

The models considered in this paper can be divided into two families, namely, symmetric (Lotka-Volterra, Verhulst and Gompertz) and asymmetric (Host-Parasite and SIR) models. The SLs derived for these models are different than the Lagrangians previously obtained [22,28]. The main difference is the explicit time-dependence of those previously found Lagrangians as compared to the SLs derived in this paper.

Let us point out that the derived standard Lagrangians for the population dynamics models are the most basic as there are also alternative Lagrangians, which may have much more complicated forms and yet they give the same equations of motion [2,3,5]. Typically, the alternative Lagrangians are obtained when a null Lagrangian (see [30,32] or Section 2.1 is added to a standard Lagrangian. The functions $F_i({\dot{w}}_i,w_i)$ are null Lagrangians and therefore their contributions to the equations of motion cannot be described directly by the SLs. Instead, the E–L equation must be modified as it is shown by equation (12).

For each considered population dynamics model, its two standard Lagrangians can be written in the following form

$$L_i({\dot{w}}_i,w_i)=E_{kin,i}({\dot{w}}_i,w_i)-E_{pot,i}\left(w_i\right),$$

(27)

where $i=1$ and 2, and $E_{kin,i}({\dot{w}}_i,w_i)$, $E_{pot,i}\left(w_i\right)$ are the kinetic and potential energy terms, respectively. These terms are easy to identify in the derived SLs (see Equations 15a, 15b, 18a, 18b, 21a, 21b, 23a, 23b, 26a and 26b) and they can be used to make comparisons between the Lagrangians and models they represent.

The $E_{kin,i}({\dot{w}}_i,w_i)$ terms in all four models have the same factor ${({\dot{w}}_i/w_i)}^2/2$, where i = 1 and 2, which represents the ratio at which the population changes with respect to its value at a given time. However, for the Verhulst and Host-Parasite models, this ratio is modified by the other factors that depend on the concentration of species at a given time. It is interesting that the $E_{kin,i}({\dot{w}}_i,w_i)$ terms in the Lotka-Volterra, Gompertz and SIR models are independent from any constant parameters, but for the other two models they are; in the case of the Host-Parasite models only the variable $w_1$ shows such a dependence.

The $E_{pot,i}\left(w_i\right)$ terms of the Lotka-Volterra model depends linearly on the concentration of species; however, the Verhulst, Gompertz and Host-Parasite models also have nonlinear (second-order) terms in the concentration of species. The SIR model is exceptional as its SL for the variable $w_1$ does not depend on any potential energy-like term. On the other hand, the SL for the variable $w_2$ does depend on the potential energy-like term that is linear in this variable. In all models, the $E_{pot,i}\left(w_i\right)$ terms depend on the constant parameters that appear in the derived second-order ODEs for these models. An interesting result is the presence of logarithmic terms in the Lotka-Volterra and Gompertz models and the exponential integral $Ei$ for the variable $w_1$ for the Host-Parasite model. It must be also noted that the form of the $E_{pot,i}\left(w_i\right)$ term for the SIR model is the simplest among all the models considered here.

Now, the functions $F_i({\dot{w}}_i,w_i)$ may depend only on ${\dot{w}}_i\left(t\right)$, or on ${\dot{w}}_i\left(t\right)w_i\left(t\right)$, and the constant parameters, or may depend on higher powers of these variables, and be even the arguments of logarithmic functions. As the presented results demonstrate, the form of $F_i({\dot{w}}_i,w_i)$ significantly differs for different models, with the simplest being for the SIR and Lotka-Volterra models, and then with the increasing complexity for the Host-Parasite and Verhulst models. The most complex form of $F_i({\dot{w}}_i,w_i)$ is found for the Gompertz model.

4. New ecological and environmental implications

4.1. Conserved and nonconserved quantities

The derived standard Lagrangians can be used to gain new insights into the population dynamics models. One of the most important tasks in studying any time-evolving models is to find quantities that remain unchanged during the evolution; we refer to such quantities as conserved. Since the derived SLs do not depend explicitly on time, one may expect that the total energy $E_{tot}=E_{kin}+E_{pot}$ is conserved. To demonstrate this, we calculate the energy function $E_{fun,i}$ for the Lagrangian given by the equation (27) and obtain

$$E_{fun,i}({\dot{w}}_i,w_i)={\dot{w}}_i\frac{\partial L_i}{\partial \dot{x}}-L_i({\dot{w}}_i,w_i),$$

(28)

which gives $E_{fun,i}({\dot{w}}_i,w_i)=E_{tot,i}({\dot{w}}_i,w_i)=E_{kin,i}({\dot{w}}_i,w_i)+E_{pot,i}\left(w_i\right)$ = constant for all considered models. The energy function can be used to derive equations of motion by using [2]

$$\frac{d{E}_{fun,i}}{dt}=-\frac{\partial L_i}{\partial t},$$

(29)

which also gives $E_{fun,i}({\dot{w}}_i,w_i)=E_{tot,i}({\dot{w}}_i,w_i)$ = constant because the derived SLs do not depend explicitly on time. This makes our results significantly different from those previously obtained [22,28].

The existence of the conserved quantity $E_{tot,i}$ is easy to understand because the derived SLs describe only the homogeneous parts of the ODEs that represent the models, namely, the SLs give

$${\ddot{w}}_i+{\alpha }_i\left(w_i\right){\dot{w}}_i^2+{\gamma }_i\left(w_i\right)w_i=0,$$

(30)

which describes oscillatory systems. From an ecological point of view, this means that the populations of interacting species oscillate in time, and despite the presence of the quadratic damping $\left(w_i^2\right)$ term the system is conservative [9,33] and the resulting oscillations are periodic in time.

However, we must keep in mind that in our approach the linear damping terms are separated in the force-like functions denoted by $F_i({\dot{w}}_i,w_i)$, which significantly vary for different models. The separation was natural because these functions represent null Lagrangians [7,30,31,32], which means that no standard Lagrangian can properly account for them [11,12] because the presented Lagrangian formalism is valid only for conservative systems [4,5].

Let us now investigate the effects of $F_i({\dot{w}}_i,w_i)$ on the oscillatory parts of the models described by equation (30). The force-like function will modify this oscillatory behavior by causing the systems to reach equilibrium faster (damping) or diverge from it (driven). We begin by pointing out that $F_i({\dot{w}}_i,w_i)$ for all models, except the Gompertz model, is linear in $\dot{x}$, which means that we may write $F_i({\dot{w}}_i,w_i)=f_i\left(w_i\right)\dot{x}$, where $f_i\left(w_i\right)$ accounts for all terms that depend exclusively on $w_i$. This allows us to define the Rayleigh dissipative function [2] as

$$R_i({\dot{w}}_i,w_i)=\frac{1}{2}{f}_{e,i}\left(w_i\right){\dot{w}}_i^2,$$

(31)

where $f_{e,i}\left(w_i\right)=f_i\left(w_i\right)e^{2I_{\alpha }\left(w_i\right)}$, and write the E–L equation in the following form

$$\frac{d}{dt}\left(\frac{\partial L_i}{\partial {\dot{w}}_i}\right)-\frac{\partial L_i}{\partial w_i}=\frac{\partial R_i}{\partial {\dot{w}}_i}.$$

(32)

It must be noted that the sign of the term on the RHS of this equation is determined by the sign of the function $F_i({\dot{w}}_i,w_i)$ or $f_i\left(w_i\right)$, and that this sign is typically different for different population dynamics models (see Section 3; the minus sign means a ’damping force’ and the plus sign means a ’driving force’. For the Gompertz model, Rayleigh’s force cannot be defined because the dependence on ${\dot{w}}_i$ is nonlinear and logarithmic, which does not allow for separation of ${\dot{w}}_i$ from the remaining expressions.

For the four models for which the Rayleigh’s function can be defined, we follow [2] and find the energy function is given by

$$\frac{d{E}_{fun,i}}{dt}=-\frac{\partial L_i}{\partial t}+R_i({\dot{w}}_i,w_i),$$

(33)

or simply as

$$\frac{d{E}_{fun,i}}{dt}=R_i({\dot{w}}_i,w_i),$$

(34)

because $L_i({\dot{w}}_i,w_i)\ne L_i\left(t\right)$. Ecologically, this is an important result as integrating equation (34) in time allows finding changes of $E_{fun,i}({\dot{w}}_i,w_i)$ and their effects of damping or driving on otherwise oscillatory behavior of the population dynamics models. Obviously, in order to perform the integration, the solutions for $w_i\left(t\right)$ must be known.

4.2. Period of oscillations

The above results show that the derived standard Lagrangians give equations that describe oscillatory behavior of the considered population dynamics models. Moreover, our approach demonstrates that the oscillatory behavior and the driving or damping functions can be separated. From an ecological point of view this means that they can be studied independently from each other, and then equation (34) can be used to determine the effects of damping or driving on the model oscillations.

Our approach can also be used to find the period of oscillations the population dynamics models undergo, and how this period of oscillations would change due to the presence of the damping or driving function. We follow [22,40] and consider the minimization of the following quantity

$$Z_i\left[w_i\left(t\right)\right]={\int }_0^P[L_i({\dot{w}}_i,w_i)+E_{tot,i}({\dot{w}}_i,w_i)]dt,$$

(35)

where P is the period of the oscillations, which is an unknown, and it is required that $E_{tot,i}({\dot{w}}_i,w_i)]$ is conserved. After the integration is performed, for which knowledge of $w_i\left(t\right)$ is necessary, the value of P can be determined and it’s theoretically predicted value can be tested by comparing it to experimental data.

The presented approach and obtained standard Lagrangians allow us to gain new biological insights into population dynamics models. The following three main insights have been identified: first, the oscillatory and damping or driving functions of the models can be treated separately and studied independently; second, the period of oscillations is theoretically predicted and to be verified by comparison to experimental data; and third, our approach shows how to account for the damping or driving effects on the oscillatory behavior of the model.

Despite the fact that the SLs were derived only for the population dynamics models, the developed method to find the SLs, and our approach to using the SLs to gain insights can be easily applied to other population dynamics models as well as a broad range of different ecological and biological systems.

5. Conclusions

We developed Lagrangian formalism for the following population dynamics models: Lotka-Volterra, Verhulst, Gompertz, Host-Parasite and SIR models. For ODEs that represent these models, we solved the inverse (Helmholtz) variational calculus problem and derived standard Lagrangians for both dependent variables for each model. The main characteristic of these Lagrangians is that their kinetic and potential energy-like terms are easy to be identified and that they can be used to find conserved quantities in these models.

The derived standard Lagrangians are used to demonstrate that the oscillatory behavior described by these Lagrangians can be treated separately and studied independently from the damping or driving functions of the models, and that the period of the oscillations can be predicted theoretically and verified by experimental data. The presented approach also shows how to account for the effects caused by the presence of the damping or driving functions in the models.

Our method of solving the inverse calculus of variation problem and deriving standard Lagrangians is applied to the models of population dynamics. However, the presented results show that the method can be easily extended to other population dynamics models as well as different ecological and biological as well as other natural science systems, whose evolution equations are known.