\subsection{\@startsection                        %
  {subsection}{2}{0mm
  {1.5\baselineskip
  {1\baselineskip
  {\indent\normalfont\normalsize\itshape}
\renewcommand*{\@seccntformat}[1]
  \normalfont\bfseries\csname the#1\endcsname\;-               %
}                                                              %
\renewcommand\subsubsection{\@startsection                     %
  {subsubsection}{2}{0mm
  {1.5\baselineskip
  {1\baselineskip
  {\indent\normalfont\normalsize\texttt}
\makeatother                                                   %
                                                               %

\begin{document}
\thispagestyle{empty}


$ $ \vspace {-2.2cm}
\begin{center}
\rule{8.8cm}{0.5pt}
\\[-5pt]
\rule{8.8cm}{0.5pt}
\end{center}
\vspace {1.7cm}


\begin{center}
{\sc\large Giuseppe Toscani},  \
{\sc\large Mattia Zanella}   \
\end{center}
\vspace {1.1cm}

\centerline{\large{\textbf{On a kinetic description of Lotka-Volterra dynamics}}}


\renewcommand{\thefootnote}{\fnsymbol{footnote}}


\renewcommand{\thefootnote}{\arabic{footnote}}
\setcounter{footnote}{0}

\vspace{0.6cm}
\begin{center}
\begin{minipage}[t]{11cm}
\small{
\noindent \textbf{Abstract.}
Owing to the analogies between the problem of wealth redistribution with taxation in a multi-agent society, we introduce and discuss a kinetic model describing  the statistical distributions in time of the sizes of groups of biological systems with prey-predator dynamic. While the evolution of the mean values is shown to be driven by a classical Lotka-Volterra system of differential equations, it is shown that the time evolution of the probability distributions of the size of groups of the two interacting species is heavily dependent both on a kinetic redistribution operator and the degree of randomness present in the system. Numerical experiments are given to clarify the time-behavior of the distributions of groups of the species.
\medskip

\noindent \textbf{Keywords.}
Kinetic models, Boltzmann equation, Wealth distribution, Taxation and redistribution, Lotka-Volterra equations.
\medskip

\noindent \textbf{Mathematics~Subject~Classification:}
91B60, 82C40, 35B40.

}
\end{minipage}
\end{center}

\bigskip


\section{Introduction}

In one of the scientific collaborations between one of the present authors with Professor Giampiero Spiga in the area of kinetic theory of multi-agent systems, the problem of taxation and consequent redistribution of taxed resources in a western society had been studied  \cite{BST}. The purpose of the work was the possibility of understanding how, in the case of a fixed mean wealth, redistribution operators can affect the tails of the equilibrium density of a kinetic model for wealth exchanges.

In this area, there have been several attempts to model simple economies in which a wealth exchange process produces in time a steady profile of wealth distribution similar to that observed in real economies. While the most popular approach  relies on methods borrowed from statistical mechanics for interacting particle systems \cite{Cha,ChCh,CCM,CCS,CPT,DY,DMT,DMT2,Hay,IKR,MT,TTZ,Sla} (cf. also the books \cite{CYC,PT}), other modelling settings that might be considered to constitute a mesoscopic approach is based on generalized Lotka-Volterra models \cite{MBR,SR}. These dynamical systems, also known as the predator-prey equations,  are frequently used to describe the evolution of biological systems in which two species interact, one as a predator and the other as prey \cite{Murray}. In particular, denoting by $x(t)\ge 0$ and $y(t)\ge0$  the size of the populations at time $t\ge 0$, their time evolution is supposed to be influenced by predation, growth and death processes. Therefore, this dynamics is classically described by the following pair of ordinary differential equations
\begin{equation} \label{LV}
\begin{cases}
\frac{d}{dt}x(t) = \alpha x(t) -\beta \,x(t)y(t) ,\\
 \frac{d}{dt} y(t)= - \delta\, y(t) + \gamma \,x(t)y(t) ,
 \end{cases}
\end{equation}
where the preys' population is characterized by the rate of Malthusian growth $\alpha\ge 0$ and by a predation rate $\beta\ge 0$. Furthermore, the predators' population is subject to losses due to death  with rate  $\delta\ge 0$ and on the growth due to predation defined by the parameter $\gamma\ge0$.

While less popular, the approach proposed in \cite{MBR,SR} is capable to establish an interesting parallel between prey-predator dynamics and economic modelling. This analogy, as documented by \cite{Gan}, goes back to the Italian economist Giuseppe Palomba, who first had used  Lotka-Volterra equations \fer{LV} in the 1939 book \cite{Pal}, and to Goodwin, who considered applications of Lotka-Volterra system to economy in the 1967 book \cite{Good}.

In this paper, we provide a kinetic interpretation of the introduced dynamics by means of a system of Boltzmann-type equations. The interacting species are characterized from a statistical perspective according to their size. Hence, the temporal evolution of the statistical distributions of two species, classically denoted as prey and predators,  is subject to variation due to birth and death as well as to interactions between the two families (predation). The building blocks of the model are given by the redistribution operator introduced  in \cite{BST}, which allows for a general representation of the birth and death process of the involved species, and a kinetic bilinear operator of Boltzmann type which quantifies the interaction between species. 
Interestingly, the structure of the interaction is similar to that recently introduced by the authors in connection with the study of kinetic models suitable for describing several patterns in compartmental epidemiology connected with the transition between susceptible and infected individuals \cite{DPaTZ, DPTZ}.

The linear operator introduced in \cite{BST} is sufficiently flexible to redistribute to agents a constant amount of wealth independently of the wealth itself, or to redistribute proportionally (or inversely proportionally) to their wealth. Hence, in the new application to a biological system, it allows the number of births  to be quantified in various ways and in terms of the size of the groups of the species.

In more details, the paper is organized as follows: in Section \ref{sect:2} we introduce a kinetic model describing the evolution of the sizes of groups of biological systems obeying a prey-predator dynamic. The transitions depend on elementary interactions determining the evolution in time of probability densities described by a system of non-Maxwellian Boltzmann-like equations. A redistribution operator depending only on the first two moments of the distribution is then considered. In Section \ref{sec:3} we recall basic results on the grazing limit in kinetic theory and we provide a Fokker-Planck description of the dynamics.  Finally, in Section \ref{sec:4} we provide several numerical tests to verify the consistency of the approach with the classical Lotka-Volterra model and we observe the emergence of fat-tailed distributions, in the group of preys, for suitable choices of the redistribution operator. 

\section{The model}\label{sect:2}

Following the general approach of \cite{PT} we consider a multiagent system composed by two species, identified as preys and predators,  that are characterized by their number of individuals $x\ge 0$ (respectively $y\ge0$). We denote by $f(x,t)$ (respectively $g(y,t)$) be the statistical distributions at a certain time $t \ge 0$ of the group of preys (respectively predators), classified in terms of their sizes $x,y \in \mathbb R_+$. 
Therefore, $f(x,t)dx$ indicates the fraction of the population preys with size in $[x,x + dx]$ at time $t\ge0$. Similarly, $g(y,t)dy$ is the fraction of the populations of predators characterized by a size in $[y,y+dy]$ at time $t\ge0$. Moreover, without loss of generality we fix the initial masses of the distributions $f$ and $g$ equal to one. Therefore, the mean sizes of the populations are given by 
\[
m_{f}(t) = \int_{\mathbb R_+}xf(x,t)dx,\qquad m_{g}(t) = \int_{\mathbb R_+}yg(y,t)dy.
\]
The goal is to characterize the time evolution of the pair of probability densities $f(x,t)$, $g(y,t)$ subject to changes including  births and deaths events, and obeying to objective mechanical and behavioral principles when interacting each other.

Let us first characterize the interaction rules between the two species. To achieve this, let us fix a certain time $t= t_0$, considered as the initial time of the observations,  and let $x_0,y_0 $ be the numbers of preys and predators in two groups at time $ t_0 $. Then, we evaluate the statistical distribution of the size of the groups at regular time intervals $ t_{n+1} $, with $n$ nonnegative integer, $ n\in \mathbb{N}_+$, consequent to  the elementary variations $x_n \to x_{n+1}$ and  $y_n \to y_{n+1}$ of the groups in each time interval. Since a prey-predator  interaction leads to a decrease in the number of preys, and, in the presence of a sufficiently high number of preys, to an increase of the number of predators, this change can be modeled as follows
\begin{equations}\label{coll1}
 x_{n+1} &= x_n  -  \Phi_\epsilon(y_n) x_n +   \eta_\epsilon(y_n) x_n. \\
 y_{n+1} &= y_n  + \Psi_\epsilon(x_n) y_n +  \tilde\eta_\epsilon(x_n)  y_n.
\end{equations}
Thus, at each time step, the number  of preys and predators is modified by two different mechanisms, expressed in mathematical terms by two multiplicative terms, both parameterized by a small positive parameter $\epsilon \ll 1$, measuring the intensities of the single variations.

To clarify, the  function $\Phi_\epsilon(\cdot)$ in the first equation of \fer{coll1} quantifies the predation rate, and is consistently assumed to be increasing with respect to the size of the predator population. Likewise, the function $\Psi_\epsilon(\cdot)$ in the second equation of \fer{coll1} quantifies the growth of the predators, increasing with respect to the prey population at least when the number of preys is sufficiently consistent.  In reason of the previous remarks, in what follows we  fix 
\begin{equation}\label{Phi}
\Phi_\epsilon(y) =\epsilon\beta \frac y{1+y}; \quad 
\Psi_\epsilon(x) = \epsilon \gamma \frac {x-\mu}{1+x}.
\end{equation}
In \fer{Phi}, $\beta$, $\gamma$ and $\mu$ are positive constants. The constants $\beta$ and $\gamma$ measure the maximal predation rate, and, respectively, the maximal growth rate of predators. Last, the constant $\mu \ge 1$ measures the lowest size of the group of preys which is compatible with a growth of the group of predators. These constants are required to satisfy the bounds  
\begin{equation}\label{bb}
\epsilon\beta < 1; \quad \epsilon\gamma\mu < 1.
\end{equation}
These conditions ensure that the predation rate is always smaller than one, and that even in absence of preys, the number of predators does not collapse to zero in a single interaction.

In \fer{coll1}, random fluctuations of the groups due to unknown factors  in the environment, like immigration, natural events,  epidemics, etc.  are  expressed by  the independent random variables $\eta_\epsilon(x), \tilde\eta_\epsilon(x)$,  of zero mean and bounded variance, given by 
\begin{equation}\label{var1}
\langle \eta_\epsilon^2(y) \rangle  = \epsilon\sigma\frac y{1+y}; \quad  \langle \tilde\eta_\epsilon^2 (x)\rangle =    \epsilon\tilde\sigma \frac {x}{1+x}. 
\end{equation}
Note that, consistently with the deterministic variations, the variance of the random fluctuations of one class  is an increasing function of the size of the  other class.  

Starting from the elementary interactions \fer{coll1} we can easily write the kinetic equations describing the evolution in time of the pair of population densities. 
In a suitable scaling, the evolution of the densities  is
quantitatively described by a system of non-Maxwellian Boltzmann-like equations \cite{MTZ}
\begin{equation}
  \label{kine}
  \begin{split}
  \frac{\partial f}{\partial t} = Q_1(f,g), \\  \frac{\partial g}{\partial t} = Q_2(g,f),
  \end{split}
\end{equation}
where the bilinear operators on the right-hand side of both equations quantify the variations of the population densities due to the elementary interactions of type \fer{coll1} where the instantaneous microscopic  transitions $x\to x^\prime$ and $y\to y^\prime$ are defined as follows 
\[
\begin{split}
x^\prime &= x - \Phi_\epsilon(y)x + x \eta_\epsilon(y) \\
y^\prime &= y + \Psi_\epsilon(x)y+y\eta_\epsilon(x).
\end{split}
\]
The homogeneous Boltzmann equations \eqref{kine} can be fruitfully written in weak form. It corresponds to say that the
solution to \fer{kine} satisfies, for all smooth functions $\varphi$  the integro-differential equations
 \begin{equations}
  \label{kine-w}
\frac{d}{dt}\int_{\mathbb{R}_+}\varphi(x)f(x,t)\,dx &=
 \int_{\mathbb{R}_+^2} \kappa_1^\epsilon(x,y)  \left\langle \varphi(x^\prime)-\varphi(x)\right\rangle
f(x,t) g(y,t) \,dx \, dy, \\ 
\frac{d}{dt}\int_{\mathbb{R}_+} \varphi(y)g(y,t)\,dy &=
 \int_{\mathbb{R}_+^2} \kappa_2^\epsilon(x,y)  \left\langle \varphi(y^\prime)-\varphi(y)\right\rangle
f(x,t) g(y,t) \,dx \, dy.
 \end{equations} 
 In \fer{kine-w} the expression $\langle\cdot\rangle$ denotes mathematical expectation with respect to the random variables $\eta_\epsilon$, $\tilde \eta_\epsilon$, and the functions $\kappa_i^\epsilon(x,y)$, $i =1,2$ are the interaction frequencies. In what follows, to maintain the connection between the kinetic model and the Lotka--Volterra system, we fix these frequencies in the form
 \begin{equation}\label{fre}
 \kappa_1^\epsilon(x,y) =\kappa^\epsilon(y)= \frac{1+y}\epsilon; \quad \kappa_2^\epsilon(x,y) = \kappa^\epsilon(x) = \frac{1+x}\epsilon.
 \end{equation}
 Within this choice, the frequency of interactions in the first equation grows with the number of predators, while in the second equation grows with the number of preys.
 
 By choosing $\varphi(x) =1$ in \fer{kine-w} we can easily observe that  $f(x,t)$ and $g(y,t)$ {remain probability densities} for any  time $t \ge0$, if they are so initially
 \begin{equation}\label{m0}
\int_{\mathbb{R}_+} f(x, t)\,dx = \int_{\mathbb{R}_+} g(y,t)\,dy = 1 .
 \end{equation}


Then, by choosing $ \varphi(x) = x$ in both equations \fer{kine-w} one easily obtains that the evolution of the mean values is defined by the coupled system
\begin{equation}\label{mean2}
\begin{split}
\frac{dm_f(t)}{dt} &= -\beta \,m_f(t)m_g(t); \\
 \frac{dm_g(t)}{dt} &= - \gamma\mu\, m_g(t) + \gamma \,m_f(t)m_g(t).
\end{split}
\end{equation}
Thus, if the functions $\Phi_\epsilon$ and $\Psi_\epsilon$ in \fer{coll1} are given by \fer{Phi}, and the collision frequencies in \fer{kine-w} are equal to \fer{fre}, the mean values of the densities satisfy the Lotka--Volterra system \fer{LV} in which the rate of birth $\alpha$ of the preys is set equal to zero.

We may obtain a kinetic interpretation of the birth rate by complementing the kinetic system \fer{kine} with an additional operator  which accounts for the variation of densities due to births of individuals in the groups. As already discussed in the introduction, we will resort to the redistribution operator introduced in \cite{BST}. In this paper, wealth redistribution mimics the action of a government that use the proceeds of taxation to improve the financial situation of agents based on their income.

According to \cite{BST}, for a given probability density function $h(x,t)$, with mean $m(t)$, we  assume the birth operator of the form
 \begin{equation}\label{redis}
 R_\chi^\alpha(h)(x,t) = \alpha \frac\partial{\partial x} \left[ \left(\chi x - (\chi + 1)m(t)\right)h(x,t)\right].
  \end{equation}
 The weight factor multiplying the distribution function inside the square brackets in \fer{redis}  involves in the mechanism the moments of order zero and one. Such a weight function adds one real parameter $\chi$ to the dynamics. This new parameter characterizes the type of redistribution.  The parameter $\alpha>0$ is determined by the constraint that, while the redistribution operator preserves the total mass, the mean number of individuals increases by a factor $\alpha$. In fact, for each value of the constant $\chi\in \mathbb R$, we have 
 \begin{equation}\label{res1}
  \int_{\mathbb{R}_+}x R_\chi^\alpha(h)(x,t) \, dx = \alpha\, m(t),
 \end{equation}
and, provided $h(x,t)$ satisfies in addition the condition $h(0,t) = 0$, we also get
 \begin{equation}\label{res0}
  \int_{\mathbb{R}_+}R_\chi^\alpha(h)(x,t) \, dx = 0.
 \end{equation}
We can observe that the condition $h(0,t) = 0$ is not necessary in the special case $\chi = -1$, which reproduces an
anti-drift operator introduced by Slanina in  \cite{Sla} to mimic wealth redistribution effects. 

As shown in \cite{BST}, the effects of the operator $R_\chi^\alpha$ on the birth process on a group described by the density function $f$ are closely related to the values of the parameter $\chi$. The case in which the percentage of births is equally distributed independently of the size of the group  is achieved with the special option $\chi=0$. 
In all other cases the births are selective, and may correspond to some partition strategy.   {In the present case, as far as the preys are concerned,  it is suitable to take into account positive values of the parameter $\chi$, which correspond to impose that the percentage  of births is inversely proportional to the size of the group.  On the contrary, we assume $\chi=0$ for the predators group.}

Taking into account the previous discussion, the full system is obtained by adding to equations \fer{kine-w} the birth operators defined in \fer{redis}. Therefore, we have
 \begin{equations}
  \label{kine-wd}
\frac{d}{dt}\int_{\mathbb{R}_+}\varphi(x)f(x,t)\,dx = & \int_{\mathbb{R}_+}\varphi(x)R_\chi^\alpha(f)(x,t)\,dx \, +\\
+ &  \int_{\mathbb{R}_+^2} \kappa_1^\epsilon(x,y) \left\langle \varphi(x^\prime)-\varphi(x)\right\rangle
f(x,t) g(y,t) \,dx \, dy, \\ 
\frac{d}{dt}\int_{\mathbb{R}_+}\varphi(y)g(y,t)\,dy =& \int_{\mathbb{R}_+}\varphi(y)R_{0}^\nu(g)(y,t)\,dy\, + \\
+ &  \int_{\mathbb{R}_+^2} \kappa_2^\epsilon(x,y) \left\langle \varphi(y^\prime)-\varphi(y)\right\rangle
f(x,t) g(y,t) \,dx \, dy.
 \end{equations} 
In \fer{kine-wd}, the constant $\nu$ which defines the percentage of the mean number of births of the group of predators, is assumed to satisfy the bound $ \nu < \gamma\mu$. In this case, resorting to property  \fer{res1}, it is immediate to conclude that the evolution of the mean numbers $m_f(t)$ and $m_g(t)$ satisfies the Lotka-Volterra system \fer{LV}, where now $\delta = \gamma\mu -\nu >0$.

\section{The limit of grazing predation rate}\label{sec:3}

We investigate now the situation in which most of the interactions between species correspond to a very small variation of the sizes of the groups, but at the same time it is possible to keep trace at the macroscopic level of all phenomena affecting predation rules. This kind of asymptotic analysis is referred to as \emph{quasi-invariant limit}, and it allows to prove that under proper assumptions on the interaction parameters,  the kinetic Boltzmann system may be approximated by a system of Fokker-Planck equations, allowing in various cases an analytical investigation of the shape of the solution. This limit corresponds to consider in \fer{kine-w} a Taylor expansion up to the second order of the differences 
\[
\varphi(x^\prime)-\varphi(x), \quad \varphi(y^\prime)-\varphi(y),
\]
and subsequently to consider $\epsilon \to 0$, see \cite{PT,T} for further details. It can be proven that, in this asymptotic regime, the kinetic system \fer{kine-w} is well-approximated by the following system of  Fokker-Planck equations
 \begin{equations}
  \label{FP}
\frac{\partial f(x,t)}{\partial t} &= \frac{\partial}{\partial x}\Bigg\{ \frac{\sigma m_g(t)}2 \frac{\partial(x^2 f(x,t))}{\partial x}+ \\
 & +  \left[ (\beta m_g(t) +\alpha\chi) x - \alpha(\chi +1)m_f(t)\right]f(x,t)\Bigg\} \\
\frac{\partial g(y,t)}{\partial t} &= \frac{\partial}{\partial y}\Bigg\{\frac{\tilde\sigma m_f(t)}2 \frac{\partial(y^2 g(y,t))}{\partial y} + \\ & +  \left[ \gamma(\mu -m_f(t))y -\nu m_g(t)\right]g(y,t)\Bigg\}.
 \end{equations}
complemented by no-flux boundary conditions \cite{FPTT}. 
We can observe that the evolution of the mean values of \fer{FP} is consistent with the one obtained at the Boltzmann level  in \fer{kine-w}, so that the mean values $m_f(t)$, $m_g(t)$ satisfy the Lotka--Volterra system \fer{LV}.

\section{Numerical tests}\label{sec:4}

In this section we provide several numerical tests to show the consistency of the proposed approach with the classical Lotka-Volterra dynamical system. Furthermore, we will show numerically how the kinetic approach allows additional information to be obtained on the large time behaviour of the densities of preys and predators in dependence of the parameters in the birth operators $R_\chi^\alpha(\cdot)$, $R_0^\nu(\cdot)$. For the approximation of the system of Fokker-Planck equations in \fer{FP} we will adopt the structure preserving numerical strategy developed in \cite{PZ}. These methods are capable to reproduce large time statistical properties of the solution density with arbitrary accuracy together with the preservation of its main physical properties, like its positivity.  We complement  the model \fer{FP}  with no-flux boundary conditions.

In all the following tests we will fix the parameter as reported in Table \ref{tab:1}.

\begin{table}
\begin{center}
\begin{tabular}{ |c| c| c| }
\hline
Parameter & Value & Meaning \\
\hline\hline  
 $\alpha$ & 1.00 & Growth rate preys \\ 
 $\beta$ & 0.50 & Predation rate preys \\  
 $\mu$ & 10.00 & Lowest size preys   \\
 $\gamma$ &0.15 & Growth rate predators \\
  $\nu$ & 1.00 & Birth rate predators\\
 $\delta$ & $\gamma\mu-\nu>0$ & Death rate predators\\
 $\sigma^2$ & $10^{-3}$ & Diffusion coefficient preys \\
 $\tilde \sigma^2$ & $10^{-3}$ & Diffusion coefficient predators \\
 \hline 
\end{tabular}
 \caption{Parameters chosen in the numerical tests. }
 \label{tab:1}
 \end{center}
\end{table}

\subsection{Test 1. Evolution of the mean values}

We consider the Fokker-Planck model \fer{FP} with initial distribution 
\begin{equation}\label{eq:fg0}
\begin{split}
f(x,0) = 
\begin{cases}
C_f \exp\left (-\dfrac{(x-m_f^0)^2}{2\kappa_f}\right) & x \ge 0 \\
0 & x<0
\end{cases} \\
g(y,0) = 
\begin{cases}
C_g \exp\left (-\dfrac{(y-m_g^0)^2}{2\kappa_g}\right) & y \ge 0 \\
0 & y<0,
\end{cases} 
\end{split}
\end{equation}
with $C_f,C_g>0$ normalization constants and with $m_f = 4$, $m_g = 3$ and $\kappa_f = \kappa_g = \frac{1}{10}$.

In Figure \ref{fig:evolution} we show the evolution of the mean values $m_f$, $m_g$ of \fer{FP} over the time horizon $[0,15]$ for an increasing number of gridpoints $N = 201,401,801$ for the discretization of the interval $[0,L]$, $L= 10$. The choice of parameters has been defined in Table \ref{tab:1} and we further considered $\chi = 0$. We highlight in red the numerical solution of the Lotka-Volterra model obtained in the same regime of parameters. We can observe how, consistently with the analytical considerations in Section \ref{sec:3}, we can obtain numerically a good approximation of the original dynamical system in terms of the mean values of the distributions of preys $f(x,t)$ and predators $g(y,t)$. 

\begin{figure}
\centering
\includegraphics[scale = 0.325]{figure/evo_mf}
\includegraphics[scale = 0.325]{figure/evo_mg}
\caption{\textbf{Test 1}. Comparison of the evolution of the mean values for the populations of preys $f(x,t)$ and predators $g(y,t)$ with the Lotka-Volterra dynamical system \fer{LV}. The Fokker-Planck system has been approximated by means of a 4th order SP methods with RK4 time integration and an increasing number of gridpoints $N = 201,401,801$ over the domain $[0,10]$. We considered $\Delta t = \Delta x^2= \Delta y^2$. }
\label{fig:evolution}
\end{figure}


\subsection{Test 2. Emergence of fat-tailed distributions}

In this test we show the evolution of the approximated distributions $f(x,t)$ and $g(y,t)$ solution to \fer{FP} for varying values of the parameters $\chi \in \mathbb R$. All the other parameters have been chosen in agreement with Table \ref{tab:1}. We consider the initial distributions defined in \fer{eq:fg0} and we solve the system of Fokker-Planck equations \fer{FP} with $N = 801$ gridpoints in the interval $[0,10]$, and time frame $[0,15]$ with $\Delta t = \Delta x^2$. 

In Figure \ref{fig:surf} we provide the total evolution of the distribution functions in the cases $\chi = -1$ (top row) and $\chi = 1$ (bottom row). As specified in Section \ref{sect:2} the two choices are coherent with two different redistribution policies. Indeed, for all $\chi >0$ the number of new births is distributed in small values of $x\ge 0$, whereas for $\chi<0$ the new births are concentrated in high values of $x \ge 0$. Coherently with \cite{BST}, we may observe how the preys population may exhibit a rather different behaviour for large times. Indeed, fat-tail-type distributions appears in the case $\chi = -1$. 

To clarify this behaviour, we present in Figure \ref{fig:3times} a plot of the distributions of preys and predators at three fixed times $t = 0,1,8$. We can observe how for $\chi= -1$ the distribution of preys develops fat tails in finite time. This behaviour is further observable in loglog scale, see Figure \ref{fig:loglog}. 

\begin{figure}
\centering
\includegraphics[scale = 0.325]{figure/f_chim1}
\includegraphics[scale = 0.325]{figure/g_chim1} \\
\includegraphics[scale = 0.325]{figure/f_chip1}
\includegraphics[scale = 0.325]{figure/g_chip1} 
\caption{\textbf{Test 2}. Evolution of the distributions $f(x,t)$ and $g(y,t)$ over the time interval $[0,15]$ and for $x,y \in [0,10]$. Approximation of the system of Fokker-Planck equations obtained with 4th order SP scheme with RK4 time integration, $N = 801$ gridpoints and $\Delta t = \Delta x^2 = \Delta y^2$. We considered $\chi = -1$ (top row) or $\chi = 1$ (bottom row) in \fer{FP}.}
\label{fig:surf}
\end{figure}

\begin{figure}
\includegraphics[scale = 0.325]{figure/f_time_chim1}
\includegraphics[scale = 0.325]{figure/g_time_chim1}\\
\includegraphics[scale = 0.325]{figure/f_time_chi1}
\includegraphics[scale = 0.325]{figure/g_time_chi1}
\caption{\textbf{Test 2}. Distributions of preys $f(x,t)$ (left) and predators $g(y,t)$ (right) at three different time steps $t = 0,1,8$, $x,y \in [0,10]$. We considered $\chi = -1$ (top row) or $\chi = 1$ (bottom row) in \fer{FP}.}
\label{fig:3times}
\end{figure}


\begin{figure}
\centering
\includegraphics[scale = 0.325]{figure/semilog_fg_chim1}
\includegraphics[scale = 0.325]{figure/semilog_fg_chip1}
\caption{\textbf{Test 2}. Distributions of preys $f(x,t)$ and predators $g(y,t)$ at time $T = 15$ in the cases $\chi = -1 $ (left) or $\chi = 1$ (right). }
\label{fig:loglog}
\end{figure}

\section{Conclusions}

Motivated by the similarities between the problem of wealth redistribution with taxation in a multi-agent society, we introduced in this paper a kinetic description of  the statistical distributions in time of the sizes of groups of biological systems obeying  a prey-predator dynamic. The kinetic system is built in such a way that the evolution of the mean values is driven by a classical Lotka-Volterra system of differential equations. Among other characteristics, it is shown that the time evolution of the probability distributions of the size of groups of the two interacting species is heavily dependent on the  kinetic redistribution operator quantifying the births. Numerical experiments clarify the time-behavior of the distributions of groups of the species, and highlight the possibility that the distribution of the sizes of preys could exhibit fat tails.
 
 \vskip 2cm
\vspace{0.5cm} \indent {\it
A\,c\,k\,n\,o\,w\,l\,e\,d\,g\,m\,e\,n\,t\,s.\;} This work has been written within the activities of GNFM group of INdAM (National Institute of
High Mathematics). GT wish to acknowledge partial support by IMATI, Institute for Applied Mathematics and Information Technologies ``Enrico Magenes'', Via Ferrata 5 Pavia, Italy. MZ acknowledges the support of MUR-PRIN2020 Project No.2020JLWP23 (Integrated Mathematical Approaches to Socio–Epidemiological Dynamics).

\bigskip
\begin{center}

