Dynamics of a stochastic predator-prey model
with fear effect and hunting cooperation

H. Qi and X. Meng,
Dynamics of a stochastic predator-prey model,
with fear effect and hunting cooperation,
in Journal of Applied Mathematics and Computing, vol. 69, 2023

Exam Presentation
Stochastic Modelling and Simulation Course
SDIC Master Degree, University of Trieste (UniTS)

Marco Tallone
September 2025

Table of Contents

1. Model Description
  • Variables and Parameters
  • Deterministic model
  • Stochastic Model
2. Simulation methods
  • Euler-Maruyama
  • Milstein
  • Runge-Kutta
3. Results and
     Model Behavior
  • Asymptotic behavior
  • Equilibrium points
  • Noise-induced effects

Model Description

Deterministic Model

To understand the deterministic model, let's start from a classic Lotka-Volterra model:

\[ \begin{cases} \frac{du}{dt} = r_0 u - a u v\\ \frac{dv}{dt} = a u v - m v\ \end{cases} \]

with:
- $u(t)$ prey population density
- $v(t)$ predator population density
- $r_0$ prey birth rate
- $a$ predator attack rate
- $m$ predator death rate

Deterministic Model

As prey density $u$ rises, individuals compete for limited food/space, thus we add an intraspecific competition term introducing a prey competition strength parameter $\beta$:

\[ \begin{cases} \frac{du}{dt} = r_0 u - \beta u^2 - a u v\\ \frac{dv}{dt} = a u v - m v\ \end{cases} \]

Deterministic Model

Predators need time to chase, hunt and eat preys. The number of preys per unit time a predator can hunt ($h$) is limited and described by a Holling type II functional response: \[ \text{Holling II:}\quad \frac{a u}{1 + hau} \] thus the model becomes:

\[ \begin{cases} \frac{du}{dt} = r_0 u - \beta u^2 - \frac{a u v}{1 + hau}\\ \frac{dv}{dt} = \frac{a u v}{1 + hau} - m v\ \end{cases} \]

Deterministic Model

In this model, predators can cooperate with rate $\alpha$ to improve their baseline attack rate $p$. Hence the overall attack rate is not constant: \[ a(v) = p + \alpha v \] which results in:

\[ \begin{cases} \frac{du}{dt} = r_0 u \mathbf{- \beta u^2} - \frac{(p + \alpha v) u v}{1 + h (p + \alpha v) u}\\ \frac{dv}{dt} = \frac{(p + \alpha v) u v}{1 + h (p + \alpha v) u} - m v\ \end{cases} \]

Deterministic Model

Moreover, cooperation affect the prey birth rate in the equation. This is further influenced by the fear $k$ that preys perceive from predators: \[ r_0 \rightarrow \frac{r_0}{1 + k \alpha v} \] Hence:

\[ \begin{cases} \frac{du}{dt} = \frac{r_0 u}{1 + k \alpha v} - \beta u^2 - \frac{(p + \alpha v) u v}{1 + h (p + \alpha v) u}\\ \frac{dv}{dt} = \frac{(p + \alpha v) u v}{1 + h (p + \alpha v) u} - m v\ \end{cases} \]

Deterministic Model

Finally we introduce a conversion efficiency $\theta$ which models how efficiently the energy from consumed prey is converted into new predator births or growth.

\[ \begin{cases} \frac{du}{dt} = \frac{r_0 u}{1 + k \alpha v} - \beta u^2 - \frac{(p + \alpha v) u v}{1 + h (p + \alpha v) u}\\ \frac{dv}{dt} = \frac{\theta (p + \alpha v) u v}{1 + h (p + \alpha v) u} - m v\ \end{cases} \]

Deterministic Model

The end result is a prey-predator model with fear effect and hunting cooperation:

\[ \begin{cases} \frac{du}{dt} = \frac{r_0 u}{1 + k \alpha v} - \beta u^2 - \frac{(p + \alpha v) u v}{1 + h (p + \alpha v) u}\\ \frac{dv}{dt} = \frac{\theta (p + \alpha v) u v}{1 + h (p + \alpha v) u} - m v \end{cases} \tag{1} \]

The state space is therefore define as: $\quad \mathcal{S} = \{x = (u, v) \in \mathbb{R}^2 | u(t) \geq 0 \land v(t) \geq 0, \forall t \geq 0\}$

Parameters of the model

here we summarize the parameters of the model introduce so far:

Parameter Description What it models
$r_0$ Prey birth rate Reproduction rate of preys in absence of predators
$\alpha$ Predator hunting cooperation rate How much predators cooperate to hunt preys
$\beta$ Competition strenght betwen preys How much preys compete among themselves for limited resources
$\theta$ Conversion efficiency of the predator How efficiently the energy from consumed prey is converted into new predator births or growth
$k$ Fear level induced by predator How much the presence of predators affects the prey birth rate
$h$ Predator handling time of prey Time needed by a predator to chase, hunt and eat a prey
$m$ Predator death rate Natural death rate of predators in absence of preys
$p$ Predator baseline attack rate Attack rate of a non-cooperating predator

Stochastic Model

Populations subject to environmental disturbances $\Rightarrow$ noise indispensable modelling factor
Hence, a stochastic prey-predator model with fear effect and hunting cooperation factor is proposed:

\[ \begin{cases} \frac{du}{dt} = \left[\frac{r_0 u}{1 + k \alpha v} - \beta u^2 - \frac{(p + \alpha v) u v}{1 + h (p + \alpha v) u}\right]dt + u(\sigma_{11} + \sigma_{12} u) dW_1(t)\\ \frac{dv}{dt} = \left[\frac{\theta (p + \alpha v) u v}{1 + h (p + \alpha v) u} - m v\right]dt + v(\sigma_{21} + \sigma_{22} v) dW_2(t) \end{cases} \tag{2} \]

where $\sigma_{ij} > 0$ are the intensity weights of the white noise and $dW_1(t), dW_2(t)$ are mutually independent standard Brownian motions.

Drift and diffusion functions

We shall also formalize the vectorial drift function $f: \mathbb{R}^2 \to \mathbb{R}^2$ and
diffusion function $g: \mathbb{R}^2 \to \mathbb{R}^2$ associated with the stochastic model:

\[ f(u, v) = \begin{bmatrix} \frac{r_0 u}{1 + k \alpha v} - \beta u^2 - \frac{(p + \alpha v) u v}{1 + h (p + \alpha v) u}\\ \frac{\theta (p + \alpha v) u v}{1 + h (p + \alpha v) u} - m v \end{bmatrix} \quad\text{and}\quad g(u, v) = \begin{bmatrix} u(\sigma_{11} + \sigma_{12} u)\\ v(\sigma_{21} + \sigma_{22} v) \end{bmatrix} \]

as well as the first derivatives of the components of the latter as:

\[ g'(u, v) = \begin{bmatrix} \frac{dg_u(u,v)}{du}\\ \frac{dg_v(u,v)}{dv}\\ \end{bmatrix} = \begin{bmatrix} \sigma_{11} + 2 \sigma_{12} u\\ \sigma_{21} + 2 \sigma_{22} v \end{bmatrix} \]

Numerical Simulation Methods

Numerical Simulation Methods

$3$ numerical simulation methods have been implemented to simulate the stochastic model:

Euler-Maruyama Algorithm
Milstein Algorithm
(used in the paper)
Stochastic Runge-Kutta Scheme

In all of them, we're approximating the real solution $x(t) = (u(t), v(t))$ on a temporal interval $[0, T]$ by partitioning it in $N$ discrete time-steps and iteratively computing the discretized model solution $x_k = (u_k, v_k)$ using an iterative update rule.

Euler-Maruyama Algorithm

The iterative update rule for the Euler-Maruyama Algorithm is:

\[ x_{k+1} = x_k + f(x_k) \Delta t + g(x_k) \sqrt{\Delta t} \xi_k \]

where $\Delta t = T / N$ is the uniform length of the single time-step and
$\xi_k$ for $k = 1, 2, \dots$ are independent Gaussian random variables $\mathcal{N}(0, 1)$.

Milstein Algorithm

The iterative update rule for the Milstein Algorithm is:

\[ x_{k+1} = x_k + f(x_k) \Delta t + g(x_k) \sqrt{\Delta t} \xi_k + \frac{1}{2} g(x_k) g'(x_k) (\zeta_k^2 - 1) \Delta t \]

where $\Delta t = T / N$ is the uniform length of the single time-step and
$\xi_k, \zeta_k$ for $k = 1, 2, \dots$ are independent Gaussian random variables $\mathcal{N}(0, 1)$.

Stochastic Runge-Kutta Scheme

The iterative update rule for the Stochastic Runge-Kutta Scheme is:

\[ x_{k+1} = x_k + f(x_k) \Delta t + g(x_k) \sqrt{\Delta t} \xi_k + \frac{1}{2} \left( g(x_k + f(x_k) \Delta t + g(x_k) \sqrt{\Delta t }) - g(x_k) \right) (\zeta_k^2 - 1) \sqrt{\Delta t} \]

where $\Delta t = T / N$ is the uniform length of the single time-step and
$\xi_k$ for $k = 1, 2, \dots$ are independent Gaussian random variables $\mathcal{N}(0, 1)$.

Comparison between simulation methods

Algorithms comparison

Convergence of simulation methods

$ m_{Euler} \approx 0.6 \quad\quad m_{Milstein} \approx 0.9 \quad\quad m_{SRK} \approx 0.9 $

Algorithms convergence

Simulation Results & Model Behavior

Numerical Simulation Results

Parameters:
$r_0 = 2.0$
$p = 0.5$
$k = 1.5$
$h = 1.2$
$\beta = 0.5$
$\theta = 0.6$
$\alpha = 0.24$
$m = 0.34$

Noise intensities:
\[ \sigma = \begin{bmatrix} 10^{-2} & 10^{-2}\\ 10^{-2} & 10^{-2} \end{bmatrix} \]

Normal Simulation

Simulation close to deterministic equilibrium

Deterministic equilibrium:
\[ f(u^*, v^*) = (0, 0) \] \[ \Downarrow \] \[ (u^*, v^*) = (2.3075, 1.1143) \]

Equilibrium Simulation

Simulations from random IC

Multiple IC simulations

Asymptotic Behavior Theorems

The paper states two important theorems:

$\text{\textbf{THEOREM 3.1}}:$ Extinction of the prey population \[ \text{when} \quad R_1 = r_0 - \frac{1}{2}\sigma_{11}^2 < 0 \quad \Rightarrow \quad\lim_{t \to \infty} u(t) = 0 \] Hence, prey population $u(t) \to 0$ will be extinct.
$\text{\textbf{THEOREM 3.2}}:$ Extinction of the predator population \[ \text{when} \quad R_2 = \frac{\theta}{h} - \left(m + \frac{1}{2}\sigma_{21}^2\right) < 0 \quad \Rightarrow \quad\lim_{t \to \infty} v(t) = 0 \] Hence, predator population $v(t) \to 0$ will be extinct.

In the above simulation: $\quad R_1 \approx 2 \quad \text{and} \quad R_2 \approx 0.16 \quad$, indeed both populations survived.

Predator Extinction

Parameters:
\[ \begin{array}{ll} r_0 = 2.0 & \beta = 0.5\\ p = 0.5 & \theta = 0.6\\ k = 1.5 & \alpha = 0.24\\ h = 1.2 & m = 0.34 \end{array} \] Noise intensities:
\[ \sigma = \begin{bmatrix} 10^{-2} & 10^{-2}\\ \mathbf{0.6} & 10^{-2} \end{bmatrix} \] R-values:
\[ \begin{array}{ll} R_1 = 2.0 & R_2 = -0.02 \end{array} \]

Predator Extinction Simulation

Predator Extinction

⚠️ NOTE!
Notice that the extinction of the predator $\Rightarrow$ $v(t) \to 0$
As a consequence, prey population $u(t)$ grows as a simple linear logistic model: \[ \frac{du}{dt} = r_0 u - \beta u^2 = r_0 u \left(1 - \frac{u}{\frac{r_0}{\beta}}\right) \] with equilibria $u_{eq.} = 0$ and $u_{eq.} = \frac{r_0}{\beta}$.

With previous parameters $r_0 = 2.0$ and $\beta = 0.5 \Rightarrow$ $u(t) \to 4$ as $t \to \infty$.

Predator Extinction

Predator Extinction Animation

Prey Extinction

Parameters:
\[ \begin{array}{ll} r_0 = 0.17 & \beta = 0.5\\ p = 0.5 & \theta = 0.6\\ k = 1.5 & \alpha = 0.24\\ h = 1.2 & m = 0.34 \end{array} \] Noise intensities:
\[ \sigma = \begin{bmatrix} \mathbf{0.6} & 10^{-2}\\ 10^{-2} & 10^{-2} \end{bmatrix} \] R-values:
\[ \begin{array}{ll} R_1 = -0.01 & R_2 = 0.16 \end{array} \]

Prey Extinction Simulation

Steady State PDF in the Predator Extinction Case

Notice that as the predator goes extinct $v(t) \to 0$, the evolution of the prey $u(t)$ is approximated by: \[ \frac{du}{dt} = r_0 u - \beta u^2 + u(\sigma_{11} + \sigma_{12} u) dW \]

Hence, assuming $\sigma_{12} \ll 1$, it can be demonstrated that:

  • if $\sigma_{11}^2 > 2 r_0$ $\quad \Rightarrow \quad$ steady state $P(u) = \delta(u)$
    (in agreement with THEOREM 3.1 since it corresponds to $R_1 < 0$)
  • if $\sigma_{11}^2 < r_0$ $\quad \Rightarrow \quad$ steady state $P(u)$ is unimodal around equilibria $u_{eq.} = \frac{r_0}{\beta}$
    (in agreement with the previous simulation)
  • if $r_0 < \sigma_{11}^2 < 2 r_0$ $\quad \Rightarrow \quad$ steady state $P(u)$ is monotone decreasing
    (not presented in the paper)

Steady State PDF in the Predator Extinction Case

Case: $r_0 < \sigma_{11}^2 < 2 r_0$ with $\sigma_{11} = 0.77$ and $r_0 = 0.5$

Predator Extinction Animation 2

Hunting Cooperation $\;\alpha\;$ effect

Fear factor $\;k\;$ effect

Combined $\;\alpha\;$ and $\;k\;$ effect on Prey Population

Combined $\;\alpha\;$ and $\;k\;$ effect on Predator Population

White-Shot Noise Model

Extension: white shot noise $=$ acts as an impulsive force with a stochastic intensity withing some determined bounds $\Rightarrow$ series of randomly distributed impulses $\sim$ exponential distribution with rate $\lambda$, which amplitude is bounded in a given range.
In this case, the model is similarly described by euations $(2)$, with a small difference:

\[ \begin{cases} \frac{du}{dt} = \left[\frac{r_0 u}{1 + k \alpha v} - \beta u^2 - \frac{(p + \alpha v) u v}{1 + h (p + \alpha v) u}\right]dt + u(\sigma_{11} + \sigma_{12} u) \sum_i \xi_i \delta(t - t_i)\\ \frac{dv}{dt} = \left[\frac{\theta (p + \alpha v) u v}{1 + h (p + \alpha v) u} - m v\right]dt + v(\sigma_{21} + \sigma_{22} v) \sum_i \zeta_i \delta(t - t_i) \end{cases} \tag{3} \]

where $\delta(t-t_i)$ is the Dirac delta function indicating the presence of an impulsive event at time $t_i$ and $\xi_i, \zeta_i$ are the random amplitudes of the impulses uniformly sampled within a given range.

Simulation White-Shot Noise

Parameters:
\[ \begin{array}{ll} r_0 = 2.0 & \beta = 0.5\\ p = 0.5 & \theta = 0.6\\ k = 1.5 & \alpha = 0.24\\ h = 1.2 & m = 0.34 \end{array} \] Noise intensities:
\[ \sigma = \begin{bmatrix} 10^{-2} & 10^{-2}\\ 10^{-2} & 10^{-2} \end{bmatrix} \] White-Shot:
\[ \begin{array}{ll} \lambda = 2.0 & \text{range: } (-1.5, 1.5) \end{array} \]

White-Shot Noise Predator extinction

Parameters:
\[ \begin{array}{ll} r_0 = 2.0 & \beta = 0.5\\ p = 0.5 & \theta = 0.6\\ k = 1.5 & \alpha = 0.24\\ h = 1.2 & m = 0.34 \end{array} \] Noise intensities:
\[ \sigma = \begin{bmatrix} 10^{-2} & 10^{-2}\\ \mathbf{0.6} & 10^{-2} \end{bmatrix} \] White-Shot:
\[ \begin{array}{ll} \lambda = 2.0 & \text{range: } (-1.5, 1.5) \end{array} \]

White-Shot Noise Predator Repopulation Attempt

Parameters:
\[ \begin{array}{ll} r_0 = 2.0 & \beta = 0.5\\ p = 0.5 & \theta = 0.6\\ k = 1.5 & \alpha = 0.24\\ h = 1.2 & m = 0.34 \end{array} \] Noise intensities:
\[ \sigma = \begin{bmatrix} 10^{-2} & 10^{-2}\\ \mathbf{0.6} & 10^{-2} \end{bmatrix} \] White-Shot:
\[ \begin{array}{ll} \lambda = 1.25 & \text{range: } (\mathbf{0}, 0.5) \end{array} \]

Thank you for your attention!