Quantifying population dynamics via a geometric mean predator–prey model

The Lotka–Volterra (LV) predator–prey model, independently proposed in Refs. 1 and 2, and its several variants put forward afterward^3,4 have been widely used and still represent a paradigm⁵ in the nonlinear statistical analysis of population dynamics, with applications transcending the field of ecology,^6–17 as, for example, in plasma physics,¹⁸ machine learning,¹⁹ and astrophysics.²⁰ Population dynamics analysis based on the LV description is increasingly gaining importance in strategic planning and policy decision-making in relation to sustainability issues.

A Hamiltonian predator–prey model, shortly referred to as geometric mean (GM) predator–prey model, capturing the main features of the standard LV model and admitting exact solutions in terms of elementary functions has been recently proposed.²³ In the GM model, the growth and death rates are proportional to the square root of the populations $x i$ rather than the populations themselves $x i$⁠. Similar behavior is also assumed for the predating rate on prey whose food supply depends on the square root of the prey population. The interaction between predator and prey populations is described by the geometric mean of the two populations $x 1 x 2$⁠, implying that the growth rate of the predator population are $c 12 x 1 x 2$ and $− c 21 x 1 x 2$⁠, respectively.

This work is aimed at comparing the geometrical mean model to the standard Lotka–Volterra model by numerically solving the differential equations and highlighting similarities and differences in the underlying dynamics. In particular, we show that the evolution dynamics’ timescale of the LV model is amplified in the GM model, which is driven by the inter-species rather than the intra-species interactions. Hence, while the characteristic time of the LV dynamics simply depends on the predator mortality, the characteristic time of the GM model depends on the joint effect of the growth rate of the predator population and the death rate of the prey population. The difference is summarized in $r$⁠, the only parameter which appears in the rescaled GM and LV models. While in the latter $r$ depends only on the demographic parameters, in the GM model, $r$ depends also on the interaction strength between the populations. These results open new scenarios in the modeling of population dynamics to analytically describe the predator–prey time evolution when the dynamics is driven by the interaction between the two species rather than by the demographic features of each species. As stated above, the main purpose of this work is to highlight the structural differences in the mean field description of the two models hence the current analysis is kept limited to the deterministic forms of the prey–predator dynamic equations. However, one should bear in mind that real-world ecosystems are featured by stochastic rather than deterministic interactions as already suggested in the early literature dealing with oscillatory phenomena in ecological and biological networks^24,25 and still under active scrutiny. Additive and multiplicative noise sources affecting the interaction processes have been considered in Refs. 26–29, and 30 The onset of the oscillatory behavior and its offset have been related to the fluctuating interactions about the network entries accounted for in the framework of random matrix approaches.^31–36 The eigenvalue spectrum between arbitrary pairs of the random matrix elements is investigated to ascertain the situations when off diagonal or circular other than diagonal elements might impact the system stability. Overall, the effects of fluctuations, in particular, those related to the inter-species terms (off diagonal entries of the random matrices) are highly relevant to the aging and extinction of a large class of interacting systems. Hence, though not included in the present discussion, one could expect that fluctuations emerging from asymmetric terms might play an important role in the geometric mean model where the inter-species dominate over the intra-species interactions.

The paper is organized as follows. In Sec. II, the main equations describing the geometric mean (GM) predator–prey dynamics model and its scaled form are briefly recalled. In Sec. III, we perform a comparative analysis from a computational viewpoint and discuss the outcomes of the analytical solutions against numerical ones. Then, we compare the two dynamics by analyzing the respective Hamiltonians. In Sec. IV, a few remarks on the outcomes and implications of the novel predator–prey model beyond current work are addressed.

The functions $X 1(T)$ and $X 2(T)$ in Eqs. (20a) and (20b) are the analytical solutions of the GM model.²³ In Sec. III, the behavior of the geometric mean predator–prey model to the standard LV will be compared. As the LV model does not admit an analytical solution, various computational steps are required for numerically evaluate the related differential equations. These issues will be addressed in Sec. III.

The first step to numerically solve the LV differential equation system is concerned with the choice of the numerical method to be adopted among the several available approaches. We consider seven different numerical methods for the solution of ordinary differential equation (ODE) well established in the literature.³⁸ The methods are labeled by the acronyms: ODE-15-S, ODE-23, ODE-23-S, ODE-23-T, ODE-23-TB, ODE-45, and ODE-113. In the acronyms, the figures refer to the order of accuracy of the numerical integration scheme, while capital letters refer to performed upgrades in these schemes. The ODE-15-S method is a variable-order (from first to fifth) numerical solver widely employed for solving differential-algebraic equations and stiff ODEs.³⁹ The ODE-23 is a low-order integration method commonly used to solve nonstiff ODEs based on the second- and third-order Runge–Kutta (RK) technique, also known as the Bogacki–Shampine scheme. The ODE-23-S is an alternative to the classical ODE-23 approach, which is used in solving stiff ODEs based on a modified Rosenbrock triple scheme related to implicit RK methods.³⁹ The ODE-23-T and ODE-23-TB integration methods are widely employed to solve stiff ODEs based on an implicit RK formula using the trapezoidal rule for approximating numerical computing of a definite integral. In the ODE-23-TB case, the numerical calculation of the definite integral is carried out with a second-order backward differentiation formula.³⁹ The ODE-45 is a medium-order integration method widely employed to solve nonstiff ODEs based on the fourth and fifth-order RK method, also known as the Dormand–Prince scheme. Finally, the last one is a high-precision solver for solving nonstiff differential equations,³⁹ namely, the ODE-113.

To address the accuracy issues with the above integration schemes in solving the GM model evolution Eqs. (8a) and (8b), we consider 5000 time samples with $T∈[0,30]$⁠, $r=1.2$ and $X 10= X 20=2$⁠. Then, we compare the obtained numerical solutions with the analytical ones given by Eqs. (20a) and (20b). To quantify the accuracy of each numerical result, we use the normalized root mean square (NRMS) defined through $NRMS= ∑ i ( x i t r u e − x i n s o l ) 2 / ∑ i ( x i t r u e ) 2$⁠, where $x t r u e$ and $x n s o l$ are the exact analytical and numerical solutions, respectively. The range of NRMS values varies from 0 (minimum error) to $∞$ (maximum error). The NRMS and the runtime for the numerical simulations are summarized in Table I. We notice that no numerical method perfectly solves the system of ODE under analysis. In addition, the fastest method (ODE-23-TB) performs with a runtime 51.40 slower compared to the analytical solution. The most accurate numerical integration method (ODE-113) solves the novel predator–prey equations with a runtime 138.44 times slower than the exact analytical solution. The numerical experiments confirm that solving ODEs is time-consuming, making it difficult to use them to deal with real issues such as, for instance, modern big data analysis problems. Thus, simulating population dynamics by solving ODEs entirely depends on a good trade-off between runtime and tolerated numerical errors. These problems do not exist when the analytical solution is known.

In this section, we discuss the similarities and differences between the dynamics of the standard LV and the GM models. It is worth remembering that both dynamics follow a population growth and decline cycle of predators and prey, as discussed in Sec. II.*** However, the standard LV population model is exponential, while the GM model’s growth rate is quadratic. We simply eliminate the iteration term between the two species to analyze these rates. In this regard, from Eqs. (3a) and (3b), note that $d X 1 d T∝ X 1$ and $d X 2 d T∝r X 2$ give the growth rates associated with the standard model, leading to $X 1(T)= X 10exp⁡( α 1T)$ and $X 2(T)= X 20exp⁡(r α 2T)$⁠, where $∝$ denotes the proportionality and $α 1$ and $α 2$ proportionality constants. At the same time, from Eqs. (8a) and (8b), we notice that $d X 1 d T∝ X 1$ and $d X 2 d T∝r X 2$ give the growth rates associated with the novel predator–prey model, leading to $X 1(T)= ( 1 2 α 1 T + X 10 ) 2$ and $X 2(T)= ( r 2 α 2 T + X 20 ) 2$ for $X 10, X 20>0$⁠. A significant difference between the models is the impact of the predator–prey pair on each other. To interpret the behavior of the predation rate, we consider only the first terms of the two predator–prey models, i.e., $d X 1 d T= ( X 1 X 2 ) α$ and $d X 2 d T=− ( X 1 X 2 ) α$⁠, in which the cases $α=1$ and $α=1/2$ represent the predator–prey interaction term related to the standard LV and GM models, respectively. In this regime, a monotonically growing population of predators is expected while the prey population decreases until extinction. Such a regime is described so that the scaled predator populations are given by $X 1(T)= X 0− X 2(T)$⁠, with $X 0= X 10+ X 20$⁠, in both cases, while the scaled prey populations are, respectively, $X 2(T)= X 0/[1−(1− X 0/ X 20)exp⁡( X 0T)]$ for $T∈[0,∞)$ and $X 2(T)= X 0 sin 2 ( − ω T + arcsin ( X 20 / X 0 ) )$ for $T∈[0, T max]$ in standard and new dynamics, where $T max=2arcsin( X 20 / X 0)$⁠. The standard LV model predicts infinite dynamics over time, even after prey extinction, i.e., $X 2(T→∞)=0$⁠. On the other hand, predators grow following a logistic curve $X 1(T)= X 0/(1+( X 0/ X 20−1)exp⁡(− X 0T))$⁠, in which their population remains constant after a lack of food (prey extinction), i.e., $X 1(T→∞)= X 0$⁠. In contrast, the GM model predicts a cessation in predation dynamics when the prey population is extinct at time $T max$⁠. As a result of the resistance of the environment, represented in this instance by the absence of prey, the GM model suggests the extinction of predation dynamics when $X 2( T max)=0$ and $X 1( T max)= X 0$⁠.

Figure 1 shows the time evolution of the scaled predator $X 1$ (purple) and prey $X 2$ (red) populations as a function of the scaled time $T$ under the same initial conditions (⁠ $X 10=1.0$ and $X 20=1.0$⁠) for the LV and the GM model. Solid curves refer to $r=0.8$⁠, while dashed lines refer to $r=1.2$⁠. A small change in the $r$-value from $0.8$ to $1.2$ noticeably affect the predator–prey dynamics. Since the oscillation frequency of the standard LV system is unknown, its dynamics may be explored by evaluating a Jacobian matrix⁴⁰ in an equilibrium point (that is, through a linearization process using first-order partial derivatives). When the standard LV model is scaled according to Eqs. (2a) and (2b), the oscillation frequency, in first-order approximation, around the equilibrium point ensues $ω LV= c 22 / c 11$⁠. At the same time, the oscillation frequency associated with the GM model in its scaled version is constant with a period of $4π$⁠, as expected due to the factor $ω=1/2$ in the arguments of the sines and cosines functions in Eqs. (20a) and (20b). Such contrasts are noteworthy when we compare the predator–prey dynamics of the two models for different r-values for the GM model, as depicted in Fig. 2.

To quantitatively compare the fluctuations in population sizes shown in Fig. 2, we consider the Pearson correlation coefficient (⁠ $ρ$⁠) as a similarity metric. In this regard, the Pearson coefficient statistically measures the degree of the linear relationship between two time series to assess the similarity between their fluctuations over time. This parameter ranges from $−$1 to 1, in which a value close to 1 indicates a systematic association between changes in the fluctuations of the two signals, i.e., a strong positive correlation. On the other hand, a value close to $−$1 indicates that the amplitudes of the fluctuations tend to vary in opposite directions, i.e., a strong negative correlation. A value close to 0 indicates the absence of a systematic relationship between the fluctuations of the two signals, i.e., a weak correlation (or null when $ρ=0$⁠). We summarize in Table II the $ρ$-values between the colored and black curves depicted in Fig. 2. In all cases, Pearson coefficients assume values close to zero, indicating that the predator–prey models predict different fluctuations between them.

Figure 3 reports the predator–prey trajectories with initial conditions $X 10=1.0$ and $X 20=1.0$⁠, and $r=0.8$ (orange) and $r=1.2$ (blue), for (a) the standard LV model described by Eqs. (3a) and (3b) and (b) the GM model described by Eqs. (8a) and (8b). The black arrows point to the direction of the increasing scaled time. It is remarkable that the frequency of the evolution cycles is different for the two models. While in the standard LV model the frequency is a function of the parameter $r$⁠, in the GM model, it is a constant $ω=1/2$⁠. By considering $X 10=1$ and $X 20=1$⁠, we notice that the GM model predicts a time evolution with population sizes more prominent than those predicted by the LV model for $r>1$ (represented in Fig. 3 by the case $r=1.2$⁠; depicted by the blue lines) and smaller when $r<1$ (represented in Fig. 3 by the case $r=0.8$⁠; depicted by the orange lines). This feature is because the $r$-parameter governs the insertion of energy, associated with environmental aspects that promote prey population growth, into the predator–prey scaled system. Indeed, as the prey and predator populations initiate to interact, their populations increase when $r>1$⁠, since the prey scaled-population variation rate (⁠ $d X 2/dT$⁠) is positive, which induces an increase in the first term of Eqs. (3a) and (8b), leaving it to be more significant than the second term of these equations, conducting to the positivity of the predator scaled-population variation rate (⁠ $d X 1/dT$⁠). When $r<1$⁠, the prey scaled-population variation rate, at the commencement of the dynamics, is less than zero, causing a decline of the prey population, which in turn induces the negativity of the predator scaled-population variation rate; therefore, both populations decrease in this regime at the beginning of the predator–prey iteration, as depicted by the orange lines in Fig. 3.

Figure 4 presents the phase portrait of the two models to assess the effect of the initial conditions on the predator-predator dynamics. The colored lines represent the phase plane trajectories corresponding to $r=1.2$ and different values of the Hamiltonian (a) Eq. (4) for the standard LV model and (b) Eq. (9) for the GM model. The black dot represents the value of the predator and prey populations whose minimum Hamilton value is reached, i.e., when $( X 1=r, X 2=1)$ and $( X 1= r 2, X 2=1)$ in the LV model and in the GM model, respectively. In both dynamics, the phase portraits exhibit closed trajectories that represent the possible population dynamics of predator and prey species. Although common patterns are observed, note that the minimum energy points (black dot in Fig. 4) are positioned differently along the $X 1$-axis; they can be obtained in the steady state from Eqs. (3) and (8) for the LV and GM models, respectively. Such equilibrium points occur when $X 2=1$ in both dynamics, but it ensues for different values of $X 1$⁠, being $X 1=r$ in the LV case and $X 1= r 2$ in the GM case. The existence of this equilibrium point suggests that if both populations are extinct, they will remain extinct until an outside factor can change that.⁴ Also, predator and prey populations’ maximum and minimum possible values differ in both dynamics following Eq. (5) for the LV model and Eq. (12) for the GM dynamics.

Figure 6 illustrates the orbits, when the $r$-parameter ranges from 0.6 to 1.5, in the phase space for the two models [panel (a) for the LV model and (b) for the GM model] in the special case where the two nodes collapse in a unique one, i.e., $( X 10, X 2 ∗= X 20)$⁠. This node is given by $( X 10=r, X 20=1)$ and $( X 10= r 2, X 20=1)$⁠, for the LV and the GM model, respectively. For the orbits on the right-hand side of the node, corresponding to $r> X 10$ for the GM model and $r> X 10$ for the LV model, the maximum predator $X 1 max$ and prey $X 2 max$ populations are increasing functions of the $r$-parameter in both the models. The minimum prey population $X 2 min$ is a decreasing function of the $r$-parameter, while the minimum predator population does not vary with the $r$-parameter and assumes the value $X 1 min= X 10$⁠. For the orbits on the left-hand side of the node, corresponding to $r< X 10$ for the GM model and $r< X 10$ for the LV model, the maximum predator population remains fixed at the value $X 1 max= X 10$⁠, while the maximum prey population $X 2 max$ is a decreasing function of the $r$-parameter. The minimum predator $X 1 min$ and prey $X 2 min$ populations are increasing functions of the $r$-parameter. Although the GM model predicts a dynamic with larger predator–prey populations than the LV one for the same $r$-value (prey per capita growth rate), both exhibit stable limits in the phase portrait even when the initial conditions coincide with the equilibrium point (⁠ $X 2=1$ in both models, $X 1=r$ and $X 1= r 2$ for the LV and GM models, respectively), except for the case $r=1$⁠.

Figure 7 shows the contour plot under surface basemap representing the two models’ Hamiltonian from a three-dimensional perspective for $r=0.5$⁠, $r=1.0$⁠, and $r=1.5$⁠. The left column refers to the standard LV model, while the right one refers to the novel predator–prey model. The region annotated with colors closest to dark red (representative of the smallest amplitudes among the H-values) is centered on the equilibrium point $(r,1)$ for the standard LV model and $( r 2,1)$ for the GM model. From a visual inspection, the similarities between the Hamiltonians are noticeable, although they have different mathematical expressions. When $r=0.5$⁠, it is worth noting that both population models predict rapid predator growth [Figs. 7(a) and 7(b)], while in the case of $r=1.5$⁠, the prey population’s growth is drastically impacted, as depicted in Figs. 7(e) and 7(f). In other words, the $r$-parameter directly influences which species dominate the energy content of the predator–prey dynamics. When $r$ takes on the value of $1$⁠, the 3D representation of the Hamiltonians predicts higher energies in the extreme populations [Figs. 7(c) and 7(d)]. Clearly, the $r$-parameter impacts the morphology of the H-surface and then of the phase-space orbits.

A comparative analysis of the standard LV model and the GM predator–prey dynamics recently introduced in Ref. 23 has been performed. The GM model is integrable and its closed-form solution assumes an analytical expression involving standard trigonometric cyclic functions. The expression of the oscillation frequency is obtained and only depends on the interaction coefficients, i.e., $ω= c 12 c 21$⁠. These results are related to the particular forms of the demographic (natality and mortality) rates and the interaction between the predator and prey populations. For the comparison, the solution of the LV model has been obtained by adopting the standard numerical integration techniques available in the literature. It is shown that the GM model and the LV one exhibit the same main features despite being based on different dynamics.

To the purpose of better understanding the mechanism underlying the LV dynamics, it can be observed that the LV model foresees two types of growth and decrease for each population. The rates of decrease of predators and growth of preys depend linearly on the two populations involved. Conversely, the growth rate of predators and the rate of decrease of prey are non-linear functions that depend on the product of the two populations. When only the two linear rates are present in the dynamics, they induce an exponential decrease for the predators and an exponential growth for the preys in the system. In Ref. 41, a limiting case was considered with only the non-linear interaction between the two populations, and it was demonstrated that the dynamics of the system is subject to a saturation of the two populations described by the logistic equation. The predators grow and the preys decrease, asymptotically reaching their limit values. It can, therefore, be concluded that the LV dynamics turns out to be a combination of a linear and a logistic dynamics. The combination of the two integrable dynamics generates the non-integrable LV dynamics. This type of analysis can also be applied to the GM model²³ to better understand its underlying dynamics.

A second difference is that while in the LV model, the duration of the phenomenon is infinite; in the GM model, the number of predators decreases non-linearly and the dynamical process evolves within a finite time $0≤t≤2 X 10$⁠, i.e., when their number becomes zero due to their mortality. The preys instead continue to grow with a law that asymptotically becomes linear instead of exponential as in the case of the LV model.

As in the LV model, the overall dynamics of the GM model combines two partial dynamics. The first is the demographic dynamics that accounts for the mortality of predators and the birth rate of preys, while the second dynamics describes the phenomenon of predation accounting for the increase of predators and decrease of preys following the interaction between the two populations. The above partial dynamics are completely different in the two models considered here. The overall dynamics of both models result to be oscillating even if different in the two cases.

This work received financial support from the TED4LAT project within the WIDERA initiative of the Horizon Europe Programme, Grant Agreement No. 101079206.

The authors have no conflicts to disclose.

S. L. da Silva: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal). A. Carbone: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal). G. Kaniadakis: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal).

Data sharing is not applicable to this article as no new data were created or analyzed in this study.