Numerical approach to an age-structured Lotka-Volterra model

1.   Introduction

Continuously structured population models describe the dynamics of animal or human populations that vary continuously with respect to some underlying physiological variable, such as body size or age; see ^[1,2,3,4] and ^[5]. Interestingly, Hopf bifurcations are sudden changes in the dynamics that may occur in ecological models, leading to the emergence of complex oscillatory behavior in the population dynamics. Specifically, a Hopf bifurcation occurs when a stable population steady-state becomes unstable, giving rise to a stable limit cycle (isolated periodic orbit) in the dynamics; see, e.g., ^[6,7,8,9] and also ^{[10,11,12,13]}.

The Lotka-Volterra model, originally developed in the early 20th century, is a classic mathematical model for the interaction between predator and prey populations; see the books ^[4,14,15,16]. In the age or size structured version of this model, age or size is explicitly incorporated as a key continuous variable, allowing for richer predator-prey dynamics. The model takes into account the fact that predators preferentially target certain ages of the preys, while preys may have different reproductive successes and survival probabilities at different ages; see the papers ^[17,18,19]. We approach the predator-prey problem using two different formulations: the partial differential equation (PDE) approach and the renewal equation approach, ^[1,4,20]. Renewal equations have been used to study a wide range of phenomena in population dynamics, from the spread of infectious diseases to the dynamics of predator-prey relationships. We will analyze a non-linear renewal equation for the birth-rate of the prey population, describing how the current number of newborn preys per unit of time changes over time, based on past newborn preys and the number of predators.

On the other hand, efficient numerical methods are computational techniques designed to solve mathematical problems quickly and accurately; see the book ^[21], and the recent papers ^{[22,23,24,25]} and also ^[26,27,28]. We checked the accuracy of the methods by using available analytical solutions, and we used them to uncover the behavior of the interacting populations for the cases where analytical results are not possible.

In the present paper, we consider structured population dynamics for two interacting species where age or size of the individuals has an important impact on the dynamics. If we consider size $x\geq x_0$ (e.g., body length) as structuring variable instead of (chronological) age $a\geq 0$, under the assumptions that individuals are born at the same size $x_0$ and the individual growth rate $g$ is a function of size, that is, $\frac{dx}{dt} = g(x)$, $x(0) = x_0$, then the size-dependent problem can be reduced to an age-dependent problem. The latter is achieved by making the change of variables $a = \int_{x_0}^x \frac{ds}{g(s)}$, as it is well known. So, we deal with an age-structured model but keep in mind that it represents a reducible size-structured model, which is a particular case of physiologically structured populations; see ^[1,3,22]. Specifically, we consider the following structured predator-prey model (first-order hyperbolic system of nonlocal partial differential equations with nonlocal boundary conditions):

where the first equation is for the dynamics of the prey age-density with mortality $\mu$ and species interaction $\gamma$, the second one is for the dynamics of the predator age-density with mortality $m$ and predator efficiency $\alpha$, and the third one gives the newborns per unit of time of preys and predators with fertility rates $\beta$ and $b$, respectively. Integrating the second equation over the age span, we can reduce the system and deal with the predator population size instead of its age-density.^* In this model, one population (the predator/consumer) feeds on the other population (the prey/resource), and the interactions between the two populations can have important effects on their respective dynamics. Finally, let us mention that the spatial distribution of the interacting populations is neglected.

^* With an abuse of notation, in the following, $v$ will stand for the population size of predators instead of its age-density.

2.   PDE predator-prey model

When dealing with interacting species, age or size of the individuals of the populations involved, has been identified as one of the key features for the joint evolution; see ^[17,18,19]. Let us consider the population dynamics for the interaction of predators and preys with age structure in the prey population. As explained in the introduction, we use age instead of size for simplicity, and we are implicitly assuming that the age structure of the predator population has a small impact on the dynamics. So, let $u(a, t)$ be the density of the prey population with respect to the age $a\geq 0$ at time $t\geq 0$, and let $v(t)$ be the total population of predators at time $t\geq 0$. The structured predator-prey model of Lotka-Volterra type in $X = L^1(\mathbb{R}_+) \times \mathbb{R}$ reads as:

with suitable non-negative initial conditions $u_0(a), v_0$ at time $t = 0$. Model parameters are: predator per capita mortality rate $\delta > 0$ (actually, it is a balance between fertility and mortality in the absence of preys), ingestion coefficient $0 < \alpha < 1$ (a measure of the efficiency of the predator) and age-specific per capita rates $\mu, \beta, \gamma \in L_+^\infty(\mathbb{R}_+)$ as prey mortality [1/time], prey fertility [1/time] and species interaction [1/(time $\cdot$ ind)], respectively. Moreover, we will make the following specific assumptions:

● The age-specific mortality rate is bounded from below: $\mu(a) \geq \mu_0 > 0$, i.e., there is always a minimum mortality. In the absence of interaction between species, survival probabilities for preys and predators are $e^{-\int_0^a \mu(s)\, ds}$ and $e^{- m a}$, for some $m > \delta$, respectively.

● The age-specific interaction rate is such that $\gamma(a) > 0$ in a neighborhood of $a = 0$, i.e., the interaction between species always takes place for the younger preys.

● We take the basic reproduction number (see, e.g., ^[29,30]) for the prey population, in the absence of predators, such that:

meaning that the prey population by itself has asynchronous exponential growth. In other words, the equation for the prey population in the absence of the predators, $\partial_t u(a, t) + \partial_a u(a, t) = - \mu(a) u(a, t)$, $u(0, t) = \int_0^\infty \beta(a) u(a, t)\, da$, has asynchronous exponential growth with positive Malthusian parameter $r > 0$; see ^[1,3,5,31].

● The system (2.1) can be also considered in finite age-span $[0, a^\dagger]$ with suitable assumptions on the mortality rate (unbounded case) in order to avoid (or minimize the number of) immortal individuals.

As is well known, the original Lotka-Volterra model was an ODE model for the total population of preys and predators, providing an explanation for the periodic dynamics of fish populations observed in the Adriatic sea; see the books ^[15,16]. The system (2.1) along with the assumptions above, can be cast into an abstract semilinear differential equation in $X = L^1(\mathbb{R}_+) \times \mathbb{R}$, where the linear part is the infinitesimal generator of a strongly continuous semigroup of positive linear operators. We have that there exists a positive mild solution to (2.1) which is continuous in time; see ^[18,19,32] for the details. For the theory of semigroups applied to population dynamics; see ^[33].

Previous works on system (2.1) have focused on the PDE formulation. However, here we will focus on the so-called renewal formulation, which has some advantages, and we will develop efficient numerical methods when analytical results are not available.

The novelty of this paper lies in the following: Regarding the age-structured Lotka-Volterra model, we introduce a pseudospectral numerical method to compute the (un)stability of the coexistence steady-state and a predictor-corrector numerical scheme to simulate the population dynamics of predators and preys. In addition, we have uncovered exact solutions of the steady-state and the long-term asymptotic population dynamics to check the accuracy of the aforementioned computational methods, when either the age-dependent parameters are proportional or species interaction is age-independent. For the latter, we determine the asymptotic age-profile of the preys, which depends on the Malthusian parameter. Finally, for a biologically meaningful age-dependent interaction, the rightmost complex-conjugate spectral values cross the imaginary axis with positive speed, showing a Hopf bifurcation of the steady-state of the predator-prey system.

3.   Non-linear renewal equation

In this section, we are going to introduce a new formulation for the present predator-prey model, which simplifies the mathematical analysis. First of all, the system (2.1) can be written equivalently in integrated form as follows. Using the integration along characteristic lines in the partial differential equation (PDE) and the variation of the constants formula in the ordinary differential equation (ODE), we get to

with $u(0, t-a) = \int_0^\infty \beta(s)\, u(s, t-a)\, ds$, $a < t$. Moreover, using (2.1) and (3.1) we can get a non-linear renewal equation. Indeed, defining the new variable $b(t): = u(0, t)$ as the birth rate of the prey population (number of newborns per unit of time), we get to the renewal formulation of the present predator-prey system:

where instead of assuming a known initial condition at time $t = 0$, we assume a known history in $(-\infty, 0]$ for the birth rate of the prey population, i.e., $b(\tau) = \phi(\tau) \, , \; \tau \leq 0$, and for the predator population $v(\tau) = \psi(\tau) \, , \; \tau \leq 0$. The first history is a non-negative locally integrable function, and the second one is a non-negative continuous function.

Interestingly, we can take advantage of this renewal equation for $(b(t), v(t))$, $t\geq 0$, when computing the steady-states and their linear stability— see Sections 3.1 and 5. Moreover, non-linear renewal equations like (3.3) can be readily extended to more general structured models of Lotka-Volterra type where the prey is structured by size (rather than age), and the individual growth rate and vital rates are density-dependent. See, e.g., ^[20] for size-structured consumer resource models.

3.1.   Steady states

Under the current assumptions given in Section 2, in particular, basic reproduction number (2.2) $R_0 > 1$, there always exists both the trivial equilibrium of system (3.3), $b^* = 0$, $v^* = 0$, and the coexistence equilibrium $b^* > 0$, $v^* > 0$. Indeed, steady-states of (3.3) are given by:

Therefore, the predator population at equilibrium $v^* > 0$ is the unique (real) solution of the non-linear equation:

Notice that the equation above is independent of parameters $\alpha$ and $\delta$ and involves the decreasing function $F(a): = -\int_0^a \gamma(s) \, ds$, $a\geq 0$. If $\gamma \in L^1(\mathbb{R}_+)$ then we can consider its mean value $\frac{F(0)+F(\infty)}{2} = -\frac{1}{2} \int_0^\infty \gamma(s) \, ds$, to get an approximation of the solution (3.5) as:

Once we have $v^* > 0$ from (3.5), the birth rate of the prey population at equilibrium is computed from the second equation in (3.4) as

Moreover, the age density of the prey population at equilibrium is readily recovered from (3.1) as $u^*(a) = b^* e^{-\int_{0}^a \mu(s) + \gamma(s) v^* \, ds}$, $a\geq 0$, using (3.5)–(3.6). For an analogous computation of the steady-states using the PDE formulation and relaxing the assumption of $R_0 > 1$, see ^[19].

For numerical purposes, see Sections 5 and 6, we can consider a non-trivial case with explicit solution to the coexistence equilibrium of (3.4). Namely, the age-dependent proportional case, see, e.g., ^[24]: generic age-specific mortality $\mu(a)$, fertility rate as $\beta(a) = \bar{\beta}\, \mu (a)$ with $\bar{\beta} > 1$ and interaction rate as $\gamma(a) = \bar{\gamma}\, \mu (a)$, $\bar{\gamma} > 0$. In this case we have that the basic reproduction number (2.2) is $R_0 = \bar{\beta} > 1$ and the explicit solution to (3.5)–(3.6) is $v^* = \frac{\bar{\beta}-1}{\bar{\gamma}} > 0$ and $b^* = \frac{\delta \, \bar{\beta}}{\alpha \, \bar{\gamma}} > 0$, respectively. In this particular case, the prey age-density at equilibrium is given by $u^*(a) = b^* e^{-\bar{\beta}\int_{0}^a \mu(s) \, ds}$; see Figure 1 for an illustration with an age-specific mortality rate.

Figure 1.  Numerical simulations of the age-structured Lotka-Volterra model— see Section 6. Converging trajectories from 3 different initial conditions ($\bigtriangledown$) to the coexistence equilibrium $b^* = \frac{\delta \, \bar{\beta}}{\alpha \, \bar{\gamma}} = 40$, $v^* = \frac{\bar{\beta}-1}{\bar{\gamma}} = 3$, for the age-dependent proportional case. Parameter values: mortality $\mu(a) = \mu_0(1+k\, a)$, $\mu_0 = 0.025, k = 18$, fertility $\beta(a) = \bar{\beta}\, \mu (a)$, $\bar{\beta} = 2$, interaction $\gamma(a) = \bar{\gamma}\, \mu (a)$, $\bar{\gamma} = 1/3$ and $\delta = 5, \alpha = 0.75$. Probability of immortal individuals $\sim 10^{-10}$. Prey population size $P^* = b^* \int_0^\infty e^{-\bar{\beta}\int_{0}^a \mu(s) \, ds} da = 50.7$.

DownLoad: Full-Size Img PowerPoint

4.   Asymptotic age-profile and limit system

In this section we study the asymptotic behavior of both formulations, the PDE (2.1) and the renewal equation (3.3), through the analysis of the age-profile (i.e., normalized age-density) of the prey population.

Let the total prey population be denoted by $P(t) = \int_0^\infty u(a, t)\, da$ and let us assume for a while that the prey age-profile $\frac{u(a, t)}{P(t)}$ reaches a stationary age-profile $\bar{\omega}(a)$, $a \geq 0$, with $\int_0^\infty \bar{\omega}(a), da = 1$. In symbols,

In other words, we are assuming that there exists the asymptotic age-profile for the prey population, which is independent of the initial populations.

Next, we can integrate over the age span in (2.1), $\frac{dP}{dt} = \int_0^\infty \left[\beta(a) -\mu(a) - \gamma(a) v(t)\right] \frac{u(a, t)}{P(t)} da \cdot P(t)$, and under the assumption above, we get to the following 2-dimensional Lotka-Volterra limit system:

where the brackets above $\langle \cdot \rangle_{\bar{\omega}}$ stand for a weighted mean over the age span:

This limit system is the classical predator-prey ODE Lotka-Volterra model for the sizes of the interacting populations, so the solutions $(P(t), v(t))$, $t\geq 0$, describe closed periodic orbits around the coexistence equilibrium

with period and amplitude of oscillation depending on each initial condition $(P_0, v_0)$ of (4.2). Taking initial conditions such that $u_0(a)/P_0 = \bar{\omega}(a)$, solutions of (4.2) are implicitly given, as expected, by $\left(\frac{v}{v_0} \right)^{v^*} \!\!\! \cdot \left(\frac{P}{P_0} \right)^{\alpha P^*} = \exp [{v-v_0+\alpha(P-P_0)}]$ with period of oscillation starting at $\frac{2\pi }{\sqrt{\delta \langle\beta - \mu \rangle_{\bar{\omega}}}} = \frac{2\pi }{\sqrt{\delta v^* \langle \gamma \rangle_{\bar{\omega}} }}$, i.e., the linearization of (4.2) around the equilibrium (4.3) has pure imaginary eigenvalues $\lambda = \pm \sqrt{\delta v^* \langle \gamma \rangle_{\bar{\omega}} } \, i$. See, for instance, ^[14,15,16]. So, the asymptotic behavior of (2.1) can be determined by the existence of the asymptotic age-profile $\bar{\omega}(a)$ and the dynamics of the limit system (4.2).

Moreover, from (4.1) we have that

So, $b(t) \sim \bar{\omega}(0) P(t)$ for large $t > 0$, and the asymptotic behavior of the prey birth-rate $b(t)$ and the predator population $v(t)$ is implicitly given by

with compatible initial conditions $b_0 = \bar{\omega}(0) P_0$. The prey birth-rate at equilibrium is given by

So, again, the asymptotic behavior of (3.3) can be determined by the existence of the asymptotic age-profile $\bar{\omega}(a)$ and equation (4.4).

Finally, we need to find conditions to assure the existence of the asymptotic age-profile. Taking the approach of section 2.6 in ^[1], we can write the system for the prey age-profile $\omega(a, t): = \frac{u(a, t)}{P(t)}$, $\int_0^\infty \omega(a, t)\, da = 1$, $t\geq 0$, which is in general coupled with the predator and prey population sizes. Indeed, from (2.1) we get:

with compatible and non-trivial (in the sense of ^[1], Chapter 2) initial conditions at time $t = 0$ as $\omega_0(a) = u_0(a)/P_0$, $P_0 = \int_0^\infty u_0(a)\, da > 0$ and $v_0 > 0$, where $u_0(a)$ is the initial condition of (2.1). The situation where the support of the fertility rate $\beta(a)$ lies to the left of the support of the initial age-density $u_0(a)$ is excluded in here, as in ^[1].

It is worth to mention that $\lambda(t)$ defined in (4.6) is seen as a time-dependent Malthusian parameter for the preys but depending on the predator population. However, we realize from the first equation in (4.6) that the prey age-profile is independent of the predator population if and only if the interaction rate is age-independent, $\gamma(a) \equiv \bar{\gamma}$, $a > 0$. Indeed, the interaction term is

Therefore, when interaction between preys and predators is age-independent, system (4.6) is uncoupled, and its asymptotic behavior is completely determined as follows:

Theorem 1 (Asymptotic behavior). Under the specific assumptions given in Section 2, if in addition $\gamma(a) \equiv \bar{\gamma}$, then the solution of (4.6) with non-trivial initial conditions (in the sense of Chapter 2 in ^[1]) is such that

where $r > 0$ is the unique real solution to the non-linear equation $1 = \int_0^\infty \beta(a) \, e^{-\int_0^a \mu(s) \, ds- r \, a}\, da < R_0$. Moreover, for large $t > 0$ and a suitable constant $c > 0$, $P(t)$ and $v(t)$ in (4.6) are implicitly given by the equation

describing closed periodic orbits around the steady-state $P^* = \frac{\delta}{\alpha \bar{\gamma}}$, $v^* = \frac{r}{\bar{\gamma}}$.

Proof. For $\gamma(a) \equiv \bar{\gamma}$, the equations in (4.6) for the age-profile are uncoupled from the others:

Since we are assuming a non-trivial initial condition, we can apply Theorem 2.10 in ^[1] to assure that the asymptotic age-profile exists, and it is given by $\bar{\omega}(a) = \frac{\exp[{-\int_{0}^a \mu(s) \, ds - r\, a}]}{\int_0^\infty \exp[{-\int_{0}^a \mu(s) \, ds - r \, a}]\, da }$, where $r$ is the unique real solution to $1 = \int_0^\infty \beta(a) \, e^{-\int_0^a \mu(s) \, ds- r \, a}\, da$. Moreover, $r > 0$ since $R_0 > 1$. So, the first part of the statement follows.

On the other hand, the equations in (4.6) for $P(t)$ and $v(t)$ are influenced by the age-profile:

Once we know the evolution of the age-profile, it is straightforward to see that

Finally, if we take initial condition $\omega_0(a) = \bar{\omega}(a)$ then (4.7) coincides with (4.2); otherwise (4.2) plays the role of limit system. Accordingly, the asymptotic behavior of $(P(t), v(t))$ in (4.6) is given by the solutions of (4.2) with $\langle \beta-\mu \rangle_{\bar{\omega}} = r$ and $\langle \gamma \rangle_{\bar{\omega}} = \bar{\gamma}$. As is well known, these are periodic orbits around the steady-state $(P^*, v^*) = \left(\frac{\delta}{\alpha \bar{\gamma}}, \frac{r}{\bar{\gamma}}\right)$, implicitly given by $v^{v^*} \cdot P^{\alpha P^*} = c_0 \, e^{v+\alpha P}$, $c_0 > 0$, or equivalently $v^{r} \cdot P^\delta = c_0^{\bar{\gamma}} \, e^{\bar{\gamma}(v+\alpha P)}$, and the second part of the statement follows. □

Notice that $r > 0$ in the theorem above is the Malthusian parameter (exponential growth rate) of the preys in the absence of predators, and $\bar{\omega}(a)$ is their asymptotic age-profile. Let us recall the well-known relationship with the basic reproduction number (2.2), $\text{sign}(r) = \text{sign}(R_0-1)$. For the computation of the basic reproduction number in structured models from the analytical and numerical points of view; see, e.g., ^{[23,24,29,30]}.

Under the assumptions of Theorem 1, we can recover the asymptotic behavior of both systems (2.1) and (3.3). Namely, for large $t > 0$, the prey age-density is $u(a, t) \sim P(t)\, \bar{\omega}(a)$, $a \geq 0$, and the prey birth-rate is $b(t)\sim P(t)\, \bar{\omega}(0)$, with $P(t)$ and $v(t)$ population sizes of preys and predators oscillating periodically in time. See Figure 2 for an illustration of the closed curves (4.4) oscillating around the steady-state for the constant interaction case. One can obtain similar pictures for the periodic orbits by changing $\mu(a)$ and $\beta(a)$ while keeping the interaction $\gamma$ as age-independent.

Figure 2.  Asymptotic predator-prey dynamics for the case of age-independent interaction. Prey population has asymptotic age-profile $\bar{\omega}(a) = \bar{\omega}(0) \exp[{-\int_{0}^a \mu(s) \, ds - r\, a}]$, with $\bar{\omega}(0) = \Big(\int_0^\infty \exp[{-\int_{0}^a \mu(s) \, ds - r \, a}]\, da\Big)^{-1}$ and Malthusian parameter $r > 0$; see Theorem 1. Periodic orbits around the steady-state $b^* = \frac{\delta\, \bar{\omega}(0)}{\alpha \bar{\gamma}} = 39.7$, $v^* = \frac{r}{\bar{\gamma}} = 2.98$, given by $\left(\frac{v}{v_0} \right)^{v^*} \left(\frac{b}{b_0} \right)^{\alpha P^*} = e^{v-v_0+\frac{\alpha}{\bar{\omega}(0)}(b-b_0)}$ with $P^* = \frac{\delta}{\alpha \bar{\gamma}} = 46.7$, for three different initial conditions ($\bigtriangledown$) with $\omega_0(a) = \bar{\omega}(a)$. Parameter values: constant interaction $\gamma(a)\equiv \bar{\gamma} = 1/7$, mortality $\mu(a) = \mu_0(1+k\, a)$, $\mu_0 = 0.025, k = 18$, fertility $\beta(a) = \bar{\beta}\, \mu (a)$, $\bar{\beta} = 2$ and $\delta = 5, \alpha = 0.75$. Basic reproduction number $R_0 = \bar{\beta} = 2$ and Malthusian parameter $r = 0.43$.

DownLoad: Full-Size Img PowerPoint

On the contrary, when interaction is age-dependent actually, $\gamma(a) \not \equiv$ ct., the unique stationary age-profile is then the one corresponding to the coexistence equilibrium computed in Section 3.1:

with $v^* > 0$ the solution of (3.5). The latter is a consequence of the fact that the equation for the age-profile is not uncoupled from the equations for the population sizes.

5.   Stability analysis via pseudospectral methods

In this section we study the linear stability of the steady states using the renewal formulation (3.3). In particular, an efficient numerical method is introduced to tackle the (un)stability of the coexistence equilibrium for a general age-dependent interaction rate $\gamma(a)$. This method is based on the finite-dimensional approximation of a differential operator; see the book ^[21], and the recent papers ^{[22,23,24,25]} and also ^[26,27,28].

First of all, let us point out that the trivial steady state (extinction equilibrium) is unstable. More precisely, under the assumption of $R_0 > 1$, $(0, 0)$ is always a saddle point in the predator-prey plane, since prey population has asynchronous exponential growth in absence of predators, $b(t) \sim b_0 e^{r t}$, $r > 0$, and predator population goes to extinction in absence of preys, $v(t) = v_0 e^{-\delta t}$, $\delta > 0$.

Regarding the coexistence steady-state $b^* > 0$, $v^* > 0$ given by (3.5)–(3.6) in Section 3.1, we can linearize (3.3) to study the behavior of solutions around the equilibrium, especially for the age-dependent interaction case.

Defining new variables $B(t) = b(t)-b^*$, $V(t) = v(t)-v^*$, from (3.3) we get to the following linear renewal equation:

where, for convenience, new age-specific functions are defined

with $\int_0^\infty \beta_*(a) \, da = 1$ and $\int_0^\infty \gamma_*(a) \, da = \frac{\delta}{b^*}$ using (3.5)–(3.6).

Next, putting $B(t) = B_0 e^{\lambda t}$, $V(t) = V_0 e^{\lambda t}$ in (5.1) we get to the following eigenvalue problem, $\lambda\in \mathbb{C}$ and $(B_0, V_0)\in \mathbb{C}^2 \setminus \{(0, 0)\}$:

For an analogous computation using the PDE formulation, see ^[19]. From this eigenvalue problem, one can readily see that $\lambda = 0$ is never an eigenvalue since $\int_0^\infty \beta_*(a) \, da = 1$.

Beyond the age-independent interaction case, already analyzed in Section 4, where (5.3) has an explicit solution as expected as $\lambda = \pm \sqrt{\delta v^* \bar{\gamma}} \, i$, $(B_0, V_0) = (-b^* \bar{\gamma}, \lambda)$, we focus on the age-dependent interaction case. However, instead of analyzing the characteristic equation for the eigenvalues, i.e. (5.3) as a 2-dimensional homogeneous linear system with determinant being equal to zero, we study the eigenvalues of the differential operator related to (5.3). More precisely, the right hand side of the system above has a common factor $\phi(a): = e^{-\lambda a} \left[B_0- b^*\int_{0}^a \gamma(s) e^{\lambda s} \, ds \cdot V_0 \right]$ which is the variation of the constants formula for the inhomogeneous ODE $\phi'(a) + \lambda\, \phi(a) = - b^* \gamma(a) \, V_0$, $\phi(0) = B_0$. Therefore, (5.3) is equivalent to the eigenvalue problem

Finally we have to study the linear differential operator $A\! : \!\text{dom}(A) \subset X \longrightarrow X$ defined by

whose eigenvalues / eigenfunctions, $\lambda \begin{pmatrix} \phi\\ V_0 \end{pmatrix} = A\begin{pmatrix} \phi\\ V_0 \end{pmatrix}$, can be approximated by the eigenvalues / eigenvectors of a suitable matrix, using pseudospectral methods. See ^[21].

Let us remark that thanks to the minimum mortality assumption $\mu(a) \geq \mu_0 > 0$ and the characterization of the growth bound using spectral values (see ^[5], ^[3], ^[19]) we can focus on the spectral bound $s(A)$ for the (un)stability of the coexistence equilibrium $(b^*, v^*)$. Moreover, spectral values of $A$ with real part larger than $-\mu_0$ are actually eigenvalues, i.e., solutions to system (5.4).

For the numerical implementation, we proceed as follows; see ^[24] and Chapter 3 in ^[1]. We consider polynomial interpolation for the functions in the age interval $[0, a^\dagger]$ with $a^\dagger$ the maximum age, $\phi(a_j) \simeq \Phi_j$, with $a_j$, $j = 0, \dots, N$, being the Chebyshev extremal points and $N$ being the discretization index. Let $H$ be the differentiation matrix of dimension $(N+1) \times (N+1)$, and we consider Clenshaw-Curtis quadrature rule with weights over the interval $[-1, 1]$, $w_j$, $j = 0, \dots, N$, $\frac{1}{2}\sum_{j = 0}^{N} w_j = 1$. Then, the finite-dimensional counterpart of the eigenvalue problem (5.4) is as follows:

where $\beta_j^* = \beta_*(a_j)$ and $\gamma_j^* = \gamma_*(a_j)$, functions defined in (5.2). Solving $\Phi_0$ from the first equation, we get to $\Phi_0 = \frac{a^\dagger}{2- \beta_0^* w_0 a^\dagger}\sum_{j = 1}^{N} \beta_j^* \Phi_j w_j$, and plugging into the other equations, we can reduce the linear system to dimension $N+1$:

From the right hand side above, we define the entries of the matrix $A_N$ of size $(N+1) \times (N+1)$ as

with constants $C_1 = \frac{a^\dagger}{2- \beta_0^* w_0 a^\dagger}$ and $C_2 = \frac{\gamma_0^* w_0 a^\dagger}{2- \beta_0^* w_0 a^\dagger}$, which is a finite-dimensional approximation to the linear differential operator $A$ defined in (5.5).

The finite-dimensional eigenvalue problem $\lambda \begin{pmatrix} \hat{\Phi}\\ V_0 \end{pmatrix} = A_N\begin{pmatrix} \hat{\Phi}\\ V_0 \end{pmatrix}$, with $\hat{\Phi} = (\Phi_1, \dots, \Phi_N)$, is solved either by Matlab's $\texttt{eig}$ or $\texttt{eigs}$.

We can check the accuracy of the proposed pseudospectral method by using the exact solution to (5.4) for the case of age-independent interaction, namely, age-independent eigenfunctions corresponding to purely imaginary eigenvalues; see Figure 3. One readily sees from (5.4) that $\phi(a) \equiv 1$ if and only if $\gamma(a) \equiv \bar{\gamma}$. Numerically, one readily sees from (5.6) that $\Phi_j = 1$, $j = 0 \dots N$ if and only if $\gamma_i = \bar{\gamma}$, $i = 0 \dots N$ since $\sum_{j = 0}^{N} H_{i j} = 0$.

Figure 3.  Accuracy of the pseudospectral method (5.6) for increasing discretization index $N$ and constant interaction. Exact solution to (5.4): pair of complex-conjugate roots on the imaginary axis $\lambda_1 = \pm \sqrt{\delta v^* \bar{\gamma}} \, i$, $\phi(a) \equiv 1$, spectral bound $s(A) = \text{Re}(\lambda_1) = 0$ versus rightmost roots $\lambda_{1, N}$ of (5.6) computed through $\texttt{eig}( A_N )$ with matrix $A_N$ in (5.8), which are a pair of complex-conjugate eigenvalues with constant eigenvector. Parameter values: mortality $\mu(a) = \frac{\theta}{a^\dagger -a}$, $\theta = 1/4$, fertility $\beta(a) = \bar{\beta} a(a^\dagger -a)$, $\bar{\beta} = 0.009$, interaction $\gamma(a)\equiv \bar{\gamma} = 1/9$ and $\delta = 5, \alpha = 0.75, a^\dagger = 10$. The accuracy of the method is limited by the computation from (3.5) of the steady-state $v^* = 0.4063275$. Matrix spectral bound $s(A_N) = \text{Re}(\lambda_{1, N}) \sim 10^{-9}$ for $N = 100$.

DownLoad: Full-Size Img PowerPoint

6.   Numerical scheme: Predictor-corrector

In this section we introduce a numerical method to compute the time-evolution of the interacting populations with the aim of checking the results obtained in previous sections. Among many numerical schemes (see chapter 7 in ^[1] and the references therein and also ^[19,34,35]) we work on the integrated form of the present predator-prey model (3.1)–(3.2), paying attention to avoid immortal individuals.

The starting point is to consider equal age and time step-sizes $\Delta a = \Delta t$ to get to

and (3.2) can be rewritten as

We recall that $b(t): = u(0, t) = \int_0^\infty \beta(a)\, u(a, t)\, da$, and $\lim_{a \to \infty} u(a, t) = 0$. Now, we can discretize age (in a finite span $[0, a^\dagger]$, with $a^\dagger$ as the maximal age) and time by the grid points $\{ (a_j, t_n) : 0 \leq j \leq J, 0 \leq n \leq N\}$ with $a_j = j \Delta a$, $t_n = n \Delta t$, $\Delta t = \Delta a = a^\dagger /J$ and use quadrature rules for the integrals involved: $u(0, t) = \int_0^\infty \beta(a)\, u(a, t)\, da\simeq \int_0^{a^\dagger} \beta(a)\, u(a, t)\, da\simeq \sum_{j = 0}^{J} w_j \, \beta(a_j) \, u(a_j, t) \, \Delta a$, with weights $\frac{1}{J}\sum_{j = 0}^{J} w_j = 1$. Isolating, we get

with $\Delta a$ small enough such that $w_0\, \beta(0)\, \Delta a < 1$ in order to assure positivity of the numerical solution. On the other hand, we can consider two approximations in (6.1): $v(t+\Delta t)\simeq v(t) \, e^{-\delta \Delta t + \alpha S(t)\, \Delta t}$ as predictor and $v(t+\Delta t)\simeq v(t) \, e^{-\delta \Delta t + \alpha [S(t)+S(t+\Delta t)]\, \Delta t/2}$ as corrector, both with $S(t) \simeq \sum_{j = 0}^{J} w_j \, \gamma(a_j) \, u(a_j, t) \, \Delta a$.

Finally, we get to the following predictor-corrector explicit numerical scheme for the approximation of the solution of (2.1) at the grid points $u(a_j, t_n)\simeq u_j^n$, $v(t_n)\simeq v^n$, computed from a given initial condition $u(a_j, 0) = u^0(a_j)$, $v(0) = v^0$:

where $\mu_j = \mu(a_j)$, $\gamma_j = \gamma(a_j)$, we set $u(a^\dagger, t_n)\simeq u_J^n = 0$, and we have used the mid-point rule for the integrals of length $\Delta t$. In addition, we can compute the prey population size over time as

and the prey birth-rate $b(t_n) = u(0, t_n) \simeq u_0^n$. Finally, we should choose maximal age $a^\dagger$ and mortality rate $\mu(a)$ to avoid or minimize the probability of immortal individuals, that is, from (3.1) for $t > a^\dagger$, $\frac{u(a^\dagger, t)}{u(0, t-a^\dagger)} \leq \exp[-\int_0^{a^\dagger} \mu(a) da]$ should be zero or small enough.

We can check the accuracy of the predictor-corrector numerical scheme (6.2) using the exact solution for the case of age-independent interaction, a closed periodic orbit; see Theorem 1. Taking an initial condition such that $u^0(a) = P_0\, \bar{\omega}(a) > 0$, $b_0 = P_0\, \bar{\omega}(0) > 0$ and $v_0 > 0$, we can compute the error as $\Vert (b_0, v_0) - (b_T, v_T)\Vert$, where the point $(b_T, v_T)$ is the closest point to the initial condition $(b_0, v_0)$ after a single revolution around the steady state $(b^*, v^*)$; see the inset plot in Figure 4 for an illustration.

Figure 4.  Accuracy of the predictor-corrector numerical scheme (6.2) for constant interaction between predators and preys. Maximal age $a^\dagger = 10$, number of steps in age $J = 5000$ and step size $\Delta t = \Delta a = 0.002$. Comparison between the analytical periodic orbit (4.4) (yellow curve) and the numerical simulation of the orbit (blue dotted curve). Inset plot (zoomed area) showing that the error of the computed curve is less than $10^{-4}$. Parameter values as in Figure 2.

DownLoad: Full-Size Img PowerPoint

7.   Numerical tests

In this section, we present several numerical experiments to investigate how the spectrum of the finite-dimensional operator, the matrix $A_N$ in (5.8), approximates the point spectrum and the spectral bound of the differential operator $A$, defined in (5.5), determining the (un)stability of the coexistence equilibrium.

Let us consider a more realistic case than the age-dependent proportional case (i.e., $\beta(a) = \bar{\beta} \, \mu (a)$ and $\gamma(a) = \bar{\gamma} \, \mu (a)$), by assuming a finite age-span $[0, a^\dagger]$} with unbounded mortality $\mu(a) = \frac{\theta}{a^\dagger -a}$, $\theta > 0$, to avoid immortality in the preys, and fertility $\beta(a) = \bar{\beta} a(a^\dagger -a)$, $\bar{\beta} > 0$ with a maximum at an intermediate age. In this scenario, the computations of previous sections are slightly different. Indeed, the basic reproduction number is computed as $R_0 = \int_0^{a^\dagger} \beta(a) \left(\frac{a^\dagger -a}{a^\dagger}\right)^{\theta}\, da > 1$, and the Malthusian parameter $r > 0$ is the solution to $1 = \int_0^{a^\dagger} \beta(a) \left(\frac{a^\dagger -a}{a^\dagger}\right)^{\theta} e^{-r a}\, da$. Moreover, equation (3.5) becomes $1 = \int_0^{a^\dagger} \beta(a) \left(\frac{a^\dagger -a}{a^\dagger}\right)^{\theta} e^{-\int_0^a \gamma(s) \, ds \cdot v^*}\, da$.

Finally, we can choose an age-dependent interaction rate like

modeling three different age-dependent interactions between predators and preys:

● If $k < 0$, the interaction is stronger at intermediate ages of the preys. $\gamma(a)$ has a local maximum. Predation is focused on the most fertile preys.

● If $k = 0$, the interaction is independent of the age of the prey. $\gamma(a)$ is a constant $\bar{\gamma} > 0$. For predators, the ages of the preys do not matter.

● If $k > 0$, the interaction is stronger at early and late ages of the preys. $\gamma(a)$ has a local minimum. Predation is focused on the less fertile preys.

Notice that increasing parameter $k$ in this interaction rate $\gamma$, from negative to positive, predators change their target preys, from mid-aged preys to both young and old preys.

See Figure 5 for an illustration of these three scenarios with a very different predator-prey dynamics. The numerical computation of the spectral bound of the linearized operator (5.5) around the coexistence equilibrium, is in agreement with the computation of the asymptotic population dynamics (i.e. the time-evolution of the interacting populations). In particular, when an isolated stable periodic orbit arises via a super-critical Hopf bifurcation, the equilibrium becomes unstable.

Figure 5.  Left panels: spectral values associated to the predator-prey equilibrium computed from matrix $A_N$ in (5.8) with $N = 200$ and interaction $\gamma(a) = \bar{\gamma} \big(1-4ka(a^\dagger -a)/{a^\dagger}^2)$, $\bar{\gamma} > 0$, with $k = -1, 0, 0.7$. From top to bottom, the rightmost complex-conjugate spectral values $\lambda_1$ cross the imaginary axis with positive speed, showing a Hopf bifurcation of the steady-state at $k = 0$. Predators change their target preys, from mid-aged preys ($k < 0$) to both young and old preys ($k > 0$). Numerical values: $\lambda_1^{k = -1} = -0.012901772 \pm 0.486759890 i$ $\to$ $\lambda_1^{k = 0} = \pm 0.475118407 i$ $\to$ $\lambda_1^{k = 0.7} = 0.009531135 \pm 0.441038955 i$. Right panels: time-evolution in the plane of prey birth-rate and predator population size, for each case. Parameter values as in Figure 3.

DownLoad: Full-Size Img PowerPoint

8.   Conclusions and discussion

We have used several mathematical methods, either analytical or numerical techniques, to quantify the impact of a predator population targeting selected age-classes (or size-classes) in a prey population.

The idealized ecological scenario of both interacting populations living in harmonic periodic oscillation (with amplitude and frequency depending on initial conditions), takes place only when their interaction is independent of the age/size of the preys; see Figures 2 and 4 and the middle panels in Figure 5. For this unrealistic case, i.e., predators do not take advantage of the preys age or size, PDE system (2.1) with age-dependent parameters cannot be reduced to an ODE system, although its asymptotic behavior is actually given by the limit system of ODEs (4.2), which turns out to be the original Lotka-Volterra model. The key-point of this fact is the computation of the asymptotic age-profile (normalized age-density) for the prey population; see Theorem 1 in Section 4.

However, for natural populations, it is more realistic to assume a predator feeding on preys when they are at certain ages or sizes. For example, smaller preys are easier to catch. In this likely scenario and for a large set of initial conditions, we have an oscillating transition, eventually stabilizing towards a steady-state population; see Figure 1 and the top panels in Figure 5. Yet, in this more realistic ecological scenario, we have also seen that, for instance, when predators focus on both young (smaller) and old (more vulnerable) preys, and preys are most fertile at intermediate ages, predator and prey populations converge to a specific periodic oscillation independently of the initial conditions; see the bottom panels in Figure 5.

With our analysis, we have contributed to the long-standing debate on the oscillations onset in ecological models and on which are the minimal model ingredients able to trigger a Hopf bifurcation of the coexistence steady-state of both species. In summary, we just need to include the preys age in the interaction term of Lotka-Volterra type (i.e., proportional to both populations) to get to a stable periodic dynamics of predator and prey populations.

The numerical results are reliable since we have checked the accuracy of the methods introduced in Sections 5 and 6, using the exact solution available for specific model parameters.

Use of AI tools declaration

The authors declare they have not used artificial intelligence (AI) tools in the creation of this article.

Acknowledgments

Members of the Catalan research group 2021 SGR 00113. This research was funded by Ministerio de Ciencia, Inovación y Universidades, grant numbers PID2019-104437GB-I00 and PID2021-123733NB-I00.

Conflict of interest

The authors declare there is no conflict of interest.