2. The nonconservative kinetic model

Let us consider an interacting system composed of particles which undergo stochastic interactions. The system is divided into $n \in \mathbb{N}$ functional subsystems [Reference Bellomo, Bianca and Delitala32] such that particles belonging to the same functional subsystem share the same strategy. This in particular means that they share the same state, and the microscopic variable (state variable) related to the $i$ th functional subsystem, for $i\in \{1,2, \ldots, n\}$ , is $u_i\in \mathbb{R}$ .

The distribution function on the $i$ th functional subsystem, for $i\in \{1, 2, \ldots, n\}$ , is

\begin{equation*}f_i(t)\,:\,[0,\,T]\rightarrow \mathbb {R}^+, \qquad T\gt 0,\end{equation*}

that measures the number of particles in the $i$ th subsystem, at time $t\gt 0$ . Moreover, the distribution function of the overall system writes

\begin{equation*}\mathbf {f}(t)=\big (\,f_1(t), f_2(t), \ldots,\, f_n(t)\big ).\end{equation*}

The interaction between particles is ruled by the following quantities:

• $\eta _{hk}\geq 0$ , for $h,k \in \{1,2, \ldots, n\}$ , that is the interaction rate between particles in the $h$ th subsystem and particles in the $k$ th subsystem.

• $B_{h,k}^i\geq 0$ , for $i,h,k \in \{1, 2, \ldots, n\}$ , that is the transition probability, i.e. the probability that a particle of the $h$ th subsystem falls into the $i$ th after interacting with a particle of the $k$ th subsystems. Specifically, the following property holds true:

  \begin{equation*}\sum _{i=1}^n B^i_{h,k}=1,\qquad \forall h,k \in \{1, 2, \ldots, n\}.\end{equation*}

The evolution of the $i$ th functional subsystem, for $i \in \{1, 2, \ldots, n\}$ , is described by the following nonlinear ordinary differential equation with quadratic nonlinearity (see [Reference Bianca2, Reference Bertotti and Delitala33] and references therein for derivation and analytical details)

(1) \begin{equation} \begin{split} \frac{df_i}{dt}(t)&=G_i[\mathbf{f}](t)-L_i[\mathbf{f}](t)\\& =\sum _{h,k=1}^n \eta _{hk}\,B^i_{hk}\,f_h(t)f_k(t)-f_i(t)\sum _{k=1}^n\eta _{ik}\,f_k(t). \end{split} \end{equation}

Specifically:

• $G_i[\mathbf{f}](t)\,:\!=\sum _{h,k=1}^n \eta _{hk}\,B^i_{hk}\,f_h(t)\,f_k(t)$ , for $i \in \{1,2, \ldots, n\}$ , is the gain term operator which measures the number of particles that fall into the $i$ th functional subsystem.

• $L_i[\mathbf{f}](t)=f_i(t)\sum _{k=1}^n\eta _{ik}\,f_k(t)$ , for $i \in \{1,2, \ldots, n\}$ , is the loss term operator which measures the number of particles that leave the $i$ th functional subsystem.

The kinetic framework (1) is conservative, i.e.

\begin{equation*}\frac {d\rho }{dt}(t)\,:\!=\sum _{i=1}^n\frac {df_i}{dt}=0,\end{equation*}

where $\rho (t)$ is the density of the overall system. Let $\mathbb{E}_p[\mathbf{f}](t)$ be the $p$ th-order moment of the system, for $p \in \mathbb{N}$ , which defines

\begin{equation*}\mathbb {E}_p[\mathbf {f}](t)\,:\!=\sum _{i=1}^n u_i^p\,f_i(t).\end{equation*}

Then, the density $\rho (t)$ corresponds to the $0$ th-order moment of the system, $\mathbb{E}_0[\mathbf{f}](t)$ . Bearing the definition of this moment in mind, a system is said to be conservative if

\begin{equation*}\frac {d\mathbb {E}_0[\mathbf {f}]}{dt}(t)=0, \qquad \forall t\gt 0.\end{equation*}

To get a nonconservative framework, another set of parameters has to be introduced. According to what was done in [Reference Bianca34], let $\mu _{hk}$ , for $h,k\in \{1, 2, \ldots, n\}$ , be the birth/death rate of particles of the $h$ th functional subsystems due to the encounters with particles of the $k$ th functional subsystem. Then, the nonconservative operator, also known as proliferative/destructive term, is defined, for $i \in \{1, 2, \ldots, n\}$ , as

\begin{equation*}P_i[\mathbf {f}](t)\,:\!=\,f_i(t)\sum _{k=1}^n\eta _{ik}\,\mu _{ik}\,f_k(t).\end{equation*}

Then, the evolution of the $i$ th functional subsystem, for $i \in \{1, 2, \ldots,n\}$ , is described by the following nonlinear ordinary differential equation

(2) \begin{equation} \begin{split} \frac{df_i}{dt}(t)&=G_i[\mathbf{f}](t)-L_i[\mathbf{f}](t)+P_i[\mathbf{f}](t)\\& =\sum _{h,k=1}^n \eta _{hk}\,B^i_{hk}\,f_h(t)\,f_k(t)-f_i(t)\sum _{k=1}^n\eta _{ik}\,f_k(t)+f_i(t)\sum _{k=1}^n\eta _{ik}\,\mu _{ik}\,f_k(t). \end{split} \end{equation}

By summing on $i \in \{1,2, \ldots, n\}$ the (3), one has

\begin{equation*}\frac {d\rho }{dt}(t)=\sum _{i,k=1}^n\eta _{ik}\,\mu _{ik}\,f_i(t)\,f_k(t),\end{equation*}

that is the macroscopic equation of the density of the system.

The equation (2) rewrites in the form

(3) \begin{equation} \frac{df_i}{dt}(t)=\sum _{h,k=1}^n \eta _{hk}\,B^i_{hk}\,f_h(t)\,f_k(t)+f_i(t)\sum _{k=1}^n\eta _{ik}\left (\mu _{ik}-1\right )f_k(t). \end{equation}

If $\mathbf{f}^0=\left (\,f_1^0,f_2^0, \ldots, f_n^0\right )\in (\mathbb{R}^+)^n$ is a suitable initial data, then the Cauchy problem, or initial value problem, related to the nonconservative kinetic framework (2) is assigned

(4) \begin{equation} \begin{cases} \displaystyle \frac{df_i}{dt}(t)=G_i[\mathbf{f}](t)-L_i[\mathbf{f}](t)+P_i[\mathbf{f}](t) &{ t\in [0,\, T]}\\[6pt] \mathbf{f}(0)=\mathbf{f}^0. \end{cases} \end{equation}

3. The forcing term

In the previous section, we have considered a system of interacting particles such that there is no interaction due to the external environment. Let now consider an external force field acting on the system. This force acts in a different way on each functional subsystem. Then, the $i$ th component, for $i \in \{1, 2, \ldots, n\}$ , of the external force field that acts on the system is the function

\begin{equation*}F_i[\mathbf {f}](t)\,:\, [0,\, T] \rightarrow \mathbb {R}.\end{equation*}

The external force field acting on the overall system is described by the vector function

\begin{equation*}\mathbf {F}[\mathbf {f}](t)=\left (F_1[\mathbf {f}](t), \, F_2[\mathbf {f}](t), \ldots, F_n[\mathbf {f}](t)\right ),\end{equation*}

such that

\begin{equation*}\mathbf {F}[\mathbf {f}](t)\,:\,(\mathbb {R}^+ )^n \times [0,\,T]\rightarrow \mathbb {R}^n.\end{equation*}

The first novelty of the current paper lies in the fact that the components of this force field can be negative during the time evolution and depend explicitly on the distribution function $\mathbf{f}(t)$ . Moreover, in this paper a particular shape is considered. Specifically, for $i\in \{1, 2, \ldots, n\},$

(5) \begin{equation} F_i[\mathbf{f}](t)=f_i(t)\,F_i(t), \qquad t\gt 0, \end{equation}

where

(6) \begin{equation} F_i(t)\,:\,[0,\,T]\rightarrow \mathbb{R}. \end{equation}

In the spirit of [Reference Bianca and Mogno35], bearing the previous framework (2) in mind, the evolution of the $i$ th functional subsystem in a nonconservative kinetic framework, where both proliferative/destructive events and an external force field occur, writes

(7) \begin{equation} \frac{df_i}{dt}(t)=G_i[\mathbf{f}](t)-L_i[\mathbf{f}](t)+P_i[\mathbf{f}](t)+f_i(t)\,F_i(t). \end{equation}

In vector form, this new kinetic framework (7) writes

\begin{align*} \begin{split} \frac{d\mathbf{f}}{dt}(t)&=\mathbf{G}[\mathbf{f}](t)-\mathbf{L}[\mathbf{f}](t)+\mathbf{P}[\mathbf{f}](t)+\mathbf{f}(t)\cdot \mathbf{F}(t)\\& =\mathbf{J}[\mathbf{f}](t)+\mathbf{P}[\mathbf{f}](t)+\mathbf{f}(t)\cdot \mathbf{F}(t), \end{split} \end{align*}

where $\mathbf{J}[\mathbf{f}](t)$ represents the conservative part of the framework.

The evolution of the $i$ th functional subsystem, for $i \in \{1, 2, \ldots, n\}$ , is thus influenced by the $i$ th component of the external force field, $F_i(t)$ , weighted by the distribution function $f_i(t)$ . This choice is motivated by the fact that commonly the action of an external force field may depend on the number of particles that are in that particular state, i.e. the current functional subsystem. In other words, we must expect that the effects of the force are stronger and more perceivable, the larger is the number of particles in that state, that is the larger is the value of the distribution function $f_i(t)$ . However, in this perspective, we need to give an explicit dependence of the forcing term on $\mathbf{f}$ , and we have chosen to assume a linear dependence, with $\mathbf{F}$ depending only on the time. On one hand, this seems quite plausible, at least in this initial phase of study, since it allows to grasp at least some very interesting non-mechanical phenomena. But also, from a mathematical viewpoint, the choice of the term $f_i(t)\,F_i(t)$ ensures the well-posedness of the problem. As matter of fact if $f_i(\,\bar{t}\,)=0$ and $F(\bar{t})\lt 0$ for some $\bar{t}\gt 0$ , then the distribution function $f_i(t)$ still remains nonnegative for $t\gt \bar{t}$ .

By using the equation (3), the kinetic framework (7) rewrites, for $i \in \{1,2, \ldots, n\}$ ,

(8) \begin{equation} \frac{df_i}{dt}(t)=\sum _{h,k=1}^n \eta _{hk}\,B^i_{hk}\,f_h(t)\,f_k(t)+f_i(t)\sum _{k=1}^n\eta _{ik}\left (\mu _{ik}-1\right )f_k(t)+f_i(t)\,F_i(t). \end{equation}

4. Analytical results

This section aims to prove some analytical results about the nonconservative kinetic system (7) (or in its extended form (8)). Specifically, the next result shows the existence and uniqueness of solutions of the related Cauchy Problem.

Theorem 1. Assume that

1. (i) There exists $\eta \gt 0$ such that $\eta _{hk}\leq \eta$ , $\forall h,k \in \{1, 2, \ldots, n\}$ .

2. (ii) There exists $\mu \gt 0$ such that $|\mu _{hk}|\leq \mu$ , $\forall h,k\in \{1,2, \ldots, n\}$ .

3. (iii) There exists $F\gt 0$ such that $|F_i(t)|\leq F$ , $\forall i \in \{1,2, \ldots, n\}$ and $\forall t\gt 0$ .

Let $\mathbf{f}^0=\left (\,f^0_1, f^0_2, \ldots, f^0_n\right )\in (\mathbb{R}^+)^n$ , such that $\sum _{i=1}^nf^0_i=1$ . Then, there exists a unique local positive solution $\mathbf{f}(t)\in \left (C\left ([0,\, T]\right )\right )^n$ , for $T\gt 0$ , of the Cauchy Problem

(9) \begin{equation} \begin{cases} \displaystyle \frac{df_i}{dt}(t)=G_i[\mathbf{f}](t)-L_i[\mathbf{f}](t)+P_i[\mathbf{f}](t)+f_i(t)\,F_i(t) &{ i \in \{1,2, \ldots, n\},\, t\gt 0}\\[6pt] \mathbf{f}(0)=\mathbf{f}^0. \end{cases} \end{equation}

Moreover, if

(10) \begin{equation} \sum _{i=1}^n \int _0^{+\infty } f_i(\tau )\left (\sum _{k=1}^n\eta _{ik}\mu _{ik}\,f_k(\tau )+F_i(\tau )\right )\, d\tau \lt +\infty, \end{equation}

then the solution exists globally in time.

Proof. The local existence of solutions of the Cauchy Problem (9) is ensured by using Lipschitz arguments. Specifically, let us consider the operator defined by the right-hand side of the equation (8), i.e. for $i\in \{1,2, \ldots, n\}$ ,

\begin{equation*}\mathbb {L}_i[\mathbf {f}](t)\,:\!=\sum _{h,k=1}^n \eta _{hk}\,B^i_{hk}\,f_h(t)\,f_k(t)+f_i(t)\sum _{k=1}^n\eta _{ik}\left (\mu _{ik}-1\right )f_k(t)+f_i(t)\,F_i(t).\end{equation*}

By using assumptions (i)–(iii), it is possible to conclude that this operator is Lipschitz. Precisely, let $\mathbf{f}(t),\mathbf{g}(t)\in C\left ([0,T]\right )$ , then, for $i \in \{1,2, \ldots, n\}$ ,

(11) \begin{equation} \begin{split} \big |\mathbb{L}_i[\mathbf{f}]-\mathbb{L}_i[\mathbf{g}]\big |&=\Bigg |\sum _{h,k=1}^n \eta _{hk}\,B^i_{hk}\,f_h(t)\,f_k(t)+f_i(t)\sum _{k=1}^n\eta _{ik}\left (\mu _{ik}-1\right )f_k(t)+f_i(t)\,F_i(t)\\&\quad -\sum _{h,k=1}^n \eta _{hk}\,B^i_{hk}\,g_h(t)g_k(t)+g_i(t)\sum _{k=1}^n\eta _{ik}\left (\mu _{ik}-1\right )g_k(t)+g_i(t)\,F_i(t)\Bigg |\\& \leq \Bigg | \sum _{h,k=1}^n\eta _{hk}B^i_{hk}\big (\,f_h(t)\,f_k(t)-g_h(t)g_k(t)\big )+\sum _{k=1}^n\eta _{ik}\left (\mu _{ik}-1\right )\big (\,f_i(t)\,f_k(t)-g_i(t)g_k(t)\big )\\&\quad +F_i(t)\left (\,f_i(t)-g_i(t)\right )\Bigg |\\& \leq \eta \sum _{h=1}^nf_h(t)\sum _{k=1}^nB^i_{hk}\big |\,f_k(t)-g_k(t)\big |+\eta \sum _{h=1}^ng_k(t)\sum _{k=1}^nB^i_{hk}\big |\,f_h(t)-g_h(t)\big |\\&\quad +\eta (\mu -1)\,f_i(t)\sum _{k=1}^n\big |\,f_k(t)-g_k(t)\big |+\eta (\mu -1)\big |\,f_i(t)-g_i(t)\big |\sum _{k=1}^ng_k(t)\\& \quad F\left |f_i(t)-g_i(t)\right |. \end{split} \end{equation}

By using the (11), one has

(12) \begin{equation} \begin{split} \big \|\mathbb{L}[\mathbf{f}]-\mathbb{L}[\mathbf{g}]\big \|_1&=\sum _{i=1}^n\big |\mathbb{L}_i[\mathbf{f}]-\mathbb{L}_i[\mathbf{g}]\big |\\& \leq \eta \sum _{i=1}^n\sum _{h=1}^n\,f_h(t)\sum _{k=1}^nB^i_{hk}\big |\,f_k(t)-g_k(t)\big |+\eta \sum _{i=1}^n\sum _{h=1}^ng_k(t)\sum _{k=1}^nB^i_{hk}\big |\,f_h(t)-g_h(t)\big |\\& \quad+\eta (\mu -1)\sum _{i=1}^n\,f_i(t)\sum _{k=1}^n\big |\,f_k(t)-g_k(t)\big |+\eta (\mu -1)\sum _{i=1}^n\big |\,f_i(t)-g_i(t)\big |\sum _{k=1}^ng_k(t)\\& \quad +F\sum _{i=1}^n\big |\,f_i(t)-g_i(t)\big |\\& \leq L\|\mathbf{f}(t)-\mathbf{g}(t)\|_1, \end{split} \end{equation}

where the constant $L$ depends on the parameters of the kinetic system, i.e. $F, \eta, \mu$ . Finally, the operator $\mathbb{L}[\mathbf{f}](t)$ is proved to be Lipschitz by passing to the maximum in relations (12). The local existence of a solution $\mathbf{f}(t)$ of the (9) is thus gained.

In order to prove the positivity of the solution $\mathbf{f}(t)$ , equation (8) needs to be written in the form

(13) \begin{equation} \frac{df_i}{dt}(t)=\mathcal{Q}_i[\mathbf{f}](t)+f_i(t)\mathcal{S}_i[\mathbf{f}](t), \end{equation}

where

\begin{align*} \mathcal{Q}_i[\mathbf{f}](t)&\,:\!=\sum _{h,k=1}^n \eta _{hk}\,B^i_{hk}\,f_h(t)\,f_k(t)\\ \mathcal{S}_i[\mathbf{f}](t)&\,:\!=\sum _{k=1}^n\eta _{ik}\left (\mu _{ik}-1\right )f_k(t)+F_i(t). \end{align*}

Let now

\begin{equation*}\gamma _i(t)\,:\!=\int _0^t-\mathcal {S}_i[\mathbf {f}](\tau )\,d\tau .\end{equation*}

Then, equation (13) yields

(14) \begin{equation} f_i(t)=f^0_ie^{-\gamma _i(t)}+\int _0^te^{\gamma _i(\tau )-\gamma _i(t)}\,\mathcal{Q}_i[\mathbf{f}](\tau )\,d\tau \,, \end{equation}

for $i \in \{1,2, \ldots, n\}$ .

Due to positivity of the initial data $\mathbf{f}^0$ and exponential function, positivity of the solution $\mathbf{f}(t)$ of the Cauchy Problem (9) is seen at once.

Bearing expression (8) in mind and integrating it on $[0,t]$ , for $t\gt 0$ , one has

(15) \begin{equation} f_i(t)=f_i^0+\int _0^t\left (\sum _{h,k=1}^n \eta _{hk}\,B^i_{hk}\,f_h(\tau )\,f_k(\tau )+f_i(\tau )\sum _{k=1}^n\eta _{ik}\left (\mu _{ik}-1\right )f_k(\tau )+f_i(\tau )\,F_i(\tau )\right )\, d\tau. \end{equation}

The global existence of solution would be violated if there were a time $\bar{t}\gt 0$ such that $|\mathbf{f}|\rightarrow +\infty$ as $t \rightarrow \bar{t}$ . If we sum relations (15) on $i=1,2,\ldots, n$ , then we get

(16) \begin{equation} \begin{split} \sum _{i=1}^nf_i(t)&=\sum _{i=1}^nf_i^0+\sum _{i=1}^n\left (\int _0^t\left (\sum _{h,k=1}^n \eta _{hk}\,B^i_{hk}\,f_h(\tau )\,f_k(\tau )+f_i(\tau )\sum _{k=1}^n\eta _{ik}\left (\mu _{ik}-1\right )f_k(\tau )+f_i(\tau )\,F_i(\tau )\right )\, d\tau \right )\\& =1+\sum _{i=1}^n \int _0^t f_i(\tau )\left (\sum _{k=1}^n\eta _{ik}\mu _{ik}\,f_k(\tau )+F_i(\tau )\right )\, d\tau. \end{split} \end{equation}

The assumption (10) and the equation (16) ensure that the above-mentioned condition cannot occur. Then, the global existence of a unique positive solution $\mathbf{f}(t)\in C\left (0,+\infty \right )$ of the Cauchy Problem (9) is proved.

5. Macroscopic equations and stationary state

Theorem 1 of previous section proves the well-posedness of the new kinetic framework (8). In order to recover the macroscopic equations, let sum on $i =1,2, \ldots, n$ the (8), then

(17) \begin{equation} \frac{d\rho }{dt}(t)=\sum _{i=1}^n\sum _{k=1}^n\eta _{ik}\mu _{ik}\,f_i(t)\,f_k(t)+\sum _{i=1}^n\,f_i(t)\,F_i(t), \end{equation}

where the conservative part of the system has disappeared. The macroscopic equation (17) describes the time evolution of the density of the whole system, i.e. $\rho (t)$ . Bearing assumptions of Theorem 1 in mind, one has

\begin{equation*}\rho (0)=1.\end{equation*}

Then, the macroscopic dynamics of this new kinetic framework with force term is defined by the following Cauchy Problem

(18) \begin{equation} \begin{cases} \displaystyle \dot{\rho }(t)=\sum _{i=1}^n\sum _{k=1}^n\eta _{ik}\mu _{ik}\,f_i(t)\,f_k(t)+\sum _{i=1}^nf_i(t)\,F_i(t) &{ t\gt 0}\\[6pt] \rho (0)=1. \end{cases} \end{equation}

A stationary solution of the kinetic framework (8) is any function $\mathbf{g}\in \left (\mathbb{R}^+\right )^n$ which is a solution of the related stationary problem

(19) \begin{equation} \sum _{h,k=1}^n \eta _{hk}\,B^i_{hk}\,g_hg_k+g_i\sum _{k=1}^n\eta _{ik}\left (\mu _{ik}-1\right )g_k+f_i\,F_i(t)=0. \end{equation}

In general, the explicit form of stationary solutions to the kinetic system is not explicitly known. But sometimes it is possible to prove analytical results about existence and uniqueness of stationary solutions. In this paper, the contents of these results will be shown from a numerical point of view, by providing some numerical simulations for a specific ecological problem.

6. An application inspired to ecology

This section aims at presenting some numerical experiments towards the new kinetic framework (8). In particular, a predator–prey model is considered (for details see [Reference Berryman36–Reference Wangersky39] and references therein). The strong, sometimes devastating, effects of the forcing term on the behaviours and even on the existence of two populations will appear clearly from the examples below. It is worth noting that these are numerical experiments.

In what follows, three different scenarios are investigated, with respect to particular values of interaction rates, transition probabilities, birth/death rates and force field. These three scenarios present three different shapes for components $F_i(t)$ of the external force field. As already stated before, this external action may represent, among others, resource availability, rainfall levels and climate change. For instance, a large availability of resources, from the environment, can increase populations, as well as poor water levels can reduce a population far beyond natural mortality.

Firstly, the ecological system is divided into two functional subsystems. Accordingly, two distribution are introduced:

• $f_1(t)$ for the group of prey;

• $f_2(t)$ for the group of predators.

In what follows the parameters acquire a specific meaning:

• The interaction rates $\eta _{hk}$ , for $h,k\in \{1, 2\}$ , give the number of encounters between couple of elements of each functional subsystem.

• The transition probabilities $B^i_{hk}$ , for $i,h,k \in \{1, 2\}$ , and the birth/death rates $\mu _{hk}$ , for $h,k \in \{1, 2\}$ , describe the kinds of pairwise interactions between elements of each functional subsystem.

• The external force field $F_i(t)$ , for $i \in \{1, 2\}$ , describes the action of the environment on the $i$ th population.

In what follows, the current values of these quantities will be defined for each scenario.

6.1. Exponential decay of force term

In this first scenario, the interaction rate is

(20) \begin{equation} \eta _{hk}=0.3 \qquad \forall h,k \in \{1, 2\}. \end{equation}

And the birth/date rate is

(21) \begin{equation} \mu _{hk}=\left (\begin{matrix} 0.03 &&-0.01\\ 0.01 &&-0.03 \end{matrix} \right ). \end{equation}

Specifically:

• $\mu _{11}$ is positive, then a birth rate is considered for the prey;

• $\mu _{22}$ is negative, then a death rate is considered for the predators.

• $\mu _{12}$ is negative and $\mu _{21}$ is positive ( $\mu _{12}=-\mu _{21}$ ) since predators may feed on prey.

Moreover, transition probabilities are

(22) \begin{equation} B^i_{11}=\left (\begin{matrix} 1&& 0\end{matrix}\right ), \quad B^i_{12}=\left (\begin{matrix} 0.5&& 0.5\end{matrix}\right ), \quad B^i_{21}=\left (\begin{matrix} 0&& 1\end{matrix}\right ), \quad B^i_{22}=\left (\begin{matrix} 0.2&& 0.8\end{matrix}\right ). \end{equation}

In particular:

• $B^i_{11}=\left (\begin{matrix} 1&& 0\end{matrix}\right )$ means that when a prey encounters another prey, then it does not change functional subsystem.

• $B^i_{12}=\left (\begin{matrix} 0.5&& 0.5\end{matrix}\right )$ models the fact that when a prey encounters a predator, then it can be eaten with a certain probability.

• $B^i_{21}=\left (\begin{matrix} 0&& 1\end{matrix}\right )$ means that when a predator encounters another one, then it does not change functional subsystem.

• $B^i_{22}=\left (\begin{matrix} 0.2&& 0.8\end{matrix}\right )$ models the fact that any pairwise interaction between prey may reduce the related population with a certain probability.

Finally, the force field has the following shape

\begin{equation*} \mathbf {F}[\mathbf {f}](t)=\left (e^{-t}f_1(t),\, e^{-t}f_2(t)\right ), \end{equation*}

where the related components (6) read

\begin{align*} F_1(t)&=e^{-t}\\ F_2(t)&=e^{-t}. \end{align*}

In this scenario, there is an exponential decay of the external force field. Roughly speaking, it is a very fast action on the evolution of the system of a plausibly monotonic variation of environment, for instance the exhaustion of some basic resources. In particular, this model seems to describe rather faithfully the gradual but constant expansion of the so-called ‘Sahel belt’, the region of desertification on the borders of Sahara, where the progressive depletion of the water led to the gradual but fast extinction first of herbivores and then of lions.

In fact, given the initial data $\mathbf{f}^0=\left (0.5,\, 0.5\right )$ , the dynamics of the system is defined by the Cauchy problem (9). Figure 1 shows the two solutions. In particular, after the peak of prey, there is the peak of predators. Nevertheless, the action of the environment determines the extinction of the prey and the consequent extinction of predators.

Figure 1. The solution $f_1(t)$ , red curve of prey and $f_2(t)$ , blue curve of predators, with binary interactions furnished by (20), (21) and (22), and exponential external force field $\mathbf{F}=(e^{-t}, \, e^{-t})$ with exponential decay.

Remark 3. Once we have concluded that the exhaustion of common basic resources leads to the extinction of both prey and predators, it is to be expected that the ‘velocity of exhaustion’ is immaterial. In fact, the evolution of the system is not so different with respect to the previous case of exponential decay if the force field is assumed to decay much more slowly as $t\to +\infty$ . In Figure 2, the following cases are considered:

Figure 2. The solution $f_1(t)$ , red curve of prey and $f_2(t)$ , blue curve of predators, with binary interactions furnished by (20), (21) and (22), and two different non-exponential decay external force field. From left to right: polynomial decay $\mathbf{F}=\left(\frac{t}{t^2+1}, \, \frac{t}{t^2+1}\right)$ and square root decay $\mathbf{F}=\left(\frac{1}{2+\sqrt{t}}, \, \frac{1}{2+\sqrt{t}}\right)$ .

1. (i) polynomial decay:

   \begin{align*} F_1(t)&=\frac{t}{t^2+1}\\[3pt] F_2(t)&=\frac{t}{t^2+1}. \end{align*}

2. (ii) square root decay:

   \begin{align*} F_1(t)&=\frac{1}{2+\sqrt{t}}\\[3pt] F_2(t)&=\frac{1}{2+\sqrt{t}}. \end{align*}

The shape of solutions is very similar to the one of exponential decay, though the velocity of extinction is obviously different, as in the considered cases the solutions go to zero much more slowly (compare Figure 1 with Figure 2 $_{1,2}$ ).

6.2. Periodic force term

In a second scenario, the interaction rate is

(23) \begin{equation} \eta _{hk}=0.01 \qquad \forall h,k \in \{1, 2\}. \end{equation}

Then, the birth/death rate is identically zero, i.e. in this scenario there are neither births nor deaths due to pairwise interactions between individuals of the system. Therefore, there are only conservative interactions.

The transition probabilities are now

(24) \begin{equation} B^i_{11}=\left (\begin{matrix} 1&& 0\end{matrix}\right ), \quad B^i_{12}=\left (\begin{matrix} 0.6&& 0.4\end{matrix}\right ), \quad B^i_{21}=\left (\begin{matrix} 0&& 1\end{matrix}\right ), \quad B^i_{22}=\left (\begin{matrix} 0.4&& 0.6\end{matrix}\right ). \end{equation}

That is, only interactions between a predator and a prey can modify the system. In particular, the predators can feed on prey, and the prey can only reproduce.

Furthermore, components (6) of the external force field are

\begin{equation*} F_1(t) = 0.05\,\cos \left (\frac {2\pi t}{100}\right ) \end{equation*}

and

\begin{equation*} F_2(t) = 0.05\,\cos \left (\frac {2\pi t}{100}\right ). \end{equation*}

In other words, the external force field is periodic. This may be taken as a good expression of seasonal changes, leading to alternate increases and decreases in water (and prey food) reserves; any periodic climate change may be described in this way.

Assuming the initial data $\mathbf{f}^0=\left (0.9,\, 0.1\right )$ , the evolution of the related system (9) is shown in Figure 3. In particular, there is a periodic-oscillating shape of solutions. During the first part of iterations, the peak of predator group has a time delay with respect to the peak of prey group, as one expects. But the action of a periodic external force field turns out to produce the complete overlapping of the two solutions after a certain time. This is rather surprising, especially in connection with the interpretation of the forcing term as the periodic turning of the seasons. In the context of competition between predators and prey, a time delay between the oscillations is expected under the action of any external force.

Figure 3. The solution $f_1(t)$ , red curve of prey and $f_2(t)$ , blue curve of predators, with conservative interactions furnished by (23) and (24), without birth/death rates, and a periodic external force field $\mathbf{F}=\left(0.05\,\cos \left(\frac{2\pi t}{100}\right), \, 0.05\,\cos \left(\frac{2\pi t}{100}\right)\right)$ .

Remark 4. If the external force field has a periodic and time-decay shape, then the behaviour of solutions is shown in Figure 4 .,

\begin{align*} F_1(t)&=e^{-0.005t}\,0.05\,\cos \left (\frac{2\pi t}{100}\right )\\[3pt] F_2(t)&=e^{-0.005t}\,0.05\,\cos \left (\frac{2\pi t}{100}\right ). \end{align*}

Some oscillations are preserved for a certain time, after which the solution collapses to zero, i.e. goes to extinction.

Figure 4. The solution $f_1(t)$ , red curve of prey and $f_2(t)$ , blue curve of predators, with conservative interactions furnished by (23) and (24), without birth/death rates, and a periodic external force filed, with exponential decay, $\mathbf{F}=\left(e^{-0.005t}\,0.05\,\cos \left(\frac{2\pi t}{100}\right), \, e^{-0.005t}\,0.05\,\cos \left(\frac{2\pi t}{100}\right)\right)$ .

6.3. A mixed case

The last scenario analysed is very close to the previous one, as shown in Figure 3. Nevertheless, the birth/death rate is not zero, i.e. there are nonconservative interactions between agents of this system, i.e.

(25) \begin{equation} \mu _{hk}=\left (\begin{matrix} 0.2 && -0.2\\ 0 && -0.1 \end{matrix} \right ). \end{equation}

There is a birth rate for prey and a death rate for predators. Moreover, the condition $\mu _{12}\lt 0$ means that predators can eat prey during an encounter.

The encounter rate reads

(26) \begin{equation} \eta _{hk}=0.01, \qquad \forall h,k \in \{1, 2\}. \end{equation}

The transition probability is

(27) \begin{equation} B^i_{11}=\left (\begin{matrix} 1&& 0\end{matrix}\right ), \quad B^i_{12}=\left (\begin{matrix} 0.6&& 0.4\end{matrix}\right ), \quad B^i_{21}=\left (\begin{matrix} 0&& 1\end{matrix}\right ), \quad B^i_{22}=\left (\begin{matrix} 0&& 1\end{matrix}\right ). \end{equation}

The only non-trivial probability is $B^i_{12}$ that is related to the fact that a prey can be eaten, with a certain probability, by a predator during an encounter.

Finally, the external force field preserves the periodic structure as well as in the previous example, i.e.

\begin{align*} F_1(t)&=0.05\,\cos \left (\frac{2\pi t}{100}\right )\\[3pt] F_2(t)&=0.05\,\cos \left (\frac{2\pi t}{100}\right ). \end{align*}

The initial data are $\mathbf{f}^0=\left (0.9,\,0.1\right )$ . The simulation towards the related Cauchy problem (9) is shown in Figure 5. There is a sort of mix-structure if it is compared with the one characterised by exponential decay (Figure 1) and periodic shape (Figure 3). In particular, there is the typical oscillation of ‘predator–prey model’. There is the exponential decay due to the extinction. Moreover, there is a sort of big oscillation due to the small oscillations of the system. Roughly speaking, there are two different oscillations related to different time scale.

Figure 5. The solution $f_1(t)$ , red curve of prey and $f_2(t)$ , blue curve of predators, with binary interactions furnished by (25), (26) and (27), and periodic external force field $\mathbf{F}=\left(0.05\,\cos \left(\frac{2\pi t}{100}\right), \, 0.05\,\cos \left(\frac{2\pi t}{100}\right)\right)$ .