A nonconservative kinetic framework under the action of an external force field: Theoretical results with application inspired to ecology

Abstract

The present paper deals with the kinetic-theoretic description of the evolution of systems consisting of many particles interacting not only with each other but also with the external world, so that the equation governing their evolution contains an additional term representing such interaction, called the ‘forcing term’. Firstly, the interactions between pairs of particles are both conservative and nonconservative; the latter represents, among others, birth/death rates. The ‘forcing term’ does not express a ‘classical’ force exerted by the external world on the particles, but a more general influence on the effects of mutual interactions of particles, for instance, climate changes, that increase or decrease the different agricultural productions at different times, thus altering the economic relationships between different subsystems, that in turn can be also perturbed by stock market fluctuations, sudden wars, periodic epidemics, and so on. Thus, the interest towards these problems moves the mathematical analysis of the effects of different kinds of forcing terms on solutions to equations governing the collective (that is statistical) behaviour of such nonconservative many-particle systems. In the present paper, we offer a study of the basic mathematical properties of such solutions, along with some numerical simulations to show the effects of forcing terms for a classical prey–predator model in ecology.

1. Introduction

Kinetic-theoretic models play nowadays an increasingly important role in the description and prediction (at least in stochastic terms) of the behaviour and evolution of systems consisting of a large number of subjects of any nature, see [1–5] and references therein. For the simple good reason that these systems are a generalised version of the systems of material particles that typically represent gases in the classical Boltzmann’s Kinetic Theory [6]. The main difference between classical systems (gases) and general systems is twofold. On one hand, the objects that belong to these latter need not to be material particles, but can be cells [7–9], animals [10], human beings [11, 12] or even artificial or abstract structures, for example the vehicular traffic dynamics [13, 14]. Therefore, the variables describing their state at each instant are not only position and velocity but also, in principle, biological activation [15], wealth [16, 17] and also opinion [18, 19]. On the other hand, the interactions between members of a general system are not deterministically connected to their states, and the final effect of any number of different interactions is not a linear function of their separate effects like in classical mechanics. But, once we have realised that these differences can be reduced to different interpretations of the same mathematical scheme, then the scheme itself can be applied to the widest range of practical problems till now tackled. As a consequence, the research about mathematical features and possible applications of kinetic-theoretic language has been flourishing more and more for about the last forty years. Among others, there are recent developments towards the application of this framework to mathematical epidemiology, due to COVID-19 pandemic event, see [20–24] and references therein.

The systems under consideration may be conservative or nonconservative, isolated or interacting with the external world. A system is said to be conservative if not one of its possible evolutions may result in an increase or decrease in the number of its elements, whereas it is said to be isolated if only interactions between its particles are allowed, and no interactions with objects not belonging to the system are taken into account. Most of the researches about the use of the language, the notions and the methods of Kinetic Theory to describe many-particle systems have concerned isolated conservative systems, since they represent a large number of phenomena that are interesting in themselves and often reflect serious problems whose solution is overwhelmingly important for the present human life. Nevertheless, nonconservative systems can deliver a more faithful picture of reality in many cases, since birth (entrance) and death (exit) of individuals are the rules in almost all systems. On the other hand, even in classical problems like the dynamics of interactions between predators and prey or the dynamics of buying-and-selling between human beings from different social and economic classes, external influences like the alternation of seasons and of heavy rains and droughts or emergency states (wars and epidemics) can play an important role. Accordingly, this paper aims at studying some mathematical features of the solutions to equations governing the evolution of such systems in the kinetic-theoretic scheme in the presence not only of a term describing birth/death effects but also of an additional forcing term $\mathbf{F}$ . It is worth pointing out that the action of external world on the evolution of a phenomenon has already been considered in some nonkinetic situations. Among others, the interested reader may refer to [25] and references therein, for an application towards mathematical epidemiology.

The forcing term, considered in this paper, aims at describing the effects of external influences on the evolution of the whole stochastic interacting system. This action does not depend on the binary interactions between pairs of agents. Moreover, it is important to notice that this term is not a force in the classical sense of word, but it rather implicitly expresses an external influence on the frequencies of interactions and their consequences. Accordingly, as will be shown and discussed in Section 3, it must in general depend on the distribution function $\mathbf{f}$ on the state space (i.e. the set of the values of the state variables of the system), that is $\mathbf{F}[\mathbf{f}](t)$ . So far, when an external force field, acting on a stochastically interacting system, has been considered, it has had only a dependence with respect to the time and occurred not in combination with nonconservative binary interactions; this is a limiting aspect for some applications, if one wants to gain an increasingly realistic description. These effects will be also illustrated in some special cases by means of numerical simulations in Section 6. These numerical experiments refer to an ecological systems, composed of interacting prey and predators. Specifically, the external force field, depending also on distribution functions, models the action of the external environment on these populations, as well as, among others, resource availability, rainfall levels and climate change. Definitely, this paper presents a new kinetic model where two nonconservative aspects act simultaneously on the whole system. On one hand, we consider binary and stochastic nonconservative interactions, in terms of birth/death events. On the other, an external force field occurs, which depends on distribution functions, and not only on time. This new modelling choice, with the derivation of related kinetic equations, seems to be instructive and of particular interest in the present context, since it allows to gain a rather faithful picture of real phenomena.

The introduction of an external force field $\mathbf{F}[\mathbf{f}](t)$ over a system may result in the loss of some analytical properties. For example, the solution, if it exists and it is unique, may be not positive or unbounded. Indeed, nonconservative structure may cause blow-up phenomena [26, 27]. Recently, thermostatted kinetic framework has been introduced for this aim (see [2, 28–31] and references therein). In order to avoid the loss of positivity of solution, a particular shape of the external force field $\mathbf{F}[\mathbf{f}](t)$ has been considered in this paper. Nevertheless, this shape does not ensure the boundedness too, then this new model requires some further analytical assumptions.

The contents of the paper are distributed as follows. Section 2 is devoted to recall the basic kinetic-theoretic language and to write down equations governing the evolution of nonconservative systems. The forcing term is introduced in Section 3, and some basic results about the well-posedness of the system of these equations are given in Section 4. Section 5 treats the formulation of macroscopic equations and the existence of stationary solutions. In Section 6, we perform numerical simulations inspired by an ecological system, modelled by using the new kinetic framework. In particular, we consider external actions decaying with time, at various different rates (Section 6.1) and periodic (Section 6.2). Finally, in Section 7, some conclusions are drawn from both results and simulations, with possible future developments of the present research.

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 [32] 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 [2, 33] 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 [34], 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 [35], 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.

Remark 1. The assumption on the density of the initial data $\mathbf{f}^0$ in Theorem 1, i.e. $\sum _{i=1}^nf^0_i=1$ , is necessary for the proof. However, it is not important the value 1, but the boundedness.

Remark 2. The solution $\mathbf{f}(t)$ of the Cauchy Problem (9), ensured by Theorem 1, can grow without limits during the evolution of the system. This aspect can be in contradiction with intuition since the resources are not infinity. This aspect is regulated by the proliferative term $\mathbf{P}[\mathbf{f}](t)$ and the shape of the force field $\mathbf{F}(t)$ , which represent the nonconservative part of this new kinetic framework (8).

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 [36–39] 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)$ .

7. Conclusions and perspectives

In this paper, we have presented and discussed a discrete kinetic model, i.e. the system (7), where nonconservative interactions occur. Specifically, binary nonconservative interactions, i.e. birth/death processes, are modelled by parameters $\mu _{hk}$ , for $h,k \in \{1, 2, \ldots, n\}$ , whereas the external action by a suitable function, called external force field, $\mathbf{F}[\mathbf{f}](t)$ . Some preliminary analytical results have been furnished, and numerical experiments have shown the application of this model in ecological context.

The novelty is twofold. On one hand, this model presents both nonconservative binary interactions and external force field. To the best of authors’ knowledge, this represents a first attempt in this direction. Indeed, some previous attempts have been made for birth/death processes (see, for instance, [40]), but not combined with the action of external environment. On the other hand, the external force field depends explicitly on the distribution function $f_i(t)$ , for $i\in \{1, 2, \ldots, n\}$ , even compared to what was done in previous works. Indeed, at the current state of art, there are kinetic models that consider just a time-dependent force, i.e. it does not explicitly depend on the distribution functions of the state. Therefore, the shape (5) of $\mathbf{F}[\mathbf{f}](t)$ is such that if the $i$ th functional subsystem is empty, then the external force field does not act on it, i.e. $F_i[\mathbf{f}](t)=0$ . Among others, this choice is related to the fact that consequences of action of an external environment on a system depend on the current state of the system itself. Furthermore, $F_i[\mathbf{f}](t)$ , for $i \in \{1, 2, \ldots, n\}$ , can be also negative, that is they are also allowed to lower the number of elements of each functional subsystem. It is worth stressing that numerical experiments of Section 6 show that this model is not only a theoretical curiosity, nonetheless of interest, but it is driven by the aim of a realistic depiction of an interacting system. For instance, as simulations show, the effects of external phenomena on a system composed of prey and predators depend on the number of individuals that can ‘perceive’ them.

It is however the interpretation of the model in view of a faithful description of the interactions of a system with the external world that seems to deserve a special attention. In this connection, let us have a closer look at the simulations above. We have simulated here the evolution under different scenarios: forces vanishing as $t\to 0$ , no matter at which rate, combined with nonconservative binary interactions; periodic external forces, also in the case of rapidly decreasing amplitude, with only conservative binary interactions; periodic external forces with nonconservative binary interactions. As we have seen in Section 6, the presence of the forcing term can give a good account of the behaviour of an animal system subjected to consumption of external resources (for instance, water, by global heating; in this connection, see Section 6.1) or seasonal climate changes with alternate rainy and dry periods, that increase or decrease the amount of food at disposal of prey (see [41–43] and references therein). It will be noticed that the plots obtained are in accordance with what we would have expected for nonconservative systems and/or in the cases in which the external force is assumed to vanish (at any rate) after a sufficiently long time. The other cases show instead rather unexpected evolution, in which not only the time delay between two subsequent peaks of prey and predators respectively is absent but also (and this is more plausible) no ‘extinction’ occurs. A comparison of the prescriptions on parameters in all cases shows the strong dependence of the solutions on them, [44, 45]. This is a strong thrust towards an accurate and as exhaustive as possible analysis of the dependence of solutions on parameters, which is the first important perspective of development of the present research.

It is worth pointing out that the new kinetic framework (7) may represent a general approach for modelling nonconservative interacting systems, where birth/death events and external actions cannot be neglected if a more realistic description is the aim. Therefore, it could represent a suitable model for describing nonconservative systems, that stochastically interact, in different contexts. From the authors’ purposes, as near future perspective, we trust to model the economic behaviour of a human society influenced by the fluctuations of stock markets. Indeed, the scenarios of Section 6.2, where a periodic force term occurs, may have, at least, two economic interpretations. Firstly, such periodic force is suitable to describe financial fluctuations in stock markets, leading to alternate increases and decreases in the money at disposal of commercial activities. Nevertheless, the periodicity is however too regular to depict faithfully this case: financial fluctuations vary rapidly with time in both frequency and amplitude. Furthermore, we have to assume that the two interacting populations of prey and predators are assumed to be traders and consumers. Then, results of numerical simulations show that an increased availability of money also leads to an increase of the number of consumers but also to the opening of a greater number of stores, and conversely, a smaller amount of money in circulation will result in a smaller number of both consumers and traders. Of course, these considerations need to be further investigated, both mathematically and economically. Moreover, for this aim, the use of a vector activity variable [46, 47] appears more appropriate.

The model and the mathematical results offered in the present paper have to be considered as the first step of a research which deserves to be extended and depend both from the mathematical viewpoint and in connection with possible applications to natural, social and economic behaviour of large collectivities in relation to external phenomena. It is worth stressing that the dependence of the forcing term on the distribution $\mathbf{f}$ turns out to be of the greatest relevance in order to give a description of the evolution of the system that could fit the real conditions of the phenomena. And, of course, different forms of such dependence require to be analysed in future research aiming to continue this one, along with deeper and more technical properties, like stability and periodicity of solutions and possible bifurcation effects, that however allow us to expect deeper insights in all possible applications to natural phenomena and to social and economic behaviours of our present complicated society, [48–52]. Furthermore, it may be of interest the research and study of analytical conditions towards the shape of the external force field such that the existence of a bounded and positive solution, globally in time, is ensured. We expect some dependency of results on both initial data and parameters of the system, i.e. interactions rates, transition probabilities and birth/death rates. Similarly, the study of case with a continuous activity variable $u$ is of interest. This latter situation will provide a kinetic model defined by a system of nonlinear integro-differential equations, more suitable for some specific applications.