# Novel predator-prey model admitting exact analytical solution

G. Kaniadakis<EMAIL_ADDRESS>Department of Applied Science and
Technology, Politecnico di Torino,
Corso Duca degli Abruzzi 24, 10129 Torino, Italy

###### Abstract

The Lotka-Volterra predator-prey model still represents the paradigm for the
description of the competition in population dynamics. Despite its extreme
simplicity, it does not admit an analytical solution, and for this reason,
numerical integration methods are usually adopted to apply it to various
fields of science. The aim of the present work is to investigate the existence
of new predator-prey models sharing the broad features of the standard Lotka-
Volterra model and, at the same time, offer the advantage of possessing exact
analytical solutions. To this purpose, a general Hamiltonian formalism, which
is suitable for treating a large class of predator-prey models in population
dynamics within the same framework, has been developed as a first step. The
only existing model having the property of admitting a simple exact analytical
solution, is identified within the above class of models. The solution of this
special predator-prey model is obtained explicitly, in terms of known
elementary functions, and its main properties are studied. Finally, the
generalization of this model, based on the concept of power-law competition,
as well as its extension to the case of $N$-component competition systems, are
considered.

##  I Introduction

The standard Lotka-Volterra (LV) predator-prey model is described by the
bilinear first-order coupled differential equations

$\displaystyle\frac{d\,x_{1}}{d\,t}=c_{12}\,x_{1}\,x_{2}-c_{11}\,x_{1}\ \ ,$
(1.1) $\displaystyle\frac{d\,x_{2}}{d\,t}=-c_{21}\,x_{1}\,x_{2}+c_{22}\,x_{2}\
\ ,$ (1.2)

where the positive functions $x_{1}=x_{1}(t)$ and $x_{2}=x_{2}(t)$ represent
the predator and the pray populations, respectively, while $c_{ij}$ are
positive constants describing the interaction between the two populations.

The above deterministic model was created in the third decade of the twentieth
century in the field of ecology [1, 2]. Almost a century after it was first
proposed it still attracts the interest of the scientific community, and
intense research activity has been carried out on this model up to the present
day [3, 4, 5, 6, 7, 8, 9, 10, 11, 12].

Ecology undoubtedly represents the field in which the LV model has been used
the most and in a systematic way, but its versatility has allowed it to be
used for the study of an ever increasing variety of physical, natural and
artificial systems, such as in plasma physics [13], in spin-wave patterns
[14], in the formation of crystallization fronts [15], in multimode dynamics
in optical systems [16], in neural networks [17] etc. A very large, but not
exhaustive, list of applications of the LV model can be found in Refs. [18,
19].

The undisputed success of the LV model has not prevented researchers from
going further and proposing generalizations of the model in order to describe
some empirically observed phenomenologies, especially in the field of ecology
and biology [20, 21, 22, 23, 24]. A huge debate is currently ongoing in
ecology on how to model predation, and the discussion on how the predator’s
consumption rate (known as functional response) can be influenced by the
predator and the prey population densities, has been stimulated by some works
[25, 26, 27, 28, 29]. In prey-dependent models the functional response is
independent on predator population size while in the predator-dependent models
the functional response depends of both predator and prey population sizes.
Special case of predator-dependent models are the so called ratio-depended
models (ratio of prey population size to predator population size) promoted by
Arditi and Ginzburg. A critical discussion on the principal predator-prey
models and limits of their validity can be found in ref. [26]. Mathematical
generalizations of the LV model has recently been considered in the framework
of fractional calculus [30, 31].

An old but very efficient and simple tool that is used in the construction of
ordinary differential equations, and which is suitable for describing the
evolution of populations, is undoubtedly the formalism of power-law functions.
It is worth noting that the first equation, describing the evolution of a
population, proposed by Verhulst in 1838, contains a term governing the
saturation of the population, which is a second-degree power-law of the
population function. In the last two decades of the twentieth century, the
formalism of power laws was systematically employed to propose population
evolution models described by coupled, first-order, ordinary nonlinear
differential equation systems [32]. This formulism was first introduced in
theoretical biochemistry [33, 34, 35, 36]. In ecology, an important
generalization of the LV system, involving the power-law formalism, was
proposed in 1988 by Brenig [37].

Despite its apparent simplicity, the LV model is far from being understood. It
does not admit analytical solution and, for this reason, its solution is
obtained numerically or in an approximate way by considering its linearized
version. These difficulties have not prevented interest in the LV model, which
still represents the paradigm for describing the phenomenon of predator-prey
competition. However, a common feature of all the models that generalize the
LV one is that they do not lead to analytical solutions as well.

Knowledge of the possible analytical solution, in terms of known elementary
functions of an evolution model of competing populations is important from,
and not only the practical point of view because of the possibility of
conducting quick comparison with empirical data. The existence of a simple
analytical solution is also extremely important for a model, from a
theoretical point of view, as it allows us to better understand the nature of
the competition mechanism. This ultimately gives an added value to the model
itself, as it makes it more transparent.

The aim of the present work is to propose a novel predator-prey model that
captures the main features and symmetries of the standard LV model and, at the
same time, has the additional property of admitting a simple analytical
solution, in terms of known elementary functions. Furthermore, the
generalization of this model, in the context of the power-law function
formalism, yields a second, more general, model that presents a simple
Hamiltonian form, suitable for numerical evaluation.

The paper is organized as follows: First, we introduce a general class of
Hamiltonian predator-prey models, where the competition terms are expressed
through two arbitrary functions and its main properties are studied. In
particular, we obtain the first integral of motion of this class of models and
focus on their Hamiltonian formalism. Subsequently, we identify the only
existing model, within this general class of models, admitting an exact
analytical solution and obtain it in explicit form. Then, as a first extension
of this integrable model, a class of models involving power-law competition
rates, is proposed. Finally a further extension of the formalism developed for
the here proposed predator-prey system, in the case of an arbitrary number of
interacting populations, is also studied.

##  II A general class of Predator-Prey Models

Four fundamental properties of the standard LV model immediately emerge from a
direct inspection of the two equations expressing the rate of change of the
two predator-prey populations:

i) The rate of change of each population is assumed to be the difference
between rate of growth and rate of loss.

ii) The competition term plays the growth rate role for the predator and of
loss rate role for the prey, and it is obtained by considering the product of
the contribution of the two populations.

iii) The contribution of the two populations in the expressions of the rate
terms is quantified through two characteristic functions.

iv) The two characteristic functions of the competing populations are assumed
to be linear functions of the related populations $f_{i}(x_{i})=x_{i}$.

Hereafter, we consider a class of predator-prey models for which only the
first three properties of the standard LV model continue to apply i.e. i), ii)
and iii). The fourth property is released, so that the new class of predator-
prey models is obtained, starting from the equations that define the LV model
and after substituting $x_{i}\rightarrow f_{i}(x_{i})$, in all the rate term,
in the right hand side of Eqs. (1.1) and (1.2), thus obtaining

$\displaystyle\frac{d\,x_{1}}{d\,t}=c_{12}\,f_{1}(x_{1})\,f_{2}(x_{2})-c_{11}\,f_{1}(x_{1})\
\ ,$ (2.3)
$\displaystyle\frac{d\,x_{2}}{d\,t}=-c_{21}\,f_{1}(x_{1})\,f_{2}(x_{2})+c_{22}\,f_{2}(x_{2})\
\ .$ (2.4)

The above coupled differential equations describe a class of predator-prey
models that capture many of the features of the standard LV one. This class of
models is very general, because of the arbitrariness of the two characteristic
functions $f_{i}(x_{i})$ that obey the conditions $f_{i}(0)=0$ and
$d\,f_{i}(x_{i})/d\,x_{i}>0$.

The two coupled first order differential Eqs. (2.3) and (2.4) can be easily
uncoupled obtaining the two second order differential equations for the
predator and prey populations. First we introduce the transformed populations
$w_{i}=w_{i}(t)$ defined through $w_{i}=f_{i}(x_{i})$ with $i=1,2$ and the
auxiliary functions $\Lambda_{i}(w)$ according to

$\displaystyle\frac{d\,\Lambda_{i}(w)}{d\,w}=\frac{1}{w\,f^{\prime}_{i}(f_{i}^{-1}(w))}\
\ ,$ (2.5)

$f^{\prime}_{i}(x_{i})$ and $f_{i}^{-1}(x_{i})$ indicating the derivative and
the inverse function of $f_{i}(x_{i})$ respectively. Then the introduction of
the two non-linear first order differential operators

$\displaystyle{\cal
D}_{ij}(w)=\frac{c_{ii}}{c_{ij}}+\frac{1}{c_{ij}}\,\frac{d\,\Lambda_{i}(w)}{d\,w}\,\frac{d\,w}{d\,t}\
\ ,$ (2.6)

permits us to write the predator evolution equation in the form of the
following non-linear second-order differential equation

$\displaystyle{\cal D}_{21}\\!\left({\cal
D}_{12}(w_{1})\right)+w_{1}=\lambda_{21}\ \ ,$ (2.7)

with $\lambda_{ij}=2\,c_{ii}/c_{ij}$. The prey evolution equation follows from
the predator one by exchanging the indexes $1\leftrightarrow 2$ and inverting
time $t\rightarrow-t$. It is remarkable that the dynamics of the predator-pray
system is univocally fixed by the forms of the two functions $f_{i}(x)$ or
equivalently by the forms of the two functions $\Lambda_{i}(w)$.

##  III Hamiltonian formalism

To introduce the Hamiltonian formalism of the above general class of models we
first obtain the first integral of motion of the system. By direct comparison
of Eqs. (2.3) and (2.4), the differential equation follows as

$\displaystyle\frac{d\,x_{2}}{d\,x_{1}}=-\frac{f_{2}(x_{2})}{f_{1}(x_{1})}\,\frac{c_{21}\,f_{1}(x_{1})-c_{22}}{c_{12}\,\,f_{2}(x_{2})-c_{11}}\
\ .$ (3.8)

After separation of the variables

$\left(c_{12}-\frac{c_{11}}{f_{2}(x_{2})}\right)dx_{2}+\left(c_{21}-\frac{c_{22}}{f_{1}(x_{1})}\right)dx_{1}=0\
,$ (3.9)

and integration, the first integral of motion is obtained in the form

$\displaystyle
H=c_{21}\,x_{1}-c_{22}\,F_{1}(x_{1})+c_{12}\,x_{2}-c_{11}\,F_{2}(x_{2})\ ,\ \
$ (3.10)

where $H$ is the integration constant while the $F_{i}(x_{i})$ functions are
defined according to

$\displaystyle\frac{dF_{i}(x_{i})}{dx_{i}}=\frac{1}{f_{i}(x_{i})}\ \ .$ (3.11)

The quantity $H$, can be viewed as a function of the two variables $x_{i}$
i.e. $H=H(x_{1},x_{2})$ and, after taking into account the expression of $H$
given by Eq. (3.10), the definition (3.11), the evolution equations (2.3) and
(2.4) it follows that

$\frac{d\,H}{d\,t}=\frac{\partial\,H}{\partial\,x_{1}}\,\frac{d\,x_{1}}{d\,t}+\frac{\partial\,H}{\partial\,x_{2}}\frac{d\,x_{2}}{d\,t}=0\
,$ (3.12)

so that we can conclude that $H$ is a conserved quantity whose value $H_{0}$,
depends on the initial conditions. Thus, equation $H(x_{1},x_{2})=H_{0}$
defines, in the phase space i.e. the $x_{1}x_{2}$ plane, a level curve of the
function $H(x_{1},x_{2})$ that is the orbit of the system.

By employing the above introduced functions $F_{i}(x_{i})$, the predator-prey
evolution Eqs. (2.3), (2.4) can be written in the form

$\displaystyle\frac{d\,F_{1}(x_{1})}{d\,t}=c_{12}\,f_{2}(x_{2})-c_{11}\ \ ,$
(3.13) $\displaystyle\frac{d\,F_{2}(x_{2})}{d\,t}=c_{22}-c_{21}\,f_{1}(x_{1})\
\ .$ (3.14)

The introduction of the transformation

$\displaystyle y_{i}=F_{i}(x_{i})\ ,$ (3.15)

appears natural at this point. The evolution equations for the transformed
populations $y_{i}=y_{i}(t)$ become

$\displaystyle\frac{d\,y_{1}}{d\,t}=c_{12}\,\phi_{2}(y_{2})-c_{11}\ ,$ (3.16)
$\displaystyle\frac{d\,y_{2}}{d\,t}=-c_{21}\,\phi_{1}(y_{1})+c_{22}\ ,$ (3.17)

with $\phi_{i}(y_{i})=f_{i}(\Phi_{i}(y_{i}))$ and
$\Phi_{i}(y_{i})=F_{i}^{-1}(y_{i})$ being the inverse function of
$F_{i}(x_{i})$.

As a consequence, the inverse transformations of the ones defined through Eq.
(3.15) assume the form

$\displaystyle x_{i}=\Phi_{i}(y_{i})\ ,$ (3.18)

so that the conserved quantity $H$ of the system governed by evolution Eqs.
(2.3) and (2.4), if expressed in terms of the transformed populations $y_{i}$,
takes the form

$H=c_{21}\,\Phi_{1}(y_{1})-c_{22}\,y_{1}+c_{12}\,\Phi_{2}(y_{2})-c_{11}\,y_{2}\
.$ (3.19)

An alternative expression of the function $\phi_{i}(y_{i})$ follows easily,
after taking into account Eqs. (3.11) and (3.15). It thus obtains

$\displaystyle\phi_{i}(y_{i})\\!\\!\\!\\!$
$\displaystyle=f_{i}(\Phi_{i}(y_{i}))=f_{i}(x_{i})=\frac{1}{\frac{d\,F_{i}(x_{i})}{d\,x_{i}}}=\frac{d\,x_{i}}{d\,y_{i}}\
,\ \ \ \ \ \ \ \ \ $ (3.20)

and after taking into account Eq. (3.18), it follows

$\displaystyle\phi_{i}(y_{i})=\frac{d\,\Phi_{i}(y_{i})}{d\,y_{i}}\ .$ (3.21)

The thus obtained direct differential link between the $\phi_{i}(y_{i})$ and
$\Phi_{i}(y_{i})$ functions, permits us to write evolution equations (3.16))
and (3.17), for the transformed populations $y_{i}$ in canonical form

$\displaystyle\frac{d\,y_{1}}{d\,t}=\frac{\partial\,H}{\partial\,y_{2}}\ ,$
(3.22)
$\displaystyle\frac{d\,y_{2}}{d\,t}=-\frac{\partial\,H}{\partial\,y_{1}}\ .$
(3.23)

We can then conclude that the conserved quantity $H$, is the Hamiltonian of
the transformed canonical system $\\{y_{i}(t)\\}$. In other words,
transformation (3.15) maps the real system $\\{x_{i}(t)\\}$ into the canonical
one $\\{y_{i}(t)\\}$.

It is remarkable that, according to the present formalism, any predator-prey
system described by evolution equations (2.3) and (2.4), is to associated one
and only one canonical system, whose dynamics is governed by evolution
equations (3.16) and (3.17). The nature of the predator-prey competition is
fixed by the two functions $f_{i}(x_{i})$. Furthermore, the functions
$f_{i}(x_{i})$ univocally determine the link between the real and the
canonical system by fixing the functions $F_{i}(x_{i})$ and $\Phi_{i}(x_{i})$
and also the functions $\phi_{i}(x_{i})$ which govern the competition dynamics
for the canonical system, according to the scheme

$\displaystyle f_{i}(x_{i})\Longleftrightarrow
F_{i}(x_{i})\Longleftrightarrow\Phi_{i}(y_{i})\Longleftrightarrow\phi_{i}(y_{i})\
.$ (3.24)

It is worth noting that the transformations that allows the real system to be
mapped to the canonical system, i.e. the functions $\Phi_{i}(y_{i})$ can be
obtained directly from the functions $f_{i}(x_{i})$. Indeed, Eq. (3.21), after
taking into account that $\phi_{i}(y_{i})=f_{i}(\Phi_{i}(y_{i}))$, assumes the
form

$\displaystyle\frac{d\,\Phi_{i}(y_{i})}{d\,y_{i}}=f_{i}(\Phi_{i}(y_{i}))\ .$
(3.25)

Finally, from the last differential equation and after taking into account Eq.
(3.21), it easily follows also the direct link between the functions
$f_{i}(x_{i})$ and $\phi_{i}(y_{i})$

$\displaystyle\phi_{i}(y_{i})=f_{i}\left(\,\int_{0}^{y_{i}}\\!\phi_{i}(w)\,dw\right)\
.$ (3.26)

Starting from a given predator-prey system, in which $f_{i}(x_{i})$ is fixed,
the two latter equations can be employed to obtain the associated Hamiltonian
system. For instance, in the case of the standard LV system, described by Eqs.
(1.1) and (1.2) where $f_{i}(x_{i})=x_{i}$, we obtain $\Phi_{i}(z)=\exp(z)$
and $\phi_{i}(z)=\exp(z)$, by employing Eqs. (3.25) and (3.26) respectively.
The evolution equations for the canonical counterpart of the standard LV
model, then reads

$\displaystyle\frac{d\,y_{1}}{d\,t}=c_{12}\,\exp(y_{2})-c_{11}\ ,$ (3.27)
$\displaystyle\frac{d\,y_{2}}{d\,t}=c_{22}-c_{21}\,\exp(y_{1})\ .$ (3.28)

The introduction of the Hamiltonian formalism is not of great utility for the
standard LV model, from the computational point of view. The system, also in
its canonical form, already known in literature [38], continues to not admit
any analytical solution. Interestingly, the Hamiltonian formalism turns out to
be very useful for the construction of a new, exactly solvable, predator-prey
model, as will be seen below.

##  IV The new exactly solvable model

Hereafter, we are interested in identifying the predator-prey model among the
infinity of models described by Eqs. (2.3) and (2.4), if it exists, that
admits an explicit analytical solution. This requirement can be imposed on the
Hamiltonian system so that, starting from a solvable Hamiltonian system, we
can go back to identify the corresponding real system by exploiting the
transformation that links the two systems. The existence of an analytical
solution for the Hamiltonian system described by Eqs. (3.16) and (3.17), can
be imposed by requiring it to be linear, i.e.,

$\displaystyle\phi_{i}(y_{i})=2\,y_{i}\ .$ (4.29)

The non-essential multiplicative constant has been chosen equal to $2$, in
such a way that the transformation allowing to pass from the Hamiltonian
system to the real one has the following simple form:

$\displaystyle\Phi_{i}(y_{i})=y_{i}^{2}\ ,$ (4.30)

as imposed by Eq. (3.21).

Starting from the above expression of the functions $\phi_{i}(y_{i})$ or
$\Phi_{i}(y_{i})$ and by employing Eqs. (3.26) and (3.25), respectively, the
functions $f_{(}x_{i})$ follows immediately. For instance, after substituting
the expression of $\Phi_{(}y_{i})$, as given by Eq. (4.30), into Eq. (3.25),
it obtains the functional equation $2\,y_{i}=f_{i}(y_{i}^{2})$, which yields

$\displaystyle f_{i}(x_{i})=2\,\sqrt{x_{i}}\ .$ (4.31)

From Eq. (3.11) and the expression of $f_{i}(x_{i})$, it follows that

$\displaystyle F_{i}(x_{i})=\sqrt{x_{i}}\ .$ (4.32)

At this point, after posing $\gamma_{ii}=2\,c_{ii}$, $\gamma_{ij}=4\,c_{ij}$
with $i=1,2$, $j=1,2$ and $i\neq j$, for simplicity, we can write the
evolution equations for the real system

$\displaystyle\frac{d\,x_{1}}{d\,t}=\gamma_{12}\,\sqrt{x_{1}}\,\sqrt{x_{2}}-\gamma_{11}\,\sqrt{x_{1}}\
\ ,$ (4.33)
$\displaystyle\frac{d\,x_{2}}{d\,t}=-\gamma_{21}\,\sqrt{x_{1}}\,\sqrt{x_{2}}+\gamma_{22}\,\sqrt{x_{2}}\
\ ,$ (4.34)

and for the corresponding Hamiltonian system

$\displaystyle\frac{d\,y_{1}}{d\,t}=\frac{\gamma_{12}}{2}\,y_{2}-\frac{\gamma_{11}}{2}\
\ ,$ (4.35)
$\displaystyle\frac{d\,y_{2}}{d\,t}=-\frac{\gamma_{21}}{2}\,y_{1}+\frac{\gamma_{22}}{2}\
\ .$ (4.36)

together with the quadratic transformation

$\displaystyle x_{i}=y_{i}^{2}\ ,$ (4.37)

that links the two systems.

It is remarkable that evolution equations (4.33) and (4.34) describe a new
simple but different from the standard LV predator-prey model. Interestingly,
this new model captures the main features of the LV model and additionally has
the property to admit an exact analytical solution, which can be expressed in
terms of known elementary functions. It is important to note that the
discovery of the above integrable model does not exclude the existence of
other integrable models within the class described by Eqs (2.3) and (2.4).

The conserved quantity $H$ assumes, for the two systems, the forms

$\displaystyle
H=\frac{\gamma_{21}}{4}\,x_{1}-\frac{\gamma_{22}}{2}\,\sqrt{x_{1}}+\frac{\gamma_{12}}{4}\,x_{2}-\frac{\gamma_{11}}{2}\,\sqrt{x_{2}}\
,\ \ $ (4.38)

and

$H=\frac{\gamma_{21}}{4}\,y_{1}^{2}-\frac{\gamma_{22}}{2}\,y_{1}+\frac{\gamma_{12}}{4}\,y_{2}^{2}-\frac{\gamma_{11}}{2}\,y_{2}\
.$ (4.39)

The above expressions of $H$ defines the orbits of the system in the two
$x_{1}x_{2}$ and $y_{1}y_{2}$ phase spaces. After some tedious but simple
algebra, and after posing

$y_{1c}=\frac{\gamma_{22}}{\gamma_{21}}\ \ \ ,\ \ \
y_{2c}=\frac{\gamma_{11}}{\gamma_{12}}\ ,$ (4.40)
$h=H+\frac{\gamma_{11}^{2}}{4\,\gamma_{12}}+\frac{\gamma_{22}^{2}}{4\,\gamma_{21}}\
,$ (4.41)

and

$a_{1}=\sqrt{\frac{4\,h}{\gamma_{21}}}\ \ \ ,\ \ \
a_{2}=\sqrt{\frac{4\,h}{\gamma_{12}}}\ ,$ (4.42)

Eq. (4.39), which defines the orbit of the Hamiltonian system in the
$y_{1}y_{2}$ plane, assumes the form

$\frac{(y_{1}-y_{1c})^{2}}{a_{1}^{2}}+\frac{(y_{2}-y_{2c})^{2}}{a_{2}^{2}}=1\
,$ (4.43)

thus indicating that the orbit is an ellipse.

The solution of the linear system of Eqs. (4.35) and (4.36) is given by

$\displaystyle y_{1}=y_{1c}+a_{1}\,\sin(\omega\,t+\theta_{0})\ ,$ (4.44)
$\displaystyle y_{2}=y_{2c}+a_{2}\,\cos(\omega\,t+\theta_{0})\ ,$ (4.45)

where the frequency $\omega$ is given by

$\displaystyle\omega=\frac{1}{2}\,\sqrt{\gamma_{12}\gamma_{21}}\ \ ,$ (4.46)

while $\theta_{0}$ is related to the initial conditions, according to
$\sin(\theta_{0})=(y_{10}-y_{1c})/a_{1}$ and
$\cos(\theta_{0})=(y_{20}-y_{2c})/a_{2}$.

After substitution of the expressions of $\theta_{0}$, $y_{1c}$, $y_{2c}$,
$a_{1}$ and $a_{2}$ in Eqs. (4.44) and (4.45), the time evolution of the
Hamiltonian populations $y_{i}$ can be obtained, in terms of the initial
conditions $y_{i0}$ and the model parameters $\gamma_{ij}$. Finally, from the
quadratic transformation given by Eq. (4.37), we immediately obtain the time
evolution of $x_{i}$

$\displaystyle x_{1}=\bigg{[}\\!\\!\\!\\!\\!$
$\displaystyle\frac{\gamma_{22}}{\gamma_{21}}+\left(\sqrt{x_{10}}-\frac{\gamma_{22}}{\gamma_{21}}\right)\,\cos(\omega\,t)$
(4.47)
$\displaystyle+\,\sqrt{\frac{\gamma_{12}}{\gamma_{21}}}\,\left(\sqrt{x_{20}}-\frac{\gamma_{11}}{\gamma_{12}}\right)\,\sin(\omega\,t)\bigg{]}^{2},$
$\displaystyle x_{2}=\bigg{[}\\!\\!\\!\\!\\!$
$\displaystyle\frac{\gamma_{11}}{\gamma_{12}}+\left(\sqrt{x_{20}}-\frac{\gamma_{11}}{\gamma_{12}}\right)\,\cos(\omega\,t)$
(4.48)
$\displaystyle-\,\sqrt{\frac{\gamma_{21}}{\gamma_{12}}}\,\left(\sqrt{x_{10}}-\frac{\gamma_{22}}{\gamma_{21}}\right)\,\sin(\omega\,t)\bigg{]}^{2},$

The above obtained functions, $x_{i}=x_{i}(t)$, represent the general exact
analytical solution of the predator-prey system described by evolution
equations (4.33) and (4.34). It is remarkable that the competition rates in
the present model are described through the characteristic sub-linear
functions $f_{i}(x_{i})=\sqrt{x_{i}}$ which substitute the linear functions
$f_{i}(x_{i})=x_{i}$ appearing in the standard LV model. This simple
substitution is sufficient to generate a predator-prey model admitting an
exact analytical solution.

The existence of the above constructed diffeomorphism linking the real non-
linear dynamical system (exhibiting an eliptic fixed point) and the
Hamiltonian linear system (with a center), is guaranteed by the Poincare-
Bendixson theorem. This diffeomorphism is obtained here together with the
solution of the real dynamic system and its frequency in closed form. For the
standard LV system, the frequency and period can only be estimated numerically
after evaluating a complicated integral [39]. The need to have the expression
of the period of the standard system of LV in a closed form is strongly felt
and has led towards the proposal of asymptotic analytic expressions that
provide only approximate estimations [40]. This highlights the usefulness of
having in closed form the solution and the frequency of the system described
by the model proposed here, especially in view of its applications.

##  V Power-law predator-prey model

An important predator-prey model is obtained when $f_{i}(x_{i})$ are power-law
functions i.e. $f_{i}(x_{i})\propto x_{i}^{\alpha_{i}}$, with $\alpha_{i}>0$.
In such a case, $F_{i}(x_{i})$, $\Phi_{i}(x_{i})$ and $\phi_{i}(x_{i})$ also
become power-law functions i.e.

$\displaystyle f_{i}(x_{i})=\frac{1}{1-\alpha_{i}}\,x_{i}^{\alpha_{i}}\ \ ,$
(5.49) $\displaystyle F_{i}(x_{i})=x_{i}^{1-\alpha_{i}}\ \ ,$ (5.50)
$\displaystyle\Phi_{i}(y_{i})=y_{i}^{\frac{1}{1-\alpha_{i}}}\ \ ,$ (5.51)
$\displaystyle\phi_{i}(y_{i})=\frac{1}{1-\alpha_{i}}\,y_{i}^{\frac{\alpha_{i}}{1-\alpha_{i}}}\
\ .$ (5.52)

After posing $\gamma_{ii}=c_{ii}/(1-\alpha_{i})$,
$\gamma_{ij}=c_{ij}/(1-\alpha_{i})(1-\alpha_{j})$, with $j\neq i$ and $i=1,2$,
$j=1,2$, Eqs, (4.33) and (4.34), that define the competition model, assume a
particularly simple form

$\displaystyle\frac{d\,x_{1}}{d\,t}=\gamma_{12}\,x_{1}^{\alpha_{1}}\,x_{2}^{\alpha_{2}}-\gamma_{11}\,x_{1}^{\alpha_{1}}\
,\ \ \ \ \ \ \ $ (5.53)
$\displaystyle\frac{d\,x_{2}}{d\,t}=-\gamma_{21}\,x_{1}^{\alpha_{1}}\,x_{2}^{\alpha_{2}}+\gamma_{22}\,x_{2}^{\alpha_{2}}\
,\ \ \ \ \ \ \ \ $ (5.54)

while the first integral of motion becomes

$\displaystyle H\\!\\!\\!\\!\\!\\!$
$\displaystyle=(1-\alpha_{1})(1-\alpha_{2})\gamma_{21}\,x_{1}-(1-\alpha_{2})\gamma_{22}\,x_{1}^{1-\alpha_{1}}$
$\displaystyle+\,(1-\alpha_{1})(1-\alpha_{2})\gamma_{12}\,x_{2}-(1-\alpha_{1})\gamma_{11}\,x_{2}^{1-\alpha_{2}}\
.$

The associated canonical system can be introduced by considering the
transformation $y_{i}=x_{i}^{1-\alpha_{i}}$ so that, after posing
$\beta_{i}=\alpha_{i}/(1-\alpha_{i})$ and
$\eta_{ij}=\gamma_{ij}/(1+\beta_{i})$, the evolution equations become

$\displaystyle\frac{d\,y_{1}}{d\,t}=\eta_{12}\,y_{2}^{\beta_{2}}-\eta_{11}\ ,$
(5.56)
$\displaystyle\frac{d\,y_{2}}{d\,t}=-\eta_{21}\,y_{1}^{\beta_{1}}+\eta_{22}\
,$ (5.57)

while the Hamiltonian of the system assumes the form

$H=\frac{\eta_{21}}{1+\beta_{1}}\,y_{1}^{1+\beta_{1}}-\eta_{22}\,y_{1}+\frac{\eta_{12}}{1+\beta_{2}}\,y_{2}^{1+\beta_{2}}-\eta_{11}\,y_{2}\
.$ (5.58)

The above power-law predator-prey model is described by two coupled
differential equations which do not admit any analytical solution, except for
the special previously examined case corresponding to $\alpha_{i}=1/2$. The
model is form invariant describing power-law interaction also in its
Hamiltonian version and belongs to the class of models described by Eqs (2.3)
and (2.4). As a consequence, compared to other more general and complex power-
law models considered in the literature [33, 34, 35, 36, 37], it has the
advantage of admitting a Hamiltonian formalism with its first integral easily
obtained in closed form. Last but not least, this simple model represents a
two-parameter interpolation between the standard LV model ($\alpha_{i}=1$) and
the model ($\alpha_{i}=1/2$), proposed in the previous section, which shows
the solution in closed form.

##  VI $N$-Component models

The interaction between $N$ distinct populations within the LV model is
described by means of the following $N$ coupled nonlinear differential
equation

$\displaystyle\frac{d\,x_{i}}{d\,t}=\lambda_{i}\,x_{i}+x_{i}\,\sum_{j=1}^{N}\,A_{ij}\,x_{j}\,,$
(6.59)

with $i=1,2,...,N$, while the bilinear term $A_{ij}x_{i}x_{j}$, represents the
competition between the populations $x_{i}$ and $x_{j}$. An important
generalization of the above $N$-component model was proposed in Ref. [37],
involving power-law functions, whose solution can be obtained by employing
numerical integration methods.

Hereafter, we briefly consider the extension of the previously introduced
power-law model to $N$-component systems. The evolution equations of the
system read

$\displaystyle\frac{d\,x_{i}}{d\,t}=\lambda_{i}\,x_{i}^{\alpha_{i}}+x_{i}^{\alpha_{i}}\,\sum_{j=1}^{N}\,A_{ij}\,x_{j}^{\alpha_{j}}\
,\ $ (6.60)

which can be obtained, starting from Eq. (6.59), describing the LV model, by
performing the $x_{i}\rightarrow x_{i}^{\alpha_{i}}$ substitution, on its
right-hand side.

In the following we are not interested in the Hamiltonian formalism of the
power-law competition system, but instead we focus on the transformation
$z_{i}=\,x_{i}^{1-\alpha_{i}}$, which preserves the power-law form of the
system. The evolution equations for the transformed system simplifies to

$\displaystyle(1+\beta_{i})\,\frac{d\,z_{i}}{d\,t}=\lambda_{i}+\sum_{j=1}^{N}\,A_{ij}\,z_{j}^{\beta_{j}}\
,\ \ \ $ (6.61)

with $\beta_{i}=\alpha_{i}/(1-\alpha_{i})$. The above equations have a form
suitable for numerical integration.

It is remarkable that, as in the previously considered case of two-component
systems, the $N$-component systems, for $\alpha_{i}=1/2$ which implies
$\beta_{i}=1$, also admit analytical solutions, in terms of known elementary
functions. The evolution equations of the real system become

$\displaystyle\frac{d\,x_{i}}{d\,t}=\lambda_{i}\,\sqrt{x_{i}}+\sqrt{x_{i}}\,\sum_{j=1}^{N}\,A_{ij}\,\sqrt{x_{j}}\
,\ \ $ (6.62)

while the system of the equations governing the transformed system becomes
linear i.e.

$\displaystyle\frac{d\,z_{i}}{d\,t}=\frac{1}{2}\,\lambda_{i}+\frac{1}{2}\,\sum_{j=1}^{N}\,A_{ij}\,z_{j}\,,\
\ \ $ (6.63)

so that its analytical solution is easy to obtain. Finally, the transformation
$x_{i}=\,z_{i}^{2}$ permits us to obtain the analytical and explicit solution,
of the real power-law competition model described by Eq. (6.62).

##  VII Concluding Remarks

Within the general class of Hamiltonian predator-prey models described by Eqs.
(2.3) and (2.4) an integrable model has been identified. This model described
by Eqs. (4.33) and (4.34), captures the main features of the standard LV model
with the important advantage of admitting the exact analytical solution given
by Eqs. (4.47) and (4.48).

Two extensions of this model have been considered. The first one regards its
generalization to the case of two-component systems with arbitrary power-law
competition rates, which leads to the evolution equations in the form given by
Eqs. (5.53) and (5.54). The second extension of the model regards
$N$-component systems governed by evolution equations (6.60).

Future perspectives and developments of the above models may be related to
both their application to specific fields of science where the LV model is
usually employed, as well as to their further extensions. Regarding the
possible extensions of the here considered models, one can wish and welcome
their future stochastic, fractional and diffusive generalizations.
Furthermore, by modifying some of the rate terms or by adding new rate terms
to the evolution equations of the models, it may be possible to treat
specifical empirically observed phenomenologies (Allee, Holling, etc).

## References

  * [1] Lotka, A. J. Elements of Physical Biology, (Lippincott Williams and Wilkins, Philadelphia, 1925).
  * [2] Volterra, V. Variazioni e fluttuazioni del numero d’individui in specie animali conviventi, Memoria della Reale Accademia Nazionale dei Lincei 2, 31 (1926).
  * [3] Llibrea J., Paulina MartÃnez Y., Dynamics of a family of Lotka-Volterra systems in $R^{3}$, Nonlinear Anal. 199, 111915 (2020).
  * [4] Ghasemabadi A., Doust M.H.R., Investigating the dynamics of Lotka-Volterra model with disease in the prey and predator species, Int. J. of Nonlinear Anal. Appl. 12, 633 (2021).
  * [5] Jiang J., Liang F., Wu W., Huang S., On the first Liapunov coefficient formula of 3D Lotka-Volterra equations with applications to multiplicity of limit cycles, Journal of Differential Equations 284, 183 (2021).
  * [6] Morrison R.E., Data-driven corrections of partial Lotka-Volterra models, Entropy 22, 1313 (2020).
  * [7] Palombi F., Ferriani S., Toti S., Coevolutionary dynamics of a variant of the cyclic Lotka-Volterra model with three-agent interactions, European Physical Journal B 93, 194 (2020)
  * [8] Martins A.D., Gjini E., Modeling competitive mixtures with the Lotka-Volterra framework for more complex fitness assessment between strains, Frontiers in Microbiology 11, 572487 (2020).
  * [9] Bunin G., Ecological communities with Lotka-Volterra dynamics, Physical Review E 95, 042414 (2017).
  * [10] Chen B., Dynamic behaviors of a non-selective harvesting Lotka-Volterra amensalism model incorporating partial closure for the populations, Adv. in Differ. Equ. 2018, 111 (2018).
  * [11] Lin Q., Allee effect increasing the final density of the species subject to the Allee effect in a Lotka-Volterra commensal symbiosis model, Adv. Diff. Equ. 2018. 196 (2018).
  * [12] Skvortsov A., Ristic B., Kamenev A., Predicting population extinction from early observations of the Lotka-Volterra system, Applied Mathematics and Computation 320, 371 (2018).
  * [13] Laval G. and Pellat R., Plasma Physics, Proceedings of the Summer School of Theoretical Physics, (Gordon and Breach, New York, 1975).
  * [14] Elmer F.-J., Phys. Rev. Lett. 70, 2028 (1993).
  * [15] Kurtze D. A., Nonlinear development of morphological instabilities in explosive crystallizatio, Phys. Rev. B 40 11104, 104 (1989).
  * [16] Benkert C. and Anderson D. Z., Phys. Rev. A 44, 4633 (1991).
  * [17] Noonburg V.W., A neural network modeled by an adaptive Lotka-Volterra system, SIAM J. Appl. Math. 49, 1779-1792 (1989).
  * [18] Cairo L. and Feix M. R., Families of invariants of the motion for the Lotka-Volterra equations: the linear polynomials family, J. Math. Phys. 33, 2440-2455 (1992).
  * [19] Peschel M. and Mende W., The Predator-Prey Model. Do we live in a Volterra World?, (Springer-Verlag, Wien-New York, 1986).
  * [20] Rosenzweig M.L., Paradox of enrichment: Destabilization of exploitation ecosystems in ecological time, Science 171, 385 (1971).
  * [21] Albrecht F., Gatzke H., Wax N., The limit cycles in predator-prey populations, Science 181, 1073 (1973).
  * [22] Freedman H.I., Deterministic Mathematical Models in population Ecology, (Marcel Dekker, New York, p. 254 1980), p. 254.
  * [23] Arditi R., Ginzburg R., Coupling in predator-prey dynamics: ratio-dependence. J. Theor. Biol. 139 311 (1989).
  * [24] Boudjellaba H., Sari T., Dynamic transcritical bifurcations in a class of slow-fast predator-prey models, J. Diff. Eq. 246, 2205 (2009).
  * [25] Song C., Saavedra S. Bridging parametric and nonparametric measures of species interactions unveils new insights of non-equilibrium dynamics, OIKOS 130 1027 2021.
  * [26] Abrams P.A., Ginzburg L.R. The nature of predation: prey dependent, ratio dependent or neither?, Trends Ecol. Evol. 15, 337 (2000).
  * [27] Hatton I.A., MacCann K.S., Fryxall J.M., Davies J., Smerlak M., Sinclair A.R.E., Loreau M. The predator-prey power law: Biomass scaling across terrestrial and aquatic biomes, Science 349 6252, aac6284 (2015), DOI: 10.1126/science.aac6284.
  * [28] O’Dwyer J.P. Whence Lotka-Volterra?, Theor. Ecol. 11, 441-452 (2018).
  * [29] Berryman A.A. The Orgins and evolution of predator-prey theory, Ecology 73, 1530 (1992).
  * [30] Ahmed, E., El-Sayed, A.M.A., El-Saka, H.A.A., Equilibrium points, stability and numerical solutions of fractional order predator-prey and rabies models, J. Math. Anal. Appl. 325, 542 (2007).
  * [31] Das S., Gupta P.K., A mathematical model on fractional Lotka-Volterra equations Journal of Theoretical Biology 277, 1 (2011).
  * [32] Hernandez-Bermejo B., Fairen V., Lotka-Volterra Representation of General Nonlinear Systems, Mathematical Biosciences 140, 1 (1997).
  * [33] Savageau M.A., Voit E.O., Recasting nonlinear differential equations as S-systems: a canonical nonlinear form, Mathematical Biosciences 87, 83 (1987).
  * [34] Savageau M.A., Introduction to S-systems and the underlying power-law formalism, Math1 Comput. Modelling 11, 546551 (1988).
  * [35] Voit E.O., Canonical Nonlinear Modelling: S-System Approach to Understanding Complexity, Van Nostrand Reinhold, New York, (1991).
  * [36] Voit E.O. and Savageau M.A., Equivalence between S-systems and Volterra systems, Math. Biosci. 78, 47 (1986).
  * [37] Brenig L., Complete Factorisation and analytic solutions og generalized Lotka-Volterra equations, Phys. Lett. A 133, 378 (1988).
  * [38] Verhulst F., Nonlinear Differential Equations and dynamical Systems, (Springer Science and Business Media, 2006).
  * [39] Waldvogel, J. The period in the Volterra-Lotka Predator-Prey Model, SIAM J. Numer. Anal. 20, 1264 (1983).
  * [40] Oshime, Y., Asymptotic expression of the period of the Lotka-Volterra system Japan J. of Indust. and Appl. Math. 20 353-378 (2003).