\documentclass[reqno]{amsart}
\usepackage{a4wide}
\usepackage{hyperref}
\usepackage{graphicx}

\AtBeginDocument{{\noindent\small
\emph{Electronic Journal of Differential Equations},
Vol. 2026 (2026), No. 44, pp. 1--28.\newline
ISSN: 1072-6691. URL: https://ejde.math.txstate.edu, DOI: 10.58997/ejde.2026.44}
\thanks{\copyright 2026. This work is licensed under a CC BY 4.0 license.}
\vspace{8mm}}

\begin{document}
\title[\hfilneg EJDE-2026/44\hfil Socio-environmental ODE models]
{Socio-environmental models with low-dimensional nonlinear ODEs}

\author[M. Badiale, I. Cravero \hfil EJDE-2026/44\hfilneg]
{Marino Badiale, Isabella Cravero}

\address{Marino Badiale \newline
Department of Mathematics,
University of Turin, Italy}
\email{marino.badiale@unito.it}

\address{Isabella Cravero \newline
Department of Mathematics,
University of Turin, Italy}
\email{isabella.cravero@unito.it}

\thanks{Submitted March 24, 2026. Published June 23, 2026.}
\subjclass[2020]{34A34, 34D20, 92D25}
\keywords{Ordinary differential equations; stability analysis;
 socio-environmental models}

\begin{abstract}
We study a low-dimensional system of ordinary differential equations
modeling the interactions between nature and society.
The variables are renewable and non-renewable resources
(nature), and population, wealth and pollution (society).
We find equilibria for the system and study their stability.
Also we obtain several results on the stability of trajectories,
possible collapse trajectories, and societal ``safe-harbors.''
\end{abstract}

\maketitle
\numberwithin{equation}{section}
\newtheorem{theorem}{Theorem}[section]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{claim}[theorem]{Claim}
\allowdisplaybreaks


\newcommand{\CC}{\mathbb{C}}
\newcommand{\RR}{\mathbb{R}}
\newcommand{\BigO}{\mathcal{O}}

\section{Introduction}

In this article we introduce low-dimensional  systems of ordinary differential equations
intended to model the interactions between nature and humanity, 
the former through renewable and non-renewable resources, and the latter through the 
levels of population, wealth and pollution.
Research on mathematical models of this kind started at least with the 
``World3 model'' used by the MIT group that wrote the famous book 
\emph{Limits to Growth} \cite{meadows1972limits}. 
This work generated much discussion and many analyses
that are still ongoing: see \cite{Vermeulen1976,Thissen1978} for early analysis,
\cite{meadows2004limits} for a thirty-year update, and
\cite{herrington2020update,Nebel2024,Nebel2025Corrigendum} for recent discussions.

More recently, the concern about ecological crises, and in particular climate change,
led many researchers to introduce many other different types of models for the
interactions between nature and human activities;
see for example the review articles \cite{Calvin2018,Tan2023}.
Such models may be high-dimensional, as in the case of World3, and are therefore
studied primarily by means of numerical simulations.
Another research direction concerns the development of low-dimensional models.
Their main advantage lies in their simplicity, which allows the application of the
qualitative theory of ODE systems. On the other hand, precisely because
of this simplicity, they cannot describe real-world systems with high accuracy.
Consequently, they are not suitable for quantitative predictions, but they can provide
insight into the qualitative behavior of the system under study.
A low-dimensional model of this kind is the HANDY (Human and Nature Dynamics) model,
introduced in 2014 \cite{MotRivKal14} and subsequently developed by several authors
(see, for example, \cite{AkYor20,GraWilSat20,Ton23}).
Another example of such a low-dimensional system is given in \cite{Nitzbon}.

 The starting point of our research was the HANDY model.
In our paper \cite{BaCra24a}, following \cite{ShiKad23}, we distinguished between
renewable and non-renewable resources, with the new feature of the introduction of a
term describing the replenishment of non-renewable resources due to human ingenuity.
Subsequently, we proposed several modifications of the original model, thereby moving
beyond the HANDY framework; see \cite{BaCra24b,BaCra25a,BaCra25b}.

In this article, we build upon the model developed in \cite{BaCra25b}
by introducing pollution and analyzing its impact on nature and society.
Pollution is a relevant variable in several models, starting of course with World3,
but our approach to modeling it was influenced by a paper by Perissi
(see \cite{Perissi2019}, in particular the equations on page 9) and an article by
Bardi (see \cite{Bardi2024PopulationCollapse}, in particular the ``Seneca equations"
in the Appendix). Differently from previous works, a relevant feature of the present
paper is that we consider wealth the main driver of societal dynamics, while population
becomes a variable which does not influence the others but is influenced by them.

 As in our previous papers, we denote by \( x \) the population, by \( y \)
the renewable resources, by \( z \) the non-renewable resources, and by
\( w \) the wealth derived from them.
In the present paper, we introduce a new variable \( p \) to account for pollution.

Let us now introduce the equations of the system and provide some
explanations.

 (i) The equation for $x$ is
 $$
 x' = -c x + d \frac{\mu x w}{w^2 + \mu^2 x^2} \; x - \gamma_1 p x ,
 $$
where $c,d, \mu , \gamma_1 $ are positive parameters.
 Notice that
$$
 \frac{\mu x w}{w^2 + \mu^2 x^2} = \frac{\frac{w}{\mu x}}{\frac{w^2}{\mu^2 x^2} +1 },
$$
and the function $f(t) = t/(t^2 +1 )$ tends to zero when $t\to 0$ or $t\to +\infty$.
Hence, the population decreases when the wealth per person ($w/x$) is very low or
very high, while increasing in an intermediate range. This seems to fit the evolution of
human societies. The term $-\gamma_1 p x$ means of course that pollution acts
negatively on the growth of human population. Notice that we have to assume $d>2 c$,
or even stronger hypotheses, otherwise it should be $x'(t) \leq 0$ for all
$t \geq 0$, which is not reasonable.


(ii) The equation for $y$ is
 $$
 y' = \gamma y (\lambda - y) - \delta_1 w y - \gamma_2 p y,
 $$
 where $ \gamma, \lambda, \delta_1 , \gamma_2 $ are positive parameters.
In this equation a standard logistic term is contrasted by the term $- \delta_1 w y$
that describes human exploitation based on wealth $w$, and by the term $- \gamma_2 p y$
that describes the damages caused by pollution to natural resources.

(iii) The equation for $z$ is
 $$
 z' = - \delta_2 w z + k_1 \frac{\delta_3 w }{\delta_3 w + 1} \; z - \gamma_3 p z.
 $$
where $\delta_2 , \delta_3 , k_1 , \gamma_3 $ are positive parameters.
Here the terms $- \delta_2 w z $ and $- \gamma_3 p z$ have the same meaning as in the
equation for $y$, while the term $k_1 \frac{\delta_3 w }{\delta_3 w + 1} $ means
that there is some level of replenishment of non-renewable resources due to human
ingenuity and the use of wealth. The rate of replenishment is bounded, so this term
expresses a kind of moderate optimism about non-renewable resources.
This moderate optimism takes into account some criticism raised against World3 and
similar models since their first appearance.

(iv) The equation for $p$ is
 $$
 p' = \gamma_4 w p - k_2 \frac {\delta_4 w} {\delta_4 w +1} \; p - \gamma_5 p \, y.
 $$
where $\gamma_4 , \gamma_5 , k_2 , \delta_4 $ are positive parameters.
The term $\gamma_4 w p$ means that pollution increases with wealth, while the other
two terms model counteractions to pollution: $- \gamma_5 p \, y$ stands for the
removal of pollution through the action of nature, while
$- k_2 \frac {\delta_4 w} {\delta_4 w +1} $ model human action based on wealth,
also in this case with a bounded rate of growth.

(v) The equation for $w$ is
 $$
 w' = \delta_1 wy + \delta_2 w z - \delta_3 w -\delta_4 w - s w.
 $$
The wealth increases thanks to the resources $y,z$ and decreases due to 
investments in the search for non-renewable resources and in counteracting pollution
and for direct consumption given by $-sw$ ($s$ is a positive parameter).

Putting together these equations, we obtain a 5x5 nonlinear ODE system.
The large number of parameters makes its study difficult, hence we will pass to some
simplified versions. First, we assume $\delta_3 = \delta_4$, then we apply some
standard rescaling, and after a suitable renaming of variables and parameters we
obtain the following model:
\begin{equation} \label{eq:second}
 \begin{gathered}
 x'  = -c x + d \frac{ x w}{w^2 +x^2} \; x - \gamma_0 p x, \\
 y'  = \gamma y (\lambda - y) - \delta_1 w y - \gamma_1 p y, \\
 w'  = \delta_1 wy + \delta_2 w z - w, \\
 z'  = - \delta_2 w z + k_1 \frac{ w }{ w + 1} \; z - \gamma_2 p z, \\
 p'  = \gamma_3 w p - k_2 \frac { w} { w +1} \; p - \gamma_4 p \, y
 \end{gathered}
 \end{equation}
Since the $x$-variable does not appear in the equations for the other variables,
we will study the fourth-order system consisting of the equations for $y,w,z,p$.
The equation for $x$ will be studied separately. The resulting 4x4 system is the
following
\begin{equation} \label{eq:secondnonx}
 \begin{gathered}
 y'  = \gamma y (\lambda - y) - \delta_1 w y - \gamma_1 p y, \\
 w'  = \delta_1 wy + \delta_2 w z - w, \\
 z'  = - \delta_2 w z + k_1 \frac{ w }{ w + 1} \; z - \gamma_2 p z, \\
 p'  = \gamma_3 w p - k_2 \frac { w} { w +1} \; p - \gamma_4 p \, y
 \end{gathered}
 \end{equation}

In the study of \eqref{eq:secondnonx} the numerous parameters still create
difficulties, hence some of our results will be provided for a simplified form of the
system, as follows. According to \cite{Ton23}, we assume
$\gamma= \frac 12, \, \lambda=100$, and to reduce the number of parameters we set
$\delta_1=\delta_2=\delta, \, \gamma_1 = \alpha$,
$\gamma_2 = \frac 12 \alpha, \, \gamma_3=\gamma_4= 2 \alpha, \, k_1=k_2=k$.
We then obtain the following model, depending on the three parameters
$\alpha, \delta , k$:
 \begin{equation} \label{eq:first_simply}
 \begin{gathered}
 y'  = \frac 12 y (100 - y) - \delta w y - \alpha p y, \\
 w'  = \delta wy + \delta w z - w, \\
 z'  = - \delta w z + k \frac{ w }{ w + 1} \; z - \frac 12 \alpha p z, \\
 p'  =2 \alpha w p - k \frac { w} { w +1} \; p - 2 \alpha p \, y
 \end{gathered}
 \end{equation}

 The main results we will obtain for systems \eqref{eq:secondnonx} or
\eqref{eq:first_simply} deal with the computation of equilibrium points and the study
of their stability properties. As we will see, in some cases our results can be viewed
as collapse results: it may happen that, when a trajectory of the system approaches an
equilibrium with $w=0$, wealth and population collapse to zero.

Another relevant result is the fact that, in our model, variables like $w(t)$ cannot
diverge to $+\infty$ when $t\to + \infty$. This means that a steady economic growth
is not possible, in our model. These results agree with what we obtained in our
previous papers, confirming their robustness (in a non-technical sense).
However,  we prove the existence of stable equilibrium points with positive values for
population, wealth and natural resources. This means that, if a steady growth is not
possible, there are chances to drive the society towards ``safe harbors", avoiding
collapse. Again, this result agrees with our previous results, but now the effect of
pollution seems to make it much more difficult to obtain these ``safe harbors", which
are possible only for very particular ranges of parameters.

This paper is organized as follows: after the introduction (section 1),
section 2 gives some general properties of the trajectories, concerning in
particular their asymptotics. In section 3 we compute the equilibria of system
\eqref{eq:secondnonx} or, in some cases, system \eqref{eq:first_simply}. In sections 4 and 5 we study the stability of the equilibria, obtaining almost complete results in some cases, but not in all: sometimes we have to settle for partial results, in particular for asymptotic ones. In section 6 we come back to the study of the population variable $x$, and in particular of the stability of its stationary points. In section 7 we show some numerical simulations which illustrate the theoretical results. In section 8 we summarize the main results of the paper.

 \section{General results}

In this section we study system \eqref{eq:secondnonx}.
We will work in the positive cone of $\mathbb{R}^4 $, that is
$$
 \mathcal{C} = \left\{ ( y, w, z, p) \in \mathbb{R}^4 : y > 0, \, w > 0,
 \, z > 0, \, p >0 \right\}
$$
or in its closure
$$
 \overline{\mathcal{C}}= \left\{ (y, w, z, p) \in \mathbb{R}^4
 : y \geq 0, \, w \geq 0, \, z \geq 0, \, p \geq 0 \right\}.
$$


Let us write \(X = (y, w, z, p) \in \mathbb{R}^4\) and let \(F(X)\) be the vector
field of the ODE system \eqref{eq:secondnonx}. We define
$$
 \mathcal{S} = \left\{ ( y, w, z, p) \in \mathbb{R}^4 : w > -1, \right\} 
$$
and we obtain $\overline{\mathcal{C}} \subset \mathcal{S} $ and
\(F \in C^{\infty}(\mathcal{S})\). Hence, by standard ODE theory, and applying the
same arguments as in \cite{BaCra25b}, we have the following result.

\begin{proposition}
For any \(t_0 \in \RR, \, X_0=(y_0, \, w_0, \, z_0,\, p_0) \in
\overline{\mathcal{C}}\), the problem
\begin{equation} \label{eq:sistemaprincipale}
 \begin{gathered}
  X'(t) = F(X) \\
  X(t_0) = X_0
 \end{gathered}
\end{equation}
 has a unique solution defined in a maximal interval \(J=(a,+\infty )\) with
\(a < t_0 \).
 \end{proposition}

 We will study the solutions of \eqref{eq:sistemaprincipale} with \(t_0=0\).
Arguing as in \cite{BaCra25b} it is easy to prove that $\mathcal{C} $ and
$\overline{\mathcal{C}}$ are invariant sets, and we state this result as follows.

\begin{proposition} \label{prop:due}
Let \(X(t) = (y(t),\, w(t),\, z(t),\, p(t)) \in C^1((a,+\infty))\) be the solution
of \eqref{eq:sistemaprincipale} with \( -\infty \le a < 0 \).
Let \(\zeta(t)\) be any of the four components. Then the following holds:
 \begin{itemize}
 \item[(i) ] If \(\zeta_0 = \zeta(0) > 0\), then \(\zeta(t) > 0\) for every
 \(t \in (a,+\infty )\);

 \item[(ii)] If \(\zeta_0 = 0\), then \(\zeta(t) = 0\) for every
 \(t \in (a,+\infty )\).
 \end{itemize}
 \end{proposition}


We then prove a result on the $y$-component of solutions.

 \begin{proposition}  \label{prop:esty}
Let \(X(t)\) be as in Proposition \ref{prop:due} with \(y_0 \geq 0\). Then
\[
 0 < y(t) \le \max\{y_0, \lambda\}, \quad \forall t \in [0,+\infty ).
\]
\end{proposition}

  The proof of the above proposition is the same as in \cite[Lemma 7]{BaCra24a}.
 We now study the asymptotic behavior of the solutions. Let $X(t)$ be as in
 Proposition \ref{prop:due}.

 \begin{proposition} \label{pro:w}
It cannot happen that \(w(t) \to +\infty\) as \(t \to +\infty\).
 \end{proposition}

\begin{proof}
Assume by contradiction that \(w(t) \to +\infty\) as \(t \to +\infty\).
Then the equation for $y$ implies
 \[
 y' = (\gamma(\lambda - y) - \delta_1 w - \gamma_2 p)\, y
 \le (\gamma \lambda - \delta_1 w)\, y.
 \]
Since \(\gamma \lambda - \delta_1 w \to -\infty\), we may assume
\(y' \le -y\) on \((\tau_1, +\infty)\), which implies \(y(t) \to 0\) as
\(t \to +\infty\).
 Similarly,
 \[
 z' \le (-\delta_2 w + k_1) z \le -z
 \]
on \((\tau_2, +\infty)\), hence \(z(t) \to 0\) as \(t \to +\infty\).
We have
 \[
 w' = (\delta_1 y + \delta_2 z - \sigma)\, w,
 \]
where \(\sigma = \delta_3 + \delta_4 + s\).
Since \(\delta_1 y + \delta_2 z \to 0\), we obtain
\(w' \le -\frac{\sigma}{2}\, w\) on \((\tau_3, +\infty)\), hence \(w(t) \to 0\)
as \(t \to +\infty\).
This is a contradiction with the hypothesis \(w(t) \to +\infty\), and this proves
the claim.
\end{proof}

\begin{proposition} \label{pro:z}
It cannot happen that \(z(t) \to +\infty\) as \(t \to +\infty\).
\end{proposition}

\begin{proof}
Assume by contradiction that \(z(t) \to +\infty\). Then from
\[
 w' \ge (\delta_2 z - \sigma)\, w
\]
we obtain \(w' \ge w\) on \((\tau, +\infty)\), hence \(w(t) \to +\infty\).
This is impossible by Proposition \ref{pro:w}.
 \end{proof}

\begin{proposition}
 It cannot happen that \(p(t) \to +\infty\) as \(t \to +\infty\).
\end{proposition}

\begin{proof}
Assume by contradiction that \(p(t) \to +\infty\).
Then, by Proposition~\ref{pro:w}, we have
\[
 y' \le (\gamma \lambda - \gamma_2 p)\, y \le -y
 \quad\text{and}\quad
 z' \le (k_1 - \gamma_3 p)\, z \le -z
 \quad\text{on } (\tau_1, +\infty),
 \]
hence \(y(t), z(t) \to 0\) as \(t \to +\infty\).
 Then, as above, setting $\eta =\sigma/2, $ we obtain that there exists
 $\tau_2 >\tau_1 $ such that
\[
 w'(t) \le -\eta\, w(t)
 \quad\text{on } (\tau_2, +\infty),
\]
hence \(w(t) \le C e^{-\eta t}\) on \((\tau_2, +\infty)\).
Thus, for $t> \tau_2 $, we obtain
\[
 p'(t) \le \gamma_4 C e^{-\eta t} p(t),
\]
and therefore
\[
 \frac{p'(t)}{p(t)} \le c_1 e^{-\eta t},
 \quad\text{i.e.}\quad
 \frac{d}{dt} \log p(t) \le c_1 e^{-\eta t}.
\]
Integrating on \((\tau_2, t)\) gives
\[
 \log \frac{p(t)}{p(\tau_2)}
 \le \frac{c_1}{\eta}\,(e^{-\eta \tau_2} - e^{-\eta t})
 \le c_2 e^{-\eta \tau_2},
\]
 so
\[
 p(t) \le p(\tau_2)\, \exp\!\bigl(c_2 e^{-\eta \tau_2}\bigr)
 \quad\text{for every } t > \tau_2.
\]
Hence \(p(t)\) is bounded, giving a contradiction.
\end{proof}

\begin{remark}\label{re:osc} \rm
These results do not imply that $w,z,p$ are bounded, but that,
if unbounded, they will show strong oscillations, which in general are not good
news for a society.
\end{remark}

  We now prove a link between the asymptotic of $w$ and that of $x$, assuming that
$(x, y, w, z, p)$ is a solution of system \eqref{eq:second}.

\begin{proposition} \label{prop:wx}
If \(w(t) \to 0\) as \(t \to +\infty\), then \(x(t) \to 0\) as \(t \to +\infty\).
\end{proposition}

\begin{proof}
By standard ODE theory, it is easy to prove that $x(t)$ is defined for any
$t \geq 0$, and that either $x(t)$ is the constant zero function or $x(t)>0$
for all $t\geq 0$. The thesis is trivial if $x(t)$ is the zero constant, so we assume
$x(t)>0$ for all $t \geq 0$. Assume \(w(t) \to 0\) as \(t \to +\infty\) and fix
\(\varepsilon > 0\).
Then there exists \(T_\varepsilon > 0\) such that
 \[
 w(t) < \frac12 \frac{c}{d}\,\varepsilon
 \quad \text{for every } t > T_\varepsilon,
 \]
where $c,d$ are the parameters in the equation for $x$. Now we need the following
claim.


\begin{claim}
It cannot happen that \(x(t) \ge \varepsilon\) for every \(t > T_\varepsilon\).
\end{claim}

\begin{proof}
 Indeed, if this were the case, then
 \[
 -c + d\,\frac{w}{x}
 \le -c + d\,\frac12 \frac{c}{d}\,\varepsilon\,\frac{1}{\varepsilon}
 = -\frac{c}{2}.
 \]
Hence, from the equation for $x$, we obtain
 \[
 x'
 = -cx + d\,\frac{\frac{w}{x}}{\left(\frac{w}{x}\right)^2 + 1}\,x - \gamma px
 \le -cx + d\,\frac{w}{x}\,x
 = \left(-c + d\,\frac{w}{x}\right)x
 \le -\frac{c}{2} x
 \]
for every \(t > T_\varepsilon\). Integrating from \( T_\varepsilon\) to \(t \),
we obtain
\[
 x(t) \le x(T_\varepsilon)\,\exp\!\left(-\frac{c}{2}(t - T_\varepsilon)\right),
\]
and therefore \(x(t) \to 0\) as \(t \to +\infty\), contradicting the assumption that
\(x(t) > \varepsilon\) for all \(t > T_\varepsilon\).
 This contradiction proves the claim.
 \end{proof}

We have then proved that there must exist \(s_\varepsilon > T_\varepsilon\) such that
\(x(s_\varepsilon) < \varepsilon\). The proof will be completed by a second claim,
as follows.

\begin{claim}
 One has \(x(t) < \varepsilon\) for every \(t > s_\varepsilon\).
\end{claim}

\begin{proof}
 Let
\[
 \tau_{\varepsilon}
 = \sup \left\{\, t > s_{\varepsilon} \;\middle|\; x(s)
 < \varepsilon,\ \forall\, s \in [s_{\varepsilon}, t) \right\}.
\]
We have \(x(t) < \varepsilon\) for every \(t \in [s_{\varepsilon},
\tau_{\varepsilon})\).
If \(\tau_{\varepsilon} = +\infty\), the claim is proved.
We argue by contradiction and assume \(\tau_{\varepsilon} < +\infty\).
Hence, by the definition of \(\tau_{\varepsilon}\) together with the usual continuity
arguments, we obtain
\[
 x(\tau_{\varepsilon}) = \varepsilon .
\]
Since \(x(t) < \varepsilon\) for all \(t \in [s_{\varepsilon}, \tau_{\varepsilon})\),
we must have
\[
 x'(\tau_{\varepsilon}) \ge 0.
 \]
On the other hand, \(\tau_{\varepsilon} > s_{\varepsilon} > T_{\varepsilon}\),
and we know that
 \(w(t) < \frac{1}{2}\frac{c}{d}\varepsilon\) for all \(t > T_{\varepsilon}\).
 Then, as before,
\[
 x'(\tau_{\varepsilon})
 \le \Big(-c + d\,\frac{w(\tau_{\varepsilon})}{x(\tau_{\varepsilon})}\Big)
  x(\tau_{\varepsilon})
 \le \Big(-c + d\,\frac12\frac{c}{d}\varepsilon\,\frac{1}{\varepsilon}\Big)\varepsilon
 = -\frac{c}{2}\varepsilon < 0.
\]
This contradiction shows that \(\tau_{\varepsilon} = +\infty\), completing the proof.
\end{proof}

In conclusion, we have proved that for every \(\varepsilon > 0\) there exists
\(s_{\varepsilon} > 0\) such that for every \(t > s_{\varepsilon}\) one has
\(x(t) < \varepsilon\).
 Recalling that $x(t)>0$, we obtain, by the very definition of limit, that
 \[
 \lim_{t \to +\infty} x(t) = 0.
 \]
 \end{proof}


\section{Equilibrium points}

In this section we study the equilibria of \eqref{eq:secondnonx}.
They satisfy the  system
 \begin{equation}
 \label{eq:punti_critici}
 \begin{gathered}
  y ( \gamma \lambda - \gamma y - \delta_1 w - \gamma_1 p )=0, \\
  w ( \delta_1 y + \delta_2 z - 1 ) =0, \\
  z \Big( -\delta_2 w + k_1 \frac{w}{w+1} - \gamma_2 p \Big) = 0, \\
  p \Big( \gamma_3 w - k_2 \frac{w}{w+1} - \gamma_4 y \Big) = 0.
 \end{gathered}
 \end{equation}
For each equation we have two possibilities, hence \(2^4 = 16\) possible cases.
The computations are straightforward, so we skip them and we just give a list of the
equilibria. We obtain three half-lines of equilibria and some isolated ones, as follows.
\begin{itemize}
\item[(1)] Half-line
 $\left \{ (0,0,0,\tilde{p}) \in \RR^4 : \tilde{p} \ge 0 \right \}$.

 \item[(2)] Half-line
 $\left \{ (0,0,\tilde{z}, 0 ) \in \RR^4 : \tilde{z} \ge 0 \right \}$.

 \item[(3)] Half-line
 $\left \{ (\lambda,0,\tilde{z}, 0 ) \in \RR^4 : \tilde{z} \ge 0 \right \}$.

 \item[(4)] Isolated equilibrium
 $ \left(0,\,\frac{k_1}{\delta_2} - 1,\,\frac{1}{\delta_2},\,0\right)$,
  with the condition \(k_1 \ge \delta_2\).

 \item[(5)] Isolated equilibrium
 $ \Big( 0,\, \frac{k_2}{\gamma_3} - 1,\,
 \frac{1}{\delta_2},\,
 \frac{\delta_2}{\gamma_2}\!\big(1 - \frac{k_2}{\gamma_3}\big)
 + \frac{k_1}{\gamma_2}\!\big(1 - \frac{\gamma_3}{k_2}\big)
 \Big)$.

 \item[(6)] Isolated equilibrium
$ \Big( \frac{1}{\delta_1},\; \frac{\gamma}{\delta_1^2}(\lambda \delta_1 - 1),
\; 0,\; 0 \Big)$.

 \item[(7)] Isolated equilibrium
 $ \Big( \frac{1}{\delta_1},\; w_0,\; 0,\; p_0 \Big)$,
 where
 \(\eta := \gamma_4/\delta_1\) and
 \[
 w_0 = \frac{\eta + k_2 - \gamma_3 + \sqrt{(\eta + k_2 - \gamma_3)^2
  + 4 \eta \gamma_3}}{2\gamma_3}, \quad
   p_0 = \frac{1}{\gamma_1}\Big(\gamma \lambda - \frac{\gamma}{\delta_1}
   - \delta_1 w_0\Big),
 \]


 \item[(8)] Isolated equilibrium
 $ \Big(
 \lambda - \delta_1 \frac{k_1 - \delta_2}{\gamma \delta_2},\;
 \frac{k_1}{\delta_2} - 1,\;
 z_0,\;
 0
 \Big)$,
 where
\[
 z_0 = \frac{1}{\delta_2}\big[ 1 - \delta_1 \lambda + \delta_1^2 
 \frac{k_1 - \delta_2}{\gamma \delta_2} \big]
 = \frac{1}{\delta_2} - \frac{\delta_1}{\delta_2}\lambda
 + \frac{\delta_1^2}{\delta_2^2}\frac{k_1}{\gamma}
 - \frac{\delta_1^2}{\gamma \delta_2}.
 \]


 \item[(9)] The last case is that in which we have to solve the system
 \begin{gather*}
  \gamma \lambda - \gamma y - \delta_1 w - \gamma_1 p = 0, \\
  \delta_1 y + \delta_2 z - 1 = 0, \\
  -\delta_2 w + k_1 \frac{w}{w+1} - \gamma_2 p = 0, \\
  \gamma_3 w - k_2 \frac{w}{w+1} - \gamma_4 y = 0,
 \end{gather*}
which implies
 \begin{gather*}
  \gamma \lambda - \gamma y - \delta_1 w - \gamma_1 p = 0, \\
  z = \frac{1}{\delta_2}(1 - \delta_1 y) = \frac 1 {\delta_2}
   - \frac {\delta_1} {\delta_2 \gamma_4}\Big( \gamma_3 w - k_2 \frac{w}{w+1} \Big), \\
  p = \frac{1}{\gamma_2}\Big( -\delta_2 w + k_1 \frac{w}{w+1} \Big), \\
 y = \frac{1}{\gamma_4}\Big( \gamma_3 w - k_2 \frac{w}{w+1} \Big).
 \end{gather*}
Substituting \(y\) and \(p\) into the first equation, we obtain a
second-order equation for \(w\), whose solutions determine \(y\), \(z\), and \(p\).
The resulting expressions for the equilibria depend on the parameters in a
too complicated way to allow a complete stability analysis, so we will restrict
ourselves to studying some particular cases.
 \end{itemize}

 In the next two sections we  study the stability of these equilibria.
In some cases we will be able to obtain a fairly complete description of stability,
in other cases the large number of parameters will prevent a complete study, and we
will deal only with the simplified model \eqref{eq:first_simply}, obtaining some
asymptotic results.

 To study the stability, we write \(F = (f_1, f_2, f_3, f_4)\) for the vector
field of model \eqref{eq:secondnonx}, that is
 \begin{gather*}
 f_1 = \gamma y(\lambda - y) - \delta_1 w y - \gamma_1 p y,\quad
 f_2 = \delta_1 w y + \delta_2 w z - w,\\
 f_3 = -\delta_2 w z + k_1 \frac{w}{w+1}\, z - \gamma_2 p z,\quad
 f_4 = \gamma_3 w p - k_2 \frac{w}{w+1}\, p - \gamma_4 p\, y.
 \end{gather*}
The Jacobian matrix is then as follows
\begin{equation*} %\label{eq:jacobiana}
 J_F=\begin{pmatrix}
 \gamma \lambda- 2 \gamma y - \delta_1 w - \gamma_1 p & -\delta_1 y & 0
 & - \gamma_1 y \\ \delta_1 w & \delta_1 y + \delta_2 z -1 & \delta_2 w & 0 \\
 0 & -\delta_2 z + k_1 \frac z {(w+1)^2} & - \delta_2 w + k_1 \frac{ w }{ w + 1}
 - \gamma_2 p & -\gamma_2 \, z \\
 - \gamma_4 p & \gamma_3 p- k_2 \frac p {(w+1)^2} & 0 & \gamma_3 w - k_2
 \frac { w} { w +1} - \gamma_4 y
 \end{pmatrix}.
 \end{equation*}

\section{Half-lines of equilibria}

In this section, as in the previous one, we study system \eqref{eq:secondnonx},
and we deal with cases (1)--(3).

\subsection{Case 1}
 We study here the stability of the equilibria \(P= (0, 0, 0, \tilde{p} )\) with
 \(\tilde{p} \ge 0\).
The Jacobian matrix becomes
 \begin{equation} \label{eq:jacobianaA}
 J_F(P)=\begin{pmatrix}
 \gamma \lambda - \gamma_1 \tilde{p} & 0 & 0 & 0 \\[4pt]
 0 & -1 & 0 & 0 \\[4pt]
 0 & 0 & -\gamma_2 \tilde{p} & 0 \\[4pt]
 -\gamma_4 \tilde{p} & (\gamma_3 - k_2)\tilde{p} & 0 & 0
 \end{pmatrix}.
 \end{equation}
The eigenvalues are
 \[
 \rho_1 = \gamma \lambda - \gamma_1 \tilde{p} \quad \rho_2 = -1,\quad
 \rho_3 = -\gamma_2 \tilde{p} \quad \rho_4 = 0,
 \]
so we have a first, very obvious, instability result.

\begin{proposition}
If \(\gamma_1 \tilde{p} < \gamma \lambda\), then the equilibrium point
\(P = (0, 0, 0, \tilde{p})\) is unstable.
 \end{proposition}

Instead, if \(\gamma_1 \tilde{p} > \gamma \lambda\), then the equilibrium point
\(P \) is Lyapunov-stable, but this is definitely not a trivial result, and we are
going to prove it through several lemmas.

From now on we fix $\tilde{p}$ such that \(\gamma_1 \tilde{p} > \gamma \lambda\).
Notice that this implies $\tilde{p} >0$. We will study the behavior of a trajectory
$X(t)= (y(t), w(t), z(t), p(t))$ with initial value $X_0 = (y_0, w_0, z_0, p_0)$.

\begin{lemma}  \label{lem:y}
If we fix
 \(0< \varepsilon < \frac{1}{2\gamma_1}\,(\gamma_1 \tilde{p} - \gamma \lambda)\)
and $ p_0 >\tilde{p}-\varepsilon$, then we obtain a maximal interval
$[0, b_1 )$ such that
 \[
 y(t) \le y_0 e^{-a t} \quad \text{for } t \in [0,b_1)
 \]
with \(b_1 \in (0,+\infty] \) and
\( a= \tfrac12(\gamma_1 \tilde{p} - \gamma \lambda)>0\).
\end{lemma}

\begin{proof}
If $y_0 =0$ we know that $y(t)=0$ for all $t\geq 0$ and the thesis is true with
$b_1 = + \infty$. So let us assume $y_0>0$. From
 \[
 0 < \varepsilon < \frac{1}{2\gamma_1}\,(\gamma_1 \tilde{p} - \gamma \lambda),
 \]
 we obtain
 \[
 \gamma_1 (\tilde{p} - \varepsilon) - \gamma \lambda
 > \tfrac12(\gamma_1 \tilde{p} - \gamma \lambda)
 = a.
 \]
 From the equation for \(y\) we have
\[
 y'(t) \le y(t)\,(\gamma \lambda - \gamma_1 p(t)),
\]
 hence
\[
 y'(0) \le y_0\bigl(\gamma \lambda - \gamma_1 (\tilde{p} - \varepsilon)\bigr)
 < -a\,y_0.
\]
By continuity we have \(y'(t) < -a\,y(t)\) on some left neighborhood of $0$. We then define
\[
 b_1
 = \sup\Bigl\{\tau > 0 \,\big|\,
 y'(t) < -a\,y(t),\ \forall t \in [0,\tau)\Bigr\},
\]
and we obtain \(b_1 \in (0,+\infty]\) and \(y'(t) < -a\,y(t)\) on \([0,b_1)\), and
consequently
\[
 y(t) \le y_0 e^{-a t} \quad \text{for } t \in [0,b_1).
\]
 \end{proof}

\begin{lemma} \label{lem:w}
If we fix
\(0<\varepsilon < \frac{1}{2 \left( \delta_1 + \delta_2 \right)}\)
and $0< y_0 , z_0 < \varepsilon$, then we obtain a maximal interval $[0, b_2 )$
such that
\[
 w(t) \le w_0 e^{-t/2} \quad \text{for } t \in [0,b_2), \; b_2 \in (0,+\infty].
\]
 \end{lemma}

\begin{proof}
 If $w_0 =0$ we know that $w(t)=0$ for all $t\geq 0$ and the thesis is true with $b_2 = + \infty$. Let us assume $w_0>0$. We have
 \[
 \varepsilon(\delta_1 + \delta_2) < \tfrac12,
 \]
 and if \(z_0, y_0 < \varepsilon\), then, using the equation for $w$,
 \[
 w'(0)
 = w_0\bigl(\delta_1 y_0 + \delta_2 z_0 - 1\bigr)
 < w_0\left(-\tfrac12\right).
 \]
 Arguing as in Lemma \ref{lem:y}, with a definition for $b_2$ similar to that of
 $b_1$, we obtain a maximal interval \([0,b_2)\), with \(b_2 \in (0,+\infty]\),
 on which
 \[
 w'(t) < -\tfrac12\,w(t),
 \]
 and therefore
 \[
 w(t) \le w_0 e^{-t/2}
 \quad \text{for } t \in [0,b_2).
 \]
 \end{proof}

To prove the Lyapunov-stability of $P$, we will later introduce neighborhoods of $P$
depending on a parameter $\eta >0$, and we can assume $\eta \in (0, \tilde{p})$ as in
the following lemmas.

\begin{lemma}  \label{lem:z}
Fix $\eta \in (0, \tilde{p})$ and choose $\varepsilon >0$ such that the following
inequalities are satisfied:
 $$
 (-\delta_2 + k_1 + \gamma_2)\,\varepsilon < \gamma_2 \eta , \quad \varepsilon < \frac{1}{2} (\tilde{p} - \eta ), \quad \varepsilon < \eta .
 $$
If $w_0 < \varepsilon$ and $ p_0 > \tilde{p}- \varepsilon$, then there is a maximal
\(b_3 \in (0,+\infty]\) such that
\[
 z(t) \le z_0 e^{-\gamma_2 (\tilde{p} - \eta) t}, \quad \forall t \in [0, b_3 ).
 \]
 \end{lemma}

 \begin{proof}
If $z_0 =0$ we know that $z(t)=0$ for all $t\geq 0$ and the thesis is true with
$b_3= +\infty$. Let us assume $z_0>0$. We first notice that it may happen that
\(-\delta_2 + k_1 + \gamma_2 \le 0\), in which case the hypothesis
\[
 (-\delta_2 + k_1 + \gamma_2)\,\varepsilon < \gamma_2 \eta,
\]
is satisfied by any $\varepsilon >0$. If instead \(-\delta_2 + k_1 + \gamma_2 > 0\),
then of course $\varepsilon$ must satisfy
\[
 \varepsilon < \frac{\gamma_2 \eta}{-\delta_2 + k_1 + \gamma_2}.
\]
From the equation for \(z\) we obtain
\begin{equation} \label{eq:z}
 z' \le \bigl[(-\delta_2 + k_1) w - \gamma_2 p\bigr]\,z.
\end{equation}
If \(-\delta_2 + k_1 \le 0\), then $z' \le -\gamma_2 p\,z $,
hence
 $$
 z'(0) < -\gamma_2 p_0 \,z_0 < -\gamma_2 (\tilde{p} -\varepsilon)\,z_0
 < -\gamma_2 (\tilde{p} - \eta)\,z_0.
 $$

 If instead \(-\delta_2 + k_1 > 0\), then \(-\delta_2 + k_1 + \gamma_2 > 0\), and
 from the conditions  on \(\varepsilon\) we have
\[
 (-\delta_2 + k_1)\varepsilon - \gamma_2 (\tilde{p} - \varepsilon)
 < -\gamma_2 (\tilde{p} - \eta)
\]
Since \(w_0 < \varepsilon\) and \(p_0 > \tilde{p} - \varepsilon\), it follows
from \eqref{eq:z} that
\[
 z'(0) \leq \bigl[(-\delta_2 + k_1) w_0 - \gamma_2 p_0\bigr]\,z_0
 < \bigl[(-\delta_2 + k_1)\varepsilon - \gamma_2(\tilde{p} - \varepsilon)\bigr] z_0
 < -\gamma_2(\tilde{p} - \eta)\,z_0
\]
So in any case we obtain
$$
 z'(0) < -\gamma_2(\tilde{p} - \eta)\,z_0.
$$
Arguing as in the lemmas \ref{lem:y} and \ref{lem:w}, we obtain a maximal interval
 \([0,b_3)\) on which
\[
 z'(t) < -\gamma_2 (\tilde{p} - \eta) \,z(t),
\]
and therefore
$z(t) \le z_0 e^{ -\gamma_2 (\tilde{p} - \eta) t}$.
\end{proof}

To study the behavior of $p(t)$, we start by fixing $\varepsilon_1 \in (0,1)$
satisfying all the hypotheses of the previous lemmas, that is, for a fixed
$\eta \in (0, \tilde{p} )$,
$$
\varepsilon_1 < \frac{1}{2\gamma_1}\,(\gamma_1 \tilde{p}
- \gamma \lambda), \quad \varepsilon_1
< \frac{1}{2 \left( \delta_1 + \delta_2 \right)}, \quad
(-\delta_2 + k_1 + \gamma_2)\,\varepsilon_1 < \gamma_2 \eta , \quad
\varepsilon_1 < \frac{1}{2} (\tilde{p} - \eta ), \quad \varepsilon_1 < \eta .
$$

\begin{lemma} \label{lem:p}
There exists $0 <\varepsilon_2 <\varepsilon_1$ such that, if
$0<\varepsilon <\varepsilon_2 $ and we choose initial conditions
$w_0 ,y_0 <\varepsilon^2$ and $p_0 \in ( \tilde{p} - \varepsilon^2 ,
\tilde{p} + \varepsilon^2 )$ then, setting \(b_4 = \min\{b_1,b_2\}\), we obtain
 \[
 \tilde{p} - \eta < \tilde{p }- \frac{\varepsilon}{2} <p(t)< \tilde{p }+ \frac{\varepsilon}{2}< \tilde{p }+\eta, \quad \forall t \in [0,b_4).
 \]
 \end{lemma}

\begin{proof}
Let us fix any $\varepsilon < \varepsilon_1 <1$, so that $\varepsilon^2 < \varepsilon $
satisfies all the inequalities in the hypotheses of the previous lemmas. Since
\(\frac{w}{w+1} \le w\), the equation for \(p(t)\) and the previous lemmas yield
\[
 p' \ge \bigl(-k_2 w - \gamma_4 y\bigr) p
 \ge \bigl(-k_2 w_0 e^{-t/2} - \gamma_4 y_0 e^{-a t}\bigr) p,
\]
for all \(t\in [0,b_4)\), with $a$ defined in Lemma \ref{lem:y}.
Integrating we obtain
\[
 \log\frac{p(t)}{p_0}
 \ge 2k_2 w_0 (e^{-t/2} - 1)
 + \frac{\gamma_4 y_0}{a}\,(e^{-a t} - 1)
 \ge -2k_2 w_0 - \frac{\gamma_4 y_0}{a},
\]
for \(t \in [0,b_4)\). Thus
\[
 p(t) \ge p_0 \exp\!\left(-2k_2 w_0 - \frac{\gamma_4 y_0}{a}\right)
 \quad \text{for } t \in [0,b_4).
\]
From the hypotheses on $y_0, w_0$ we then obtain
 \[
 \exp\!\left(-2k_2 w_0 - \frac{\gamma_4}{a} y_0\right)
 > \exp\!\left(-\Bigl(2k_2 + \frac{\gamma_4}{a}\Bigr)\varepsilon^2\right)
 = \exp(-C \varepsilon^2),
 \]
 with \(C = 2k_2 + \gamma_4/a\), and hence, as $ p_0 > \tilde{p} - \varepsilon^2$,
 \[
 p(t)> (\tilde{p} - \varepsilon^2)\,\exp(-C \varepsilon^2)
 = (\tilde{p} - \varepsilon^2)\bigl(1 - C \varepsilon^2 + \BigO(\varepsilon^4)\bigr)
 = \tilde{p}+ \BigO(\varepsilon^2).
 \]
 when $t \in [0,b_4)$.

Now we notice that the term $\BigO(\varepsilon^2)$ is independent from $t$,
hence it is possible to choose $\varepsilon_2 \in (0, \varepsilon_1 )$ such that,
if $\varepsilon < \varepsilon_2$, it holds
 \begin{equation}
 \label{eq:p1}
 \tilde{p} + \BigO(\varepsilon^2) > \tilde{p}- \varepsilon/2 > \tilde{p} -\eta .
 \end{equation}
To obtain the other inequality we notice that, from Lemma \ref{lem:w}, we have
 \[
 p' (t) \le \gamma_3 w(t) p(t) < \gamma_3 \varepsilon^2 e^{-t/2} p(t),
 \]
 in $[0, b_4 ) $, hence
 \[
 \frac{p'(t)}{p(t)} < \gamma_3 \varepsilon^2 e^{-t/2},
 \]
 so that
\[
 \log\frac{p(t)}{p_0}
 \le
 2\,\gamma_3 \varepsilon^2.
\]
 Hence
\[
 p(t)
 < p_0 \exp(2\gamma_3 \varepsilon^2)
 < (\tilde{p} + \varepsilon^2)\bigl(1 + 2\gamma_3 \varepsilon^2
 + \BigO(\varepsilon^4)\bigr)
 = \tilde{p}+ \BigO(\varepsilon^2).
\]
Thus it suffices to choose \(\varepsilon_2\) as to satisfy, for all
$\varepsilon < \varepsilon_2$, both \eqref{eq:p1} and
\begin{equation}
 \tilde{p} + \BigO(\varepsilon^2) < \tilde{p}+ \varepsilon/2 < \tilde{p} +\eta .
\end{equation}
\end{proof}

We now prove that $b_i = + \infty$.


\begin{lemma}\label{le:beqinfty}
 Fix $ \varepsilon_2$ as in Lemma \ref{lem:p} and $\varepsilon < \varepsilon_2$.
Assume $y_0, w_0, p_0$ are chosen as in Lemma \ref{lem:p}, then
 $$
 b_1 = b_2 = b_3 =b_4 = +\infty.
 $$
 \end{lemma}

%\label{le:infty}
\begin{proof}
We first prove that $b_4 = +\infty$ and of course, by definition, this implies
$b_1 = b_2 = +\infty$. We argue by contradiction, so we assume \(b_4 < +\infty\).
We distinguish two cases:

(i) \(b_4 = b_1\). In this case \(b_1 < +\infty\) and we have, by
 Lemma \ref{lem:p},
 \(p(t) > \tilde{p} - \frac{\varepsilon}{2}\) on \([0,b_4)\), hence
 \[
 p(b_4) \ge \tilde{p} - \frac{\varepsilon}{2} > \tilde{p} - \varepsilon.
 \]
 As $b_4 = b_1$, by continuity we have a neighborhood
 \((b_1 - \tilde{\delta}, b_1 + \tilde{\delta})\) in which it holds
$$
 p(t) > \tilde{p} - \varepsilon.
$$
This implies that it is possible to repeat the arguments of Lemma \ref{lem:y}
and to obtain \(y'(t) < -a\,y(t)\) in \((b_1 - \tilde{\delta}, b_1 + \tilde{\delta})\),
hence in all $[0, b_1 + \tilde{\delta} )$. This is in contradiction with the definition
of $b_1$, so the case $b_4 = b_1 $ is ruled out.

(ii) \(b_4 = b_2 \le b_1\). In this case we recall that the equation for $z$ implies
\[
 z'(t)
 \le \bigl((-\delta_2 + k_1)w(t) - \gamma_2 p(t)\bigr) z(t).
\]
If $-\delta_2 + k_1 \leq 0$, then $z'(t) \leq 0$ for all $t$.
If $-\delta_2 + k_1 > 0$ then, recalling Lemmas \ref{lem:w} and \ref{lem:p},
we obtain
\[
 z'(t)
 \le \bigl((-\delta_2 + k_1)w(t) - \gamma_2 p(t)\bigr) z(t)
 < \bigl((-\delta_2 + k_1)\varepsilon - \gamma_2 (\tilde{p} - \varepsilon)\bigr) z(t) < 0,
\]
because we have chosen \(\varepsilon < \varepsilon_2 \) and $\varepsilon_2$ satisfies
\[
 (-\delta_2 + k_1 + \gamma_2)\varepsilon_2 < \gamma_2 \eta <\gamma_2 \tilde{p}.
\]
Hence, in any case, \(z\) is decreasing on \([0,b_2)\), and so is \(y\),
by Lemma \ref{lem:y} and $b_4 \leq b_1$.
As \(\varepsilon < \frac{1}{2(\delta_1 + \delta_2)}\), we can choose $c_0 >0$
such that
$(\delta_1 + \delta_2)\varepsilon < \frac{1}{2} -c_0 $. We then obtain
\[
 w'(t)
 = (\delta_1 y(t) + \delta_2 z(t) - 1)\,w(t)
 < \bigl((\delta_1 + \delta_2)\varepsilon - 1\bigr) w(t)
 < - \left(\tfrac{1}{2} + c_0 \right)\,w(t).
\]
Then
\[
 w'(b_4) \leq -\left(\tfrac{1}{2} + c_0 \right)\,w(b_4) < -\tfrac{1}{2}\,w(b_4),
\]
hence, repeating the argument used for case $i)$, we obtain
\(w'(t) < -\tfrac{1}{2}w(t)\) in a neighborhood
\((b_2 - \tilde \delta, b_2 + \tilde \delta)\), in contradiction the definition
of \(b_2\).
Thus, when $\varepsilon < \varepsilon_2$, the assumption \(b_4 < +\infty\) leads
in every case to a contradiction, therefore
 \(b_4 = +\infty\), and \(b_1 = b_2 = +\infty\) as well.

Let us prove $b_3 = +\infty$. From the previous results we know that
$p(t) > \tilde{p} - \eta$ and $w(t) \leq w_0 < \varepsilon^2 < \varepsilon$ for all
$t>0$. Hence we can argue from \eqref{eq:z} exactly as in Lemma \ref{lem:z}
to obtain $z'(t) < -\gamma_2 (\tilde{p} -\eta )z(t)$ for all $t\geq 0$, which means
$b_3 = + \infty$.
 \end{proof}

The lemmas we have proved up to now imply the Lyapunov-stability of the equilibrium
$P= (0,0,0, \tilde{p})$. To state this property, we consider a particular neighborhood
basis for $P$, as follows.  For any \(\eta > 0\) define
\[
 I_{\eta} = (-\eta, \eta) \times (-\eta, \eta) \times (-\eta, \eta) \times (\tilde{p} - \eta,\, \tilde{p} + \eta),
\]
The following proposition is an easy consequence of the previous lemmas.

 \begin{proposition}  \label{prop:P}
If $P=(0,0,0, \tilde{p })$ and $\gamma_1 \tilde{p } > \gamma \lambda $, then the
equilibrium $P$ is Lyapunov-stable (in $\overline{\mathcal{C}}$), that is:
for every neighborhood $V$ of $P$ there exists a neighborhood \(U \subset V\) of \(P\)
such that, if $Q \in U \cap \overline{\mathcal{C}}$ and \(X(t,Q)\) is the trajectory with
\(X(0,Q) = Q\), then
\(X(t,Q) \in V \) for every \(t \ge 0\).
 \end{proposition}

\begin{proof}
We take $\eta >0$ such that $I_{\eta} \subset V$. We can assume without loss of
generality $\eta \in (0, \tilde{p})$. We then choose $\varepsilon_2$ and
$\varepsilon$ as in Lemma \ref{le:beqinfty}, and we take
$Q \in I_{\epsilon^2} \subset I_{\eta}$. All the hypothesis of the previous lemmas
are satisfied, so we obtain $X(t,Q) \in I_{\eta}$ for all $t\geq 0$. So, the thesis
is proved with $U=I_{\varepsilon^2}$.
\end{proof}

We cannot state the asymptotic stability of $P$ because we do not know the asymptotic
behavior of $p(t)$. However, the previous lemmas give us the asymptotics of $y,w,z$,
and we state them in the following corollary.

\begin{corollary}
 Assume the hypotheses of Proposition \ref{prop:P}. If
 \(Q = (y_0, w_0, z_0, p_0) \in I_{\varepsilon^2}\), then for the trajectory
 \(X(t,Q) = (y(t), w(t), z(t), p(t))\) it holds
 $$
 \lim_{t\to+\infty} y(t)
 = \lim_{t\to+\infty} w(t)
 = \lim_{t\to+\infty} z(t) = 0,
 $$
and $ y(t), w(t), z(t)$ are decreasing functions with exponential decay.
 \end{corollary}


 \subsection{Case 2}
 We now deal with the equilibria $P= (0,\,0,\, \tilde{z},\,0)$, with
 \(\tilde{z}\ge 0\).  In this case the Jacobian matrix is
 \begin{equation} \label{eq:jacobianaB}
 J_F(P)= \begin{pmatrix}
 \gamma \lambda & 0 & 0 & 0 \\
 0 & \delta_2 \tilde{z} - 1 & 0 & 0 \\
 0 & (k_1 - \delta_2) \tilde{z} & 0 & -\gamma_2 \tilde{z} \\
 0 & 0 & 0 & 0
 \end{pmatrix}.
 \end{equation}
The eigenvalues are
 \[
 \rho_1 = \gamma \lambda > 0,\quad
 \rho_2 = \delta_2 \tilde{z} - 1,\quad
 \rho_3 = \rho_4 = 0.
 \]
There is always a positive eigenvalue, hence every equilibrium in this family is
unstable.

\begin{proposition}
 All the equilibria $P= (0,\,0,\, \tilde{z},\,0)$, with \(\tilde{z}\ge 0\), are
 unstable.
\end{proposition}


\subsection{Case 3}
We now deal with the equilibria \( P=( \lambda, \, 0, \, \tilde{z}, \, 0)\), with
$\tilde{z} \geq 0$. At these point the Jacobian matrix is
 \begin{equation} \label{eq:jacobianaC}
 J_F(P)=\begin{pmatrix}
 -\gamma \lambda & -\delta_1 \lambda & 0 & -\gamma_1 \lambda \\
 0 & \delta_1 \lambda + \delta_2 \tilde{z} - 1 & 0 & 0 \\
 0 & (k_1 - \delta_2) \tilde{z} & 0 & -\gamma_2 \tilde{z}\\
 0 & 0 & 0 & -\gamma_4 \lambda
 \end{pmatrix}.
 \end{equation}
The eigenvalues are
 \[
 \rho_1 = -\gamma \lambda, \quad
 \rho_2 = \delta_1 \lambda + \delta_2 \tilde{z} - 1, \quad
 \rho_3 = 0, \quad
 \rho_4 = -\gamma_4 \lambda.
 \]
 We obtain an obvious instability result.

\begin{proposition}
 If \(\delta_1 \lambda + \delta_2 \tilde{z} - 1 > 0\), the equilibrium
 \(P= (\lambda, \, 0, \, \tilde{z}, \, 0)\) is unstable.
\end{proposition}


We now deal with the case
\begin{equation} \label{eq:A3}
 \delta_1 \lambda + \delta_2 \tilde{z} - 1 < 0.
\end{equation}
  We are going to prove the stability of equilibria in this case, so from now on we
assume that \eqref{eq:A3} holds. As for Case 1, we will prove the result through a
series of lemmas. First of all, from \eqref{eq:A3} we obtain that it is possible to
fix
$$
 \varepsilon_1 = \frac{1 - \delta_1 \lambda
 - \delta_2 \tilde{z}}{1 + \delta_1 + \delta_2}>0,
$$
so that for all $\varepsilon \in (0, \varepsilon_1 )$ we have
\begin{equation}  \label{eq:A3eps}
 \delta_1 (\lambda +\varepsilon )+ \delta_2 (\tilde{z} + \varepsilon ) - 1 + \varepsilon< 0.
\end{equation}
Notice that $\varepsilon_1 <1$.
As above, we will study the behavior of a trajectory
$X(t)= (y(t), w(t), z(t), p(t))$ with $X(0)= (y_0, w_0 , z_0 , p_0 )$.


\begin{lemma}  \label{le:wa3}
Assume $0< \varepsilon < \min \{ \varepsilon_1 , \lambda \} $,
$y_0 \in (\lambda - \varepsilon^2 , \lambda +\varepsilon^2 )$,
$z_0 \in (\tilde{z}- \varepsilon^2 , \tilde{z} +\varepsilon^2 ) $,
$0\leq w_0 < \varepsilon^3$. Then there is $b_1 \in (0, +\infty]$ such that
$$
 w(t) \leq \varepsilon^3 \exp({-\varepsilon t})
$$
for all $t \in [0, b_1 )$.
\end{lemma}

\begin{proof}
If $w_0 =0$ then $w(t)=0$ for all $t \geq 0$ so the thesis holds.
Assume $w_0>0$. From the equation for \(w\) we obtain
\begin{align*}
 w'(0)
 &= \delta_1 w_0 y_0 + \delta_2 w_0 z_0 - w_0
 = w_0(\delta_1 y_0 + \delta_2 z_0 - 1)\\
 &\le w_0(\delta_1(\lambda + \varepsilon^2) + \delta_2(\tilde{z}+ \varepsilon^2) - 1)\\
 &< w_0(\delta_1(\lambda + \varepsilon) + \delta_2(\tilde{z} + \varepsilon) - 1)
 < -\varepsilon w(0).
 \end{align*}
 Let us define
\[
 b_1 = \sup\{t > 0 \mid \delta_1 y(s) + \delta_2 z(s) - 1 + \varepsilon < 0,\ \forall s \in [0,t)\}.
\]
By continuity we have \(b_1 \in (0,+\infty]\), and also
\(w'(t) < -\varepsilon w(t)\) for all $t\in [0, b_1 )$, hence
\[
 w(t) \le w_0 e^{-\varepsilon t}, \quad t \in [0,b_1).
\]
By the hypothesis on $w_0$ we have
\[
 w(t) \le \varepsilon^3 e^{-\varepsilon t}, \quad t \in [0,b_1).
\]
 \end{proof}

The next lemma deals with the behavior of $z(t)$.

 \begin{lemma} \label{le:za3}
There is $\varepsilon_2 < \varepsilon_1$ such that, assuming
$\varepsilon < \varepsilon_2$ and the hypotheses of Lemma \ref{le:wa3}, it is
$$
 z(t) \leq \tilde{z} + \left( k_1 \tilde{z} +1 \right)\varepsilon^2
 +\BigO(\varepsilon^4) < \tilde{z} + \varepsilon
$$
for all $t \in [0 ,b_1 )$, where $b_1$ is defined in Lemma \ref{le:wa3}.
\end{lemma}

\begin{proof}
For $t \in [0 ,b_1 )$ we have, from the equation for $z$ and Lemma \ref{le:wa3},
 \[
 z'(t) \le k_1 \frac{w(t)}{w(t)+1} z(t) \le k_1 w (t)z (t)
 \le k_1 \varepsilon^3 e^{-\varepsilon t} z(t), \quad t \in [0,b_1).
 \]
Hence
\[
 \frac{z' (t)}{z(t)} \le k_1 \varepsilon^3 e^{-\varepsilon t},
\]
and an integration gives
\[
 \log \frac{z(t)}{z_0}
 \le k_1 \varepsilon^2(1 - e^{-\varepsilon t})
 \le k_1 \varepsilon^2,
\]
 thus
\[
 z(t) \le z_0 \exp(k_1 \varepsilon^2), \quad t \in [0,b_1).
\]
Now we have
\[
 z_0 \exp(k_1 \varepsilon^2)
 \le (\tilde{z} + \varepsilon^2)(1 + k_1 \varepsilon^2 + \BigO(\varepsilon^4))
 = \tilde{z} + \left( k_1 \tilde{z} +1 \right)\varepsilon^2 +\BigO(\varepsilon^4) .
\]
We now fix $\varepsilon_2 < \varepsilon_1$ such that, for all
$\varepsilon< \varepsilon_2$, it holds
$ \left( k_1 \tilde{z} +1 \right)\varepsilon^2 +\BigO(\varepsilon^4) <\varepsilon$,
and we have the conclusion.
\end{proof}

In the next lemma we show that \(b_1 = +\infty\).

 \begin{lemma} \label{lem:b1a3}
Assume the hypotheses of Lemmas \ref{le:wa3} and \ref{le:za3}.
There is a positive $\varepsilon_3 <\varepsilon_2$ such that, if
$ \varepsilon <\varepsilon_3$, then \(b_1 = +\infty\).
 \end{lemma}

\begin{proof}
We first notice that from Proposition \ref{prop:esty}, as we are assuming
 \(y_0 \in (\lambda - \varepsilon^2,\, \lambda + \varepsilon^2)\),
we easily obtain
\[
 y(t) \le \lambda + \varepsilon^2 < \lambda + \varepsilon
 \quad \text{for every } t \ge 0.
 \]

Now assume by contradiction that \(b_1 < +\infty\).
By continuity and the definition of $b_1$ we have
\[
 0 = \delta_1 y(b_1) + \delta_2 z(b_1) - 1 + \varepsilon.
\]
On the other hand, we know that
 $z(b_1) \le \tilde{z} + \left( k_1 \tilde{z} +1 \right)\varepsilon^2
+\BigO(\varepsilon^4)$ and $y(b_1) \le \lambda + \varepsilon^2$, hence
\begin{align*}
 0 &\le \delta_1(\lambda + \varepsilon^2)
 + \delta_2 \left[ \tilde{z} + \left( k_1 \tilde{z} +1 \right)\varepsilon^2
 +\BigO(\varepsilon^4) \right]
 -1 + \varepsilon\\
& \leq \delta_1 \lambda + \delta_2 \tilde{z} - 1 + \varepsilon
 + (\delta_1 + \delta_2 )\varepsilon^2 +\BigO(\varepsilon^4)
 < -(\delta_1 + \delta_2 )\varepsilon
 +\delta_2 )\varepsilon^2 +\BigO(\varepsilon^4)
\end{align*}
 because, from \eqref{eq:A3eps},
 \(\delta_1 \lambda + \delta_2 \tilde{z} - 1 + \varepsilon < -(\delta_1
 + \delta_2)\varepsilon\).
  Thus we obtain
\[
 (\delta_1 + \delta_2 )\varepsilon
 < \delta_2 )\varepsilon^2 +\BigO(\varepsilon^4),
\]
that is,
\begin{equation}
 \label{eq:deltas}
 1< \varepsilon +\BigO(\varepsilon^3),
\end{equation}
which of course is false for small $\varepsilon$'s. We can then state that there is
$\varepsilon_3 < \varepsilon_2$ such that, if $\varepsilon<\varepsilon_3$,
then \eqref{eq:deltas} does not hold, which implies that $b_1 <+\infty$ is impossible.
Therefore \(b_1 = +\infty\).
 \end{proof}

In the next lemma we study the behavior of $p(t)$.

 \begin{lemma}\label{le:pa3}
There are $\varepsilon_4 <\varepsilon_3$ and $b_2 \in (0, +\infty ]$ such that, if
$ \varepsilon \in (0, \varepsilon_4 )$, $0 \leq p_0 < \varepsilon^3$ and
$y_0, w_0, z_0$ are chosen as in the previous lemmas, then
$$
 p(t) \leq \varepsilon^3 \exp \Big(-\frac{\lambda \gamma_4}{2} t \Big), \quad
 \forall t \in [0, b_2 ).
$$
\end{lemma}

 \begin{proof}
If $p_0=0$ then $p(t)=0$ for all $t\geq 0$ and the thesis holds. So we assume $p_0>0$.
From the equation for $p(t)$ we have
\[
 p'(t) \le (\gamma_3 w(t) - \gamma_4 y(t)) p(t),
\]
and therefore
\[
 p'(0) \le \bigl(\gamma_3 \varepsilon^3 - \gamma_4(\lambda - \varepsilon^2)\bigr)p_0.
\]
Now we choose \(\varepsilon_4 <\varepsilon_3\) such that, for all
$\varepsilon<\varepsilon_4$, it holds
\[
 \gamma_3 \varepsilon^3 + \gamma_4 \varepsilon^2
 < \frac{\lambda \gamma_4}{2}.
\]
which means
\[
 \gamma_3 \varepsilon^3 - \gamma_4(\lambda - \varepsilon^2)
 < -\frac{\lambda \gamma_4}{2}.
\]
We then have
\[
 p'(0) < -\frac{\lambda \gamma_4}{2} p(0)
\]
and we define
\[
 b_2 = \sup\big\{t > 0 \mid p'(s) < -\frac{\lambda \gamma_4}{2} p(s),\; \forall s
 \in [0,t)\big\}.
\]
By continuity, \(b_2 > 0\), and arguing as before we obtain
 \[
 p(t) \le p_0 \exp{\Big( - \frac{\lambda \gamma_4}{2}t \Big)}
 < \varepsilon^3 \exp{\Big( - \frac{\lambda \gamma_4}{2}t \Big)} \quad 
 \text{for every } t \in [0,b_2).
 \]
 \end{proof}

We now want to study the behavior of \(y(t)\) on \([0,b_2)\). We have already
noticed in Lemma \ref{lem:b1a3} that, if $y_0 < \lambda + \varepsilon^2$, then
$y (t) < \lambda + \varepsilon^2$ for all $t\geq 0$. In the following lemma we prove
the inequality $y (t) > \lambda - \varepsilon^2$.

\begin{lemma}
There is $\varepsilon_5 \leq \varepsilon_4$ such that, if
$\varepsilon \in (0, \varepsilon_5 )$ and $y_0 ,w_0 ,z_0 , p_0$ are chosen as in the
previous lemmas, then
$$
 y(t) > \lambda - \varepsilon^2, \quad \forall t \in [0, b_2 ).
$$
\end{lemma}

\begin{proof}
From the equation for $y(t)$ and the previous results we obtain, for all
 $ \varepsilon < \varepsilon_4 $ and all $t \in [0, b_2 )$,
\[
 y'(t)
 = y(\gamma \lambda - \gamma y - \delta_1 w - \gamma_1 p)
 \ge y\bigl(\gamma \lambda - \gamma y - \delta_1 \varepsilon^3 - \gamma_1 \varepsilon^3\bigr).
\]
Choosing
$$
 \varepsilon_5 = \min \Big\{ \varepsilon_4 , \frac{\gamma}{\delta_1 + \gamma_1}
  \Big\},
$$
we have that for all $\varepsilon \in (0, \varepsilon_5 )$ it holds
\[
 (\delta_1 + \gamma_1)\varepsilon^3 < \gamma \varepsilon^2.
\]
We then obtain
\[
 y'(t) > y(t)\bigl(\gamma \lambda - \gamma y(t) - \gamma \varepsilon^2\bigr)
 = \gamma y(t)(\lambda_1 - y(t)),
\]
where \(\lambda_1 = \lambda - \varepsilon^2\).

Let \(\overline{y}(t)\) be the solution of
\begin{gather*}
 \overline{y}'(t) = \gamma \overline{y}(t)\bigl(\lambda_1 - \overline{y}(t)\bigr),\\
 \overline{y}(0) = y_0.
\end{gather*}
By comparison it is easy to obtain \(y(t) > \overline{y}(t)\) for every
\(t \in [0,b_2)\). As $y_0 > \lambda_1$ by hypothesis, from the standard properties
of the logistic equation we obtain
 \(\overline{y}(t) > \lambda_1\) for every \(t\), hence
\[
 y(t) > \lambda_1 = \lambda - \varepsilon^2
 \quad \text{for every } t \in [0,b_2).
\]
\end{proof}

 We now prove that $b_2 = +\infty$.


\begin{lemma}
Assume $\varepsilon < \varepsilon_5 $ and $y_0 , w_0 , z_0 , p_0$ as in the hypotheses
of the previous lemmas. Then $b_2 = +\infty$.
 \end{lemma}

\begin{proof}
Assume by contradiction that \(b_2 < +\infty\). From the previous lemma we then
obtain \(y(b_2) \ge \lambda - \varepsilon^2\), and arguing as in Lemma \ref{le:pa3}
we obtain
\[
 p'(b_2)
 \le \bigl(\gamma_3 \varepsilon^3 - \gamma_4(\lambda - \varepsilon^2)\bigr)p(b_2)
 < -\frac{\lambda \gamma_4}{2} p(b_2),
\]
so that, by continuity, \(p'(s) < -\frac{\lambda \gamma_4}{2} p(s)\) on some interval
\((b_2 - \delta, b_2 + \delta)\), contradicting the definition of \(b_2\).
The hypothesis \(b_2 < +\infty\) gives rise to a contradiction, so $b_2 = +\infty$.
\end{proof}

In the next lemma we complete the study of \(z(t)\).

\begin{lemma}
There is $\varepsilon_6\leq \varepsilon_5$ such that, if
$\varepsilon \in (0, \varepsilon_6 )$ and we choose $y_0 , w_0 , z_0 , p_0$ as in the
hypotheses of the previous lemmas, we obtain $z(t) > \tilde{z} - \varepsilon$ for all
$t\geq 0$.
 \end{lemma}

\begin{proof}
From the equation for $z(t)$ we obtain
 \[
 z'(t)
 \ge (-\delta_2 w - \gamma_2 p)\,z
 \ge \bigl(-\delta_2 \varepsilon^3 e^{-\varepsilon t}
 - \gamma_2 \varepsilon^3 e^{-a t}\bigr) z(t),
 \]
where $a= \frac{\lambda \gamma_4}{2} $. Thanks to the previous lemmas, this
inequality holds for all $t\geq 0$. Hence, dividing by \(z(t)\) and integrating,
\[
 \log \frac{z(t)}{z_0}
 \ge \delta_2 \varepsilon^2 (e^{-\varepsilon t} - 1)
 + \gamma_2 \frac{\varepsilon^2}{a}(e^{-a t} - 1)
 \ge -\delta_2 \varepsilon^2 - \frac{\gamma_2 \varepsilon^3}{a},
\]
and therefore
\[
 z(t) \ge z_0 \exp\!\Big(-\delta_2 \varepsilon^2 - \frac{\gamma_2
 \varepsilon^3}{a}\Big),
 \quad t \ge 0.
\]
We have assumed \(z_0 \in (\tilde{z} - \varepsilon^2,\, \tilde{z} + \varepsilon^2)\).
Then
 \[
 z(t)
 \ge (\tilde{z} - \varepsilon^2)
 \exp\!\Big(-\delta_2 \varepsilon^2 - \gamma_2 \frac{\varepsilon^3}{a}\Big)
 = (\tilde{z} - \varepsilon^2)
 \Big(1 - \delta_2 \varepsilon^2 - \gamma_2 \frac{\varepsilon^3}{a}
 + \BigO(\varepsilon^4)\Big)
 = \tilde{z} + \BigO(\varepsilon^2).
 \]
We can then choose \(\varepsilon_6 < \varepsilon_5\) such that, for all
$\varepsilon \in (0, \varepsilon_6 )$ it is
\(\BigO(\varepsilon^2) >- \varepsilon\),
which in turn implies \(z(t) \ge \tilde{z}- \varepsilon\) for all $t \geq 0$.
\end{proof}

The previous lemmas imply the Lyapunov-stability of the equilibrium
$P= ( \lambda , 0, \tilde{z} , 0 )$ in the hypothesis
$\delta_1 \lambda + \delta_ 2 \tilde{z} <1$. To state the result, we consider a
neighborhood basis for $P$ given by
 \[
 I_{\eta}(P) = (\lambda - \eta,\, \lambda + \eta) \times (-\eta, \eta)
 \times (\tilde{z} - \eta,\, \tilde{z} + \eta) \times (-\eta, \eta),
 \]
 with \(0< \eta < \lambda\).
The following proposition, which states the Lyapunov-stability of $P$, is an easy
consequence of the previous lemmas.

 \begin{proposition}  \label{prop:lyapa3}
Let $P= ( \lambda , 0, \tilde{z} , 0 )$ and assume
$\delta_1 \lambda + \delta_ 2 \tilde{z} <1$. Then $P$ is a Lyapunov-stable equilibrium,
as follows. Fix $\eta \in (0, \lambda )$ and
$\varepsilon <\min \{ \eta, \varepsilon_6 \} $. If \(Q \in I_{\varepsilon^3}(P)\) and
$X(t,Q)$ is the trajectory starting at $Q$, then \(X(t,Q) \in I_{\eta}(P)\) for every
\(t \ge 0\).
 \end{proposition}

We are not able to prove the asymptotic stability of the equilibria
$P= ( \lambda , 0, \tilde{z}, 0 )$, because we do not know if, for $Q $ near to $P$,
it holds that $z(t) \to \tilde{z}$ as $t \to +\infty$. However, we have precise
information on the asymptotic behavior of the other variables.
Regarding $w(t), p(t)$, the previous lemmas imply trivially the following corollary.

\begin{corollary}
Under the  hypotheses  in Proposition \ref{prop:lyapa3}, let
\(X(t,Q)=(y(t),w(t),z(t),p(t))\). Then \(w(t)\to0\) and \(p(t)\to0\) as
\(t\to+\infty\).
\end{corollary}

 Regarding the asymptotic behavior of $y(t)$, some more work is needed.

\begin{proposition}
 Assume the same hypotheses as in Proposition \ref{prop:lyapa3}, and let
$X(t, Q)= (y(t), w(t), z(t), p(t))$. Then $y(t) \to \lambda$ as $t \to +\infty$.
\end{proposition}

\begin{proof}
Recall that we have shown that
\(y(t) \in (\lambda - \varepsilon,\, \lambda + \varepsilon)\) for all
\(t\). To prove that \(\lim_{t \to +\infty} y(t) = \lambda\), we study separately
the \(\limsup\) and \(\liminf\).

 First, we show that \(\limsup_{t \to + \infty} y(t) \le \lambda\).
 We consider two cases:
\begin{itemize}
 \item[(i)] If $ y_0 \leq \lambda$, by Proposition \ref{prop:esty} we know that
 $y(t) \leq \lambda $ for all $t$, hence $\limsup_{t \to + \infty} y(t) \le \lambda $.

 \item[(ii)] If $y_0 > \lambda $, assume first that there exists $t_0 > 0$ such that
 $y(t_0) < \lambda$: in this case it is easy to obtain $y(t) < \lambda$ for
 $t > t_0$, so that $\limsup_{t \to + \infty} y(t) \le \lambda $. If instead
 \(y(t) \geq \lambda\) for all \(t>0\), then \(y'(t) \leq 0\) and therefore
 \(y\) is decreasing and is bounded from below. Thus \(\lim_{t \to + \infty} y(t)
 = L \ge \lambda\). Hence, using the equation for $y$ and the previous lemmas,
 we obtain \(\lim_{t \to + \infty} y'(t) = \gamma L \left( \lambda -L \right) \).
 By standard results we know that it must be \(\lim_{t \to + \infty} y'(t) = 0\),
 which implies \( \lambda =L =\lim_{t \to + \infty} y(t)\), and therefore
 \(\limsup_{t \to + \infty} y(t) \le \lambda\).
\end{itemize}
Hence, in both cases,
\[
 \limsup_{t \to + \infty} y(t) \le \lambda.
\]
We now study \(\liminf_{t \to + \infty} y(t)\) and show that it is greater than or equal to \(\lambda\).
Fix any \(\sigma \in (0, \gamma \lambda)\). Since \(w(t), p(t) \to 0\) as \(t \to + \infty\), there exists \(t_{\sigma} >0\) such that \(\delta_1 w + \gamma_1 p < \sigma\) for all \(t \ge t_{\sigma}\). Then, for \(t \ge t_{\sigma}\),
\[
 y' = y(\gamma \lambda - \gamma y - \delta_1 w - \gamma_1 p) > y(\gamma \lambda - \gamma y - \sigma)
 = \gamma y \Big(\lambda - \frac{\sigma}{\gamma} - y\Big) 
 = \gamma y (\lambda_1 - y)
\]
where \(\lambda_1 = \lambda - \frac{\sigma}{\gamma} > 0\) by hypothesis.

Define \(\overline{y}_{\sigma} = \frac 12 \min \{y(t_{\sigma}), \lambda_1\} > 0\)
and let \(\overline{y}\) be the solution to
\begin{gather*}
  \overline{y}'(t) = \gamma \overline{y}(\lambda_1 - \overline{y}) \\
  \overline{y}(t_{\sigma}) = \overline{y}_{\sigma}
\end{gather*}
 Since \(\overline{y}(t_{\sigma}) = \overline{y}_{\sigma} < y(t_{\sigma})\),
 we obtain \(y(t) > \overline{y}(t)\) for all \(t \ge t_{\sigma}\).
Also, \(0 < \overline{y}_{\sigma} < \lambda_1\), and therefore by standard results
on logistic equation \(\overline{y}(t) < \lambda_1\) for all \(t \ge t_{\sigma}\) and
\(\lim_{t \to + \infty} \overline{y}(t) = \lambda_1\). Hence,
\[
 \liminf_{t \to + \infty} y(t) \ge \lim_{t \to + \infty} \overline{y}(t)
 = \lambda_1 = \lambda - \frac{\sigma}{\gamma}.
\]
Since this holds for all \(\sigma \in (0, \gamma\lambda)\),
letting \(\sigma \to 0^+\) yields
\[
 \liminf_{t \to + \infty} y(t) \ge \lambda.
\]
Thus we have obtained
\[
 \liminf_{t \to + \infty} y(t) \ge \lambda \ge \limsup_{t \to + \infty} y(t)
\]
and therefore
\begin{equation*}
 \liminf_{t \to + \infty} y(t)
 = \limsup_{t \to + \infty} y(t)
 = \lim_{t \to + \infty} y(t)
 =\lambda. \qedhere
\end{equation*}
 \end{proof}

 \section{Isolated equilibria}

In this section we study the stability of the isolated equilibria that we have
computed above. Because  many parameters are involved, we will deal, in most cases,
with the simplified model \eqref{eq:first_simply}, proving only asymptotic results.

 \subsection{Case 4}
 The equilibrium is
 \[
 P_0 = \left(0,\,\frac{k_1}{\delta_2} - 1,\,\frac{1}{\delta_2},\,0\right),
 \]
 with the condition \(k_1 \geq \delta_2\) required to be in ${\overline{\mathcal{C}}}$. Notice that \(k_1 = \delta_2\) implies the $P_0$ is a point on the half-line of Case 2) so for the study of $P_0$ we assume \(k_1 > \delta_2\).
 \par \noindent The Jacobian matrix for \(P_0\) is
 \begin{equation} \label{eq:jacobianaP2}
 J_F(P_0)=
 \begin{pmatrix}
 \gamma \lambda + \delta_1 - \frac{\delta_1}{\delta_2} k_1 & 0 & 0 & 0 \\
 \frac{\delta_1}{\delta_2} k_1 - \delta_1 & 0 & k_1 - \delta_2 & 0 \\
 0 & -1 + \frac{\delta_2}{k_1} & 0 & -\frac{\gamma_2}{\delta_2} \\
 0 & 0 & 0 & j_{4,4}
 \end{pmatrix}
 \end{equation}
 where
 \[
 j_{4,4}
 = \frac{\gamma_3}{\delta_2} k_1 - \gamma_3 - k_2 + \frac{k_2}{k_1} \delta_2
 = \Big( \frac{k_1}{\delta_2} - 1 \Big) \Big( \gamma_3 - \frac{k_2}{k_1} \delta_2 \Big).
 \]
The eigenvalues are
 \[
 \rho_1 = \gamma \lambda + \delta_1 - \frac{\delta_1}{\delta_2} k_1, \quad
 \rho_{2,3} = \pm i \frac{|k_1 - \delta_2|}{\sqrt{k_1}}, \quad
 \rho_4 = j_{4,4}.
 \]
From these formulas for the eigenvalues we easily derive the following asymptotic
results.

 \begin{proposition}
  \begin{itemize}
 \item[(i)] If $\delta_1 \to 0^+$ and the other parameters are kept fixed,
 then $\rho_1 \to \gamma \lambda >0$, so $P_0$ is unstable for small $\delta$'s.

 \item[(ii)] If $\delta_2 /k_1 \to 0^+$ and the other parameters are kept fixed,
 then $\rho_4 = j_{4,4} \to +\infty$, so $P_0$ is unstable for small values
 of $\delta_2 /k_1$.

 \item[(iii)] If $\gamma_3 > k_2$ then $\rho_4 >0$ (recall that we are assuming
 $k_1 > \delta_2 $), so $P_0$ is unstable. This holds in particular if $k_2 \to 0^+$
 and the other parameters are kept fixed.
 \end{itemize}
 \end{proposition}

\subsection{Case 5}
Because of the difficulties in dealing with many parameters, from now on we work
with the simplified model \eqref{eq:first_simply}, and it is easy to check that the
Jacobian matrix becomes
 \begin{equation} \label{eq:jacobiana_ipersemplificata}
 J_F=\begin{pmatrix}
 50 - y - \delta w - \alpha p & -\delta y & 0 & -\alpha y \\
 \delta w & \delta y + \delta z - 1 & \delta w & 0 \\
 0 & -\delta z + k \frac{z}{(w+1)^2} & -\delta w + k \frac{w}{w+1} - \frac{\alpha}{2}p & -\frac{\alpha}{2} z \\
 -2 \alpha p & 2 \alpha p - k \frac{p}{(w+1)^2} & 0 & 2 \alpha w - k \frac{w}{w+1} - 2 \alpha y
 \end{pmatrix}
 \end{equation}
In the simplified model the equilibrium is
 \[
 P_0 = \Big( 0,\,
 \frac{k}{2\alpha} - 1,\,
 \frac{1}{\delta},\,
 \frac{1}{\alpha^2} (k - 2\alpha)(2\alpha - \delta)
 \Big),
 \]
with \(P_0 \in \overline{\mathcal{C}}\) if \(k \ge 2\alpha \ge \delta\).
If $k= 2 \alpha$ then $P_0$ reduces to a point of the half-line of Case 2), so it
is unstable. So let us study the case $k > 2 \alpha$.

 The Jacobian matrix evaluated at \(P_0\) becomes
 \begin{equation} \label{eq:jacobiana_ipersemplificataP0}
 J_F(P_0)=
 \begin{pmatrix}
  50 - \frac{1}{2\alpha}(k - 2\alpha)(4\alpha - \delta) & 0 & 0 & 0 \\[4pt]
  \frac{\delta}{2\alpha}(k - 2\alpha) & 0 &  \frac{\delta}{2\alpha}(k - 2\alpha)
  & 0 \\[4pt]
 0 &  \frac{4\alpha^2 - \delta k}{k\delta} & 0 &  -\frac{\alpha}{2\delta} \\[4pt]
  \frac{2}{\alpha}(k - 2\alpha)(\delta - 2\alpha)
   &  \frac{2}{\alpha k}(k - 2\alpha)^2(2\alpha - \delta) & 0 & 0
 \end{pmatrix}.
 \end{equation}
We have the eigenvalue
\[
 \rho_1 = 50 - \frac{1}{\alpha}(k - 2\alpha)
 \Big(2\alpha - \frac{\delta}{2}\Big).
\]
If $\rho_1 >0$, that is
$50 > \frac{1}{\alpha}(k - 2\alpha)\left(2\alpha - \frac{\delta}{2}\right) $, then
$P_0$ is unstable. If $\rho_1 \leq 0$, we need to analyze the reduced
\(3 \times 3\) matrix
 \[
 J_1(P_0) =
 \begin{pmatrix}
 0 &  \frac{\delta}{2\alpha}(k - 2\alpha) & 0 \\[4pt]
  \frac{4\alpha^2 - \delta k}{k\delta} & 0 &  -\frac{\alpha}{2\delta} \\[4pt]
  \frac{2}{\alpha k}(k - 2\alpha)^2(2\alpha - \delta) & 0 & 0
 \end{pmatrix}.
 \]
 The characteristic polynomial \(p_3(\rho) \) of \(J_1(P_0)\) is
 \[
 p_3(\rho)
 = \rho^3
 + \frac{1}{2}\frac{(k\delta - 4\alpha^2)(k - 2\alpha)}{\alpha k} \rho
 + \frac{1}{2}\frac{(k - 2\alpha)^3(2\alpha - \delta)}{\alpha k}.
 \]
Calling $\rho_2 , \rho_3 , \rho_4$ the roots of $p_3$, standard algebraic properties
give that
\[
 \rho_2 + \rho_3 + \rho_4 = 0,\quad
 \rho_2 \rho_3 \rho_4 = -\frac{1}{2}\frac{(k - 2\alpha)^3(2\alpha - \delta)}{\alpha k}.
 \]

If $2 \alpha =\delta$ we are in a degenerate case in which the Jacobian matrix does not
give information on the stability, and we are not able to state a result. If instead
it holds $2 \alpha >\delta$, we can argue as follows. We first notice that
$ \rho_2 \rho_3 \rho_4 <0$. Hence, if all three roots were real, since their sum is
zero, some of them would be positive. On the other hand, in the case of one real root,
say $\rho_2$, and two complex conjugate roots \(\rho_{3,4} = a \pm i b\), we have
$$
 \rho_2 \rho_3 \rho_4 = \rho_2 \left( a^2 + b^2 \right) <0.
$$
This implies $\rho_2 <0$, and also
 \[
 \rho_2 + \rho_3 + \rho_4 = \rho_2 + 2a = 0
 \implies  a = -\frac{\rho_2}{2}>0,
 \]
Hence, in any case, when $k > 2 \alpha$, we obtain at least an eigenvalue with a
strictly positive real part. This means that we have proved the following proposition.

\begin{proposition}\label{prop:8}
The equilibrium $P_0$ of the simplified model \eqref{eq:first_simply} is unstable
when
 \begin{itemize}
 \item[(i)] $\frac{1}{\alpha}(k - 2\alpha)\left(2\alpha
 - \frac{\delta}{2}\right)< 50$.

 \item[(ii)] $\frac{1}{\alpha}(k - 2\alpha)\left(2\alpha
 - \frac{\delta}{2}\right)\geq 50$ and $ 2 \alpha > \delta$.
\end{itemize}
\end{proposition}

\begin{remark} \rm
 We have proved Proposition \ref{prop:8} in the simplified case of system
 \eqref{eq:first_simply}, where several hypotheses are imposed in order to
 reduce the number of parameters. However, since the eigenvalues of a matrix
 depend continuously on its entries, and in our case the entries depend
 continuously on the parameters, small variations in the latter do not affect
 the conclusions of Proposition \ref{prop:8}.
 For example, in deriving system \eqref{eq:first_simply} we assumed
 \(\delta_1 = \delta_2\) and \(k_1 = k_2\). Nevertheless, the results of
 Proposition \ref{prop:8} remain valid also when \(\delta_1 \neq \delta_2\)
 and \(k_1 \neq k_2\), provided that \(|\delta_1-\delta_2|\) and
 \(|k_1-k_2|\) are sufficiently small. The same remark applies to the results
 of the following subsections.
 \end{remark}


\subsection{Case 6}

 The equilibrium point to be studied, in the simplified model, is
 \[
 P_0=\Big( \frac{1}{\delta},\; \frac{1}{2\delta^2}(100\delta - 1),\; 0,\; 0 \Big).
 \]
If \(\delta < \frac{1}{100}\) we obtain
 \(P_0 \notin \overline{\mathcal{C}}\), so we assume \(\delta > \frac{1}{100}\).
Notice that if \(\delta = \frac{1}{100}= 1/\lambda\) then we come back to Case $3)$.
The Jacobian matrix \(J_F(P_0)\) becomes
 \[
 \begin{pmatrix}
  -\frac{1}{2\delta} & -1 & 0 &  -\frac{\alpha}{\delta} \\[4pt]
  50 - \frac{1}{2\delta} & 0 &  50 - \frac{1}{2\delta} & 0 \\[4pt]
 0 & 0 &
  \frac{(1 - 100\delta)\bigl(2\delta^2 - 2k\delta + 100\delta - 1\bigr)}
 {2\delta\bigl(2\delta^2 + 100\delta - 1\bigr)}
 & 0 \\[4pt]
 0 & 0 & 0 &
\frac{196\alpha \delta^3 - 100\delta^3 k + 9798\alpha \delta^2 + \delta^2 k - 198\alpha \delta + \alpha}
 {\delta^2\bigl(2\delta^2 + 100\delta - 1\bigr)}
 \end{pmatrix}.
 \]
The eigenvalues are
\begin{gather*}
  \rho_{1,2} = \frac{-1 \pm \sqrt{1 + 8\delta - 800\delta^2}}{4\delta}, \\
  \rho_3 = \frac{(1 - 100\delta)\bigl(2\delta^2 - 2k\delta + 100\delta - 1\bigr)}
 {2\delta\bigl(2\delta^2 + 100\delta - 1\bigr)}
 = -50 + \frac{1}{2\delta} + k\,\frac{100\delta - 1}{2\delta^2 + 100\delta - 1}, \\
\begin{aligned}
 \rho_4 &=  \frac{196\alpha \delta^3 - 100\delta^3 k + 9798\alpha \delta^2
  + \delta^2 k - 198\alpha \delta + \alpha}
  {\delta^2\bigl(2\delta^2 + 100\delta - 1\bigr)} \\
 & = \frac{\alpha}{\delta^2}(100\delta - 1)
 - k\,\frac{100\delta - 1}{2\delta^2 + 100\delta - 1}
 - \frac{2\alpha}{\delta}.
 \end{aligned}
 \end{gather*}
As \(\delta > \frac{1}{100}\) we have $8\delta - 800 \delta^2 <0$, which implies
that the eigenvalues \(\rho_1\) and \(\rho_2\) are either real and strictly negative,
or a pair of complex conjugates with negative real part.

If \(k \to +\infty\), then \(\rho_3 > 0\), and the point is unstable.
If \(k \to 0\), then \(\rho_3 < 0\) (recall that \(\delta > 1/100\)) while
 \(\rho_4\) is close to
 \[
 \rho_4\big|_{k=0} = \frac{\alpha}{\delta} \Big( 98 - \frac{1}{\delta} \Big).
 \]
Therefore, since \(\delta > \frac{1}{100}\), we obtain that if \(k \to 0\) and
 \(\frac{1}{100} < \delta < \frac{1}{98}\), all four eigenvalues have negative real
 part, and the equilibrium is asymptotically stable.
 If \(k \to 0\) and \(\delta > \frac{1}{98}\), then \(\rho_4 > 0\), and the
 equilibrium is unstable.
Let us summarize these results in the following proposition

\begin{proposition}
Let
\[
 P_0=\Big( \frac{1}{\delta},\; \frac{1}{2\delta^2}(100\delta - 1),\; 0,\; 0 \Big).
\]
be the equilibrium in the case 6) of the simplified model \eqref{eq:first_simply},
with \(\frac{1}{100} < \delta\). Then we have
 \begin{itemize}
 \item[(i)] If $k\to +\infty$, $P_0$ is unstable.
 \item[(ii)] If $k\to 0^+$ and \(\delta > \frac{1}{98}\), $P_0$ is unstable.
 \item[(iii)] If $k\to 0^+$ and \(\frac{1}{100} < \delta < \frac{1}{98}\),
 $P_0$ is asymptotically stable.
 \end{itemize}
\end{proposition}

 \subsection{Case 7}

 In this case we have the equilibrium
$ P_0 = \Big( \frac{1}{\delta_1},\, w_0,\, 0,\, p_0 \Big)$,
where $w_0 , p_0$ in the simplified model are
\begin{gather*}
 w_0
 = \frac{2\alpha + k\delta - 2\alpha\delta + \sqrt{4\alpha^2 \delta^2 - 4\alpha \delta^2 k + \delta^2 k^2 + 8\alpha^2 \delta + 4\alpha \delta k + 4\alpha^2}}{4\alpha \delta}, \\
 p_0 = \frac{1}{\alpha}\Big(50 - \frac{1}{2\delta} - \delta w_0\Big).
 \end{gather*}
 We immediately observe that if \(\delta \to 0^+\) then \(p_0 < 0\) and hence
 \(P_0 \notin \overline{\mathcal{C}}\). In the same way, if $k\to +\infty$ then
 $w_0 \to +\infty$ and $p_0 <0$. So we do not consider the cases \(\delta \to 0^+\)
 or $k\to +\infty$. Let us study the case $k\to 0^+$.
The Jacobian matrix \(J_F(P_0)\) becomes
 \[
 \begin{pmatrix}
 50 - \frac{1}{\delta} - \delta w_0 - \alpha p_0 & -1 & 0 & -\frac{\alpha}{\delta} \\
 \delta w_0 & 0 & \delta w_0 & 0 \\
 0 & 0 & -\delta w_0 + k \frac{w_0}{w_0 + 1} - \frac{\alpha}{2} p_0 & 0 \\
 -2\alpha p_0 & 2\alpha p_0 - k \frac{p_0}{(w_0+1)^2} & 0 & 2\alpha w_0
 - k \frac{w_0}{w_0 + 1} - \frac{2 \alpha}{\delta}
 \end{pmatrix}.
 \]
If $k\to 0^+$ then $w_0 \to \frac{1}{\delta}$ and
$p_0 \to \frac{1}{\alpha}\left(49 - \frac{1}{2\delta}\right)$, so a necessary condition
to have the the equilibrium point in $ \overline{\mathcal{C}}$, as $k\to 0^+$, is
\(\delta \geq \frac{1}{98}\). We now compute the limit as $k\to 0^+$ of the matrix
$J_F (P_0 )$, which we call $J_0$:
 \[
 J_0=
 \begin{pmatrix}
 -\frac{1}{2\delta} & -1 & 0 & -\frac{\alpha}{\delta} \\
 1 & 0 & 1 & 0 \\
 0 & 0 & \frac{1}{4\delta}(1 - 102\delta) & 0 \\
 \frac{1}{\delta} - 98 & -\frac{1}{\delta} + 98 & 0 & 0
 \end{pmatrix}.
 \]
 We want to use the well-known fact that the eigenvalues of a matrix are continuous
 functions of the entries, so we study the eigenvalues of $J_0$.
We have an eigenvalue
\[
 \rho_3 = \frac{1}{4\delta}(1 - 102\delta) < 0,
\]
since \(\delta \geq \frac{1}{98}\).
We then consider the reduced matrix
\[
 J_1 =
 \begin{pmatrix}
 -\frac{1}{2\delta} & -1 & -\frac{\alpha}{\delta} \\
 1 & 0 & 0 \\
 \frac{1}{\delta} - 98 & -\frac{1}{\delta} + 98 & 0
 \end{pmatrix}.
 \]
Its characteristic polynomial is
 \[
 p_3(\rho)
 = \rho^3
 + \frac{1}{2\delta}\,\rho^2
 + \frac{\delta^2 - 98\alpha \delta + \alpha}{\delta^2}\,\rho
 + \frac{98\delta - 1}{\delta^2}\,\alpha.
 \]
To apply Routh-Hurwitz criterion we set
\begin{gather*}
 a_1 = \frac{1}{2\delta} > 0 \quad
 a_2 = \frac{\delta^2 - 98\alpha \delta + \alpha}{\delta^2}, \\
 a_3 = \alpha \frac{98\delta - 1}{\delta^2} > 0\quad
  a_1 a_2 - a_3 = \frac{(1 - 196\alpha)\delta^2 - 96\alpha \delta + \alpha}{2\delta^3}.
 \end{gather*}
We notice that \(a_2>0\) and \(a_1 a_2 > a_3\) when $\alpha \to 0^+$.
Hence the eigenvalues of $J_1$ have strictly negative real parts for small $\alpha$'s,
and by continuity all the eigenvalues of $J_F(P_0 )$ have strictly negative real parts
as $k\to 0^+$. So we obtain that $P_0$ is asymptotically stable for small $\alpha$'s
and $k$'s. We state our result in the following proposition, paying attention to the
fact that we have computed limits first for $k \to 0^*$ and the for $\alpha \to 0^+$.

\begin{proposition}
There is $\alpha^* >0$ such that for all $\alpha \in (0, \alpha^* )$ there is
$k(\alpha ) >0$ such that, for all $k \in (0, k(\alpha ))$ the equilibrium
$ P_0 $ is asymptotically stable.
 \end{proposition}

\subsection{Case 8}
 The equilibrium point, in the simplified model, is
$$
 P_0 =(y_0 , w_0 , z_0 , 0)= \Big(100 - 2k + 2\delta,\, \frac{k}{\delta} - 1,\,
  \frac{1}{\delta} - 100 + 2k - 2\delta,\, 0\Big).
$$
Observe that we must have \(k \le 50 + \delta\) in order for \(y_0 \ge 0\), so it is
meaningless to study the case \(k \to +\infty\).
We must also have \(k \ge \delta\) in order for \(w_0 \ge 0\), thus, if $\delta >0$
is fixed, it is also pointless to study \(k \to 0^+\).
The Jacobian matrix becomes
\[
 J_F(P_0)=
 \begin{pmatrix}
 -50 + k - \delta & -\delta y_0 & 0 & -\alpha y_0 \\
 k - \delta & 0 & k - \delta & 0 \\
 0 & - z_0 \frac{\delta}{k}(k-\delta) & 0 & -\frac{\alpha}{2} z_0 \\
 0 & 0 & 0 & w_0(2\alpha - \delta) - 2\alpha y_0
 \end{pmatrix}.
\]
The matrix \(J_F(P_0)\) has the eigenvalue
 \(\rho_1 = w_0(2\alpha - \delta) - 2\alpha y_0\),
 while the characteristic polynomial of the reduced matrix
 \[
 J_1(P_0)=
 \begin{pmatrix}
 -50 + k - \delta & -\delta y_0 & 0 \\
 k - \delta & 0 & k - \delta \\
 0 & - z_0 \frac{\delta}{k}(k-\delta) & 0
 \end{pmatrix}
 \]
 is \(p_3(\rho)= \rho^3 + a_1 \rho^2 + a_2 \rho + a_3\), where
  \begin{gather*}
 a_1 = \tfrac12 y_0, \quad
  a_2 = \frac{\delta}{k} z_0 (\delta-k)^2 + \delta y_0 (k-\delta), \\
 a_3 = \frac{\delta}{2k} y_0 z_0 (\delta-k)^2, \quad
 a_1 a_2 - a_3 = \tfrac12 \delta y_0^2 (k-\delta).
 \end{gather*}
To apply the Routh-Hurwitz criterion, we notice that these quantities are positive
 whenever  \(k>\delta\), a condition already required for \(w_0\ge0\).
 Hence all eigenvalues of \(J_1(P_0)\) have negative real part.
 Therefore, whenever \(\rho_1 \neq 0\), its sign determines the stability of \(P_0\).
 We now study the sign of \(\rho_1\) in a suitable asymptotic regime.
If \(\delta \to 0^+\) we have
 \[
 y_0 \to 100 - 2k, \quad
 w_0 \to +\infty, \quad
 z_0 \to +\infty, \quad
 \delta w_0 \to k.
 \]
 Thus
 \[
 \rho_1 = 2\alpha(w_0 - y_0) -\delta w_0\to +\infty,
 \]
and consequently the equilibrium is unstable.

Let us now see some results dealing with the parameter $\alpha$.
We have
 \[
 \rho_1
 = \Big(\frac{k}{\delta} - 1\Big)(2\alpha - \delta)
 - 2\alpha(100 - 2k + 2\delta)
 = 4\alpha k - 4\alpha \delta - 202\alpha + \delta - k + 2\alpha \frac{k}{\delta}.
 \]
 We notice that
 \(\rho_1 \to \delta - k < 0\) as \( \alpha \to 0\).
 \noindent Hence, in this case we have stability. If $\alpha \to +\infty$, then
  the sign of $\rho_1$ is the same as that of
 $$
 2(k-\delta ) + \frac{k}{\delta}-101.
 $$
We summarize our result in the following proposition.

 \begin{proposition} \label{prop:19}
 \begin{itemize}
 \item[(i)] If $\delta\to 0^+$, the equilibrium $P_0$ is unstable.
 \item[(ii)]If $\alpha \to 0^+$, the equilibrium $P_0$ is asymptotically stable.
 \item[(iii)] If $\alpha \to +\infty$ and $ 2(k-\delta ) + \frac{k}{\delta}<101$,
 the equilibrium $P_0$ is asymptotically stable.
 \item[(iv)] If $\alpha \to +\infty$ and $ 2(k-\delta ) + \frac{k}{\delta}>101$,
 the equilibrium $P_0$ is unstable.
 \end{itemize}
\end{proposition}

\subsection{Case 9}
In the simplified model the equilibria in this case are given by the solution of
the  system
 \begin{gather*}
  50 - \frac{y}{2} - \delta w - \alpha p = 0, \\
  z = \frac{1}{\delta} - y, \\
  p = \frac{2}{\alpha}\Big( -\delta w + k \frac{w}{w+1} \Big), \\
  y = w - \frac{k}{2\alpha} \frac{w}{w+1}.
 \end{gather*}
We obtain two values for \( w \), namely
\[
 w_{1,2} =
 \frac{ 8 \alpha k - 4 \alpha \delta - 198 \alpha - k
 \pm \sqrt{\Delta} }
 {4 \alpha (2 \delta - 1)},
\]
where \(\Delta =
 16 \alpha^{2}\delta^{2}
 - 64 \alpha^{2}\delta k
 + 64 \alpha^{2}k^{2}
 - 1616 \alpha^{2}\delta
 - 3168 \alpha^{2}k
 + 8 \alpha \delta k
 - 16 \alpha k^{2}
 + 40804 \alpha^{2}
 + 396 \alpha k
 + k^{2}.
 \)
 As \(k \to 0\) we obtain
 \[\lim_{k \to 0} \sqrt{\Delta} = \sqrt{16\,{\alpha}^{2}{\delta}^{2}-1616\,{\alpha}^{2}\delta+40804\,{\alpha}^{2}} = \sqrt{4 \alpha^2(2 \delta -101)^2} = 2 \alpha \vert 2 \delta -101 \vert\]
 and \( \lim_{k \to 0} w_i(k)= w_i^0\), \(i=1,2\) where
 \[
  w_1^0 =
 \frac{ - 2 \delta - 99 - \vert 2 \delta -101 \vert } {2 (2 \delta -1)}, \quad
 w_2^0 = \frac{ - 2 \delta - 99 + \vert 2 \delta -101 \vert } {2 (2 \delta -1)}.
 \]
 A straightforward computation shows that the numerator of $w_2^0$ is positive if
 and only if  \( \delta < \frac{1}{2}\),
 which implies that the denominator is negative. Hence,
 \[
 w_2^0 < 0 \quad \text{for all } \delta \neq \frac{1}{2}.
 \]
 On the other hand, the numerator of $w_1^0$ is always negative.
 Therefore, $w_1^0 > 0$ if and only if
 \(\delta < \frac{1}{2}\).
 We thus obtain the following result.

\begin{proposition}
If $\delta > \frac{1}{2}$, then both $w_i^0$, $i=1,2$, are negative.
 If $\delta \in \left(0,\frac{1}{2}\right)$, then $w_1^0>0$ and $w_2^0<0$.
 \end{proposition}

Henceforth, we assume
 \(\delta \in \left(0,\frac{1}{2}\right)\), and we consider only the value
$w_1^0= \frac{100}{1-2\delta}$. Setting for simplicity $w := w_1^0$
we obtain the equilibrium point
 \[
 P_0=\Big(
 w,\;
 w,\;
 \frac{1}{\delta}-w,\;
 -\frac{2\delta}{\alpha}w
 \Big).
 \]
The Jacobian matrix \eqref{eq:jacobiana_ipersemplificata} evaluated at this
equilibrium point is
\begin{equation}
 J_F(P_0)= \frac{1}{2\delta-1}
 \begin{pmatrix}
 50 & 100\delta & 0 & 100\alpha \\
 -100\delta & 0 & -100\delta & 0 \\
 0 & 1-102\delta & 0 & \frac{(1-102\delta)\alpha}{2\delta} \\
 -400\delta & 400\delta & 0 & 0
 \end{pmatrix}.
\end{equation}
The characteristic polynomial of this matrix,
 \(P_4(\rho)=\rho^4+d_3\rho^3+d_2\rho^2+d_1\rho+d_0\),
 is given by
\begin{align*}
 P_4(\rho)
 &= \rho^4
 - \frac{50}{2\delta-1}\rho^3
 + 100\,\frac{\delta(400\alpha+1-2\delta)}{(2\delta-1)^2}\rho^2 \\
 &\quad  + 5000\,\frac{\delta(392\delta\alpha+4\alpha+102\delta-1)}{(2\delta-1)^3}\rho
 - 10^6\,\frac{\delta\alpha(102\delta-1)}{(2\delta-1)^3}.
 \end{align*}

We notice that $d_0 <0 $ if $\delta < 1/102$, while $d_1 <0$ if
$ \delta \in [1/102 , 1/2 )$. In any case, if $\delta \in (0,1/2)$, applying
Routh-Hurwitz criterion we obtain that $P_4$ has at least a root with strictly
positive real part. This implies the instability of the equilibrium $P$ for small $k$'s
and $\delta \in (0,1/2) $.
We summarize the above analysis in the following proposition.

 \begin{proposition}
 \begin{itemize}
 \item[(i)] If $\delta > \frac{1}{2}$, then there exists $k^*>0$ such that, for all
 $k\in(0,k^*)$, we have $w_{1}, w_2<0$. Hence, the equilibria of case~9) lie outside $
\overline{\mathcal{C}}$.

 \item[(ii)] For all $\delta \in (0,1/2)$, there exists $k(\delta)>0$ such that,
for all $k\in(0,k(\delta))$, the equilibrium associated with $w_1$ is unstable.
 \end{itemize}
 \end{proposition}

\section{Equilibria for the population}

 We now take into account the behavior of the population (the \(x\)-variable).
 Let us consider the equation for $x$ in \eqref{eq:second}, that is
\begin{equation}\label{eq:x6}
 x' = -c x + d \frac{ x w}{w^2 +x^2} \; x - \gamma_0 p x.
\end{equation}
Here $w,p$ are solutions of the system \eqref{eq:secondnonx}, so they are
${\mathcal C}^1$ functions defined in $[0, + \infty )$.
Hence, we can write \eqref{eq:x6} as $x' = g(t,x)$ where $g$ is a ${\mathcal C}^1$
function on $[0, + \infty )^2$. It also satisfies $g (t,0)=0$ and, by an easy
computation, $g(t,x)\leq \frac{d}{2} x.$ It is then easy, using standard results on ODE,
to get the following proposition.

\begin{proposition}
For any $x_0 \geq 0$ there is a unique ${\mathcal C}^1$ solution $x(t, x_0 )= x(t)$
of the Cauchy problem $x' = g(t,x)$, $x(0)=x_0$. The solution is defined in
$[0, +\infty )$. If $x_0=0$ then the solution is the zero constant, $x(t, x_0)=0$.
If $x_0 >0$ then $x(t, x_0 )>0$ for all $t \geq 0$.
\end{proposition}

 We now come back to the study of equilibrium points. Let us consider an equilibrium
$ (x_0 , y_0 ,w_0 , z_0 , p_0 ) $ of the system \eqref{eq:second}. Of course, the
point $P_0 = ( y_0 ,w_0 , z_0 , p_0 )$ is an equilibrium for the system
\eqref{eq:secondnonx}. As the variables $y,w,z,p$ do not depend on $x$, it is easy to
understand that, if $P_0$ is unstable for \eqref{eq:secondnonx}, then
$(x_0 , P_0 )$ is unstable for \eqref{eq:second}. This means that we have to study
the stability of $(x_0 , P_0 )$ only for the cases in which $P_0$ is stable for
\eqref{eq:second}.

Setting $f(x,w,p)$ for the right-hand side of \eqref{eq:x6}, we compute the values
$x_0$ for which $ f(x_0 ,w_0 ,p_0 )=0$, where $w_0 , p_0$ are the components of an
equilibrium $P_0$ of system \eqref{eq:second}. The equations to be solved is
\[
 -c x_0 + d \frac{x_0 w_0}{w_0^2 + x_0^2}\, x - \gamma_0 p_0 x_0 = 0,
\]
which yields \(x_0 = 0\) and
\[
 (c + \gamma_0 p)\, x_0^2 - d w_0\, x_0 + (c + \gamma_0 p_0)\, w_0^2 = 0,
\]
hence
\begin{equation} \label{eq:x12}
 x_{1} = \frac{w_0\big(d - \sqrt{d^2 - 4(c + \gamma_0 p_0 )^2}\big)}
 {2(c + \gamma_0 p_0)}, \quad x_{2}
 = \frac{w_0\big(d + \sqrt{d^2 - 4(c + \gamma_0 p_0 )^2}\big)}
 {2(c + \gamma_0 p_0)}.
 \end{equation}

Notice that we obtain the solutions $x_{1,2}$ only with the hypothesis
 \(d\geq 2(c+ \gamma_0 p_0)\).
If \(w_0 = 0\), then \(x_1 = x_2 = 0\), but the equilibrium \(x = w=0\) is not
admissible for our equation.
 However, in order to study stability, it is not difficult to understand the behavior
of the trajectories of the system near a point $(0, P_0 )$ where $P_0$ is an equilibrium
for \eqref{eq:secondnonx} with \(w_0 = 0\).
Indeed, from our previous study, it is easy to verify such equilibria are only those of
cases (1)--(3), that is, they correspond to half-lines of equilibria.
In fact, if an equilibrium, belonging to one of the other cases, satisfies \(w = 0\),
then it must necessarily lie on one of these half-lines of equilibria.
 Moreover, as we said above, the only situations in which it is interesting to study
the stability of an equilibrium with respect to \(x\) are those cases in which stability
holds for the other variables.
 Hence, recalling that \(w(t) \to 0\) implies \(x(t) \to 0\) as \(t \to +\infty\),
from the results of Section 4 and Proposition \ref{prop:wx} we obtain the following
statement.

\begin{proposition} \label{prop:behax}
For the half-lines of equilibria, whenever stability holds, we also have
\[
 x(t) \to 0 \quad \text{as } t \to +\infty.
\]
\end{proposition}

Proposition \ref{prop:behax} answers, among others, to the question about what
happens to trajectories starting nearby equilibria with $x=w=0$. Notice that the
convergence \(w(t) \to 0\) is guaranteed only for trajectories starting sufficiently
close to the equilibrium, that is, with \(w_0\) close to \(0\), whereas the conclusion
\(x(t) \to 0\) holds for any initial value of \(x\).

 We have now to study equilibria with \(w_0 > 0\). For this, we notice that, if we
compute the Jacobian matrix of \eqref{eq:second} and its eigenvalues at the equilibrium
$(x_0 , P_0 )$, it trivially results that one eigenvalue is given by
 $\frac{\partial f}{\partial x}(x_0 ,w_0 ,p_0 )$ and the others are those of the
Jacobian matrix of \eqref{eq:secondnonx}. Now we compute
\[
 \frac{\partial f}{\partial x}(x,w,p)
 = -c - \gamma_0 p + 2 d \frac{w^3 x}{(w^2 + x^2)^2},
\]
and therefore
\[
 \frac{\partial f}{\partial x}(0,w,p)
 = -c - \gamma_0 p < 0,
 \quad \forall w,p\geq 0.
\]

On the other hand, at the points \(x_1, x_2\) defined above, we have
 \[
 \frac{d x w}{w^2 + x^2} = c + \gamma_0 p,
 \]
 where, for the sake of simplicity, we drop the indices and write $x$ for $x_i$,
$w $ for $w_0$, $p$ for $p_0$. We then obtain
 \[
 \frac{\partial f}{\partial x}(x,w,p)
 = -c - \gamma_0 p
 + 2 \frac{d w x}{w^2 + x^2} \frac{w^2}{w^2 + x^2}
 = -c - \gamma_0 p
 + 2 (c + \gamma_0 p)\frac{w^2}{w^2 + x^2} = (c + \gamma_0 p)\frac{w^2 - x^2}{w^2 + x^2}.
 \]
Let \( x_1 \) be defined as in \eqref{eq:x12}.
 Then \(x_1 < w\) if and only if
 \[
 d - \sqrt{d^2 - 4(c + \gamma_0 p)^2}
 < 2(c + \gamma_0 p),
\]
 which is equivalent to
\[
 d - 2(c + \gamma_0 p)
 < \sqrt{d^2 - 4(c + \gamma_0 p)^2}.
\]
Since \(d > 2(c + \gamma_0 p)\) by hypothesis, this inequality is equivalent to
\begin{align*}
&d^2 + 4(c + \gamma_0 p)^2 - 4d(c + \gamma_0 p)
 < d^2 - 4(c + \gamma_0 p)^2 \\
&\iff   8(c + \gamma_0 p)^2 < 4d(c + \gamma_0 p) \\
 &\iff  2(c + \gamma_0 p)  < d.
 \end{align*}
Hence, in our hypotheses, we always have $x_1 <w$.

On the other hand, for the second root \(x_2\), we have \(x_2 > w\) if and only if
\[
 \frac{w\left(d + \sqrt{d^2 - 4(c + \gamma_0 p)^2}\right)}
 {2(c + \gamma_0 p)}
 > w
 \iff
 \sqrt{d^2 - 4(c + \gamma_0 p)^2}
 > 2(c + \gamma_0 p) - d,
\]
which is true since \(d > 2(c + \gamma_0 p)\).
 Consequently,
 \[
 \frac{\partial f}{\partial x}(x_1,w,p) > 0,
 \quad
 \frac{\partial f}{\partial x}(x_2,w,p) < 0.
 \]
We can therefore state the following proposition. Remember that it deals with
isolated equilibria, in which we have obtained stability from the study of the
eigenvalues of the linearized system.

 \begin{proposition}
 Let \(P_0 = (y_0, w_0, z_0, p_0)\) be any asymptotically stable equilibrium of
\eqref{eq:secondnonx} with \(w_0 > 0\).
 Then:
 \begin{itemize}
 \item[(i)] the equilibria \((0, P_0)\) and \((x_2, P_0)\) of \eqref{eq:second}
 are asymptotically stable;

 \item[(ii)] the equilibrium \((x_1, P_0)\) of \eqref{eq:second} is unstable.
 \end{itemize}
 \end{proposition}

% page 22


 \section{Numerical results}

In this section we first report the simulations corresponding to case~1, in which we 
have a half-line of equilibria with nonnegative pollution. Both the stable and the 
unstable regimes are illustrated. We then present the results concerning stable 
configurations with positive population, arising in the specific parameter settings 
described in cases~7 and~8.

 The numerical simulations were performed using the stiff solver \texttt{ode15s},
 since for the considered parameter values explicit methods (such as \texttt{ode45})
 require an excessively large number of time steps due to the presence of widely
 separated time scales during the transient dynamics and the non-negativity
 constraints imposed on the state variables.

 In accordance with \cite{Ton23}, we set $\gamma = 0.5$ and $\lambda = 100$, and we 
simulated the scenario arising from a perturbation of the critical points 
 \( P_0 = (y_0, w_0, z_0, p_0) \) of the simplified system \eqref{eq:first_simply}, 
together with the behavior of the population \( x \) resulting from a perturbation of 
\( x_0 \), with \( x_0 = 0 \) or \( x_0 = x_2 \) given by \eqref{eq:x12}. 
Furthermore, when \( x_2 \) is considered, namely in cases~7 and ~8, in order to 
ensure that the equilibrium point satisfies the admissibility condition
\( x_2 \in \RR_+ \), we set \( c = 1 \) and choose
\( d = 2\,(c + \gamma_0 p_0) + 1 \).

\subsection{Simulation case~1}
In this case, the numerical simulations are performed with the parameter
set $c=1, \, d=3$, $\gamma_0=1$, $\delta= 0.1$,
$k= 1$, $\alpha=1$.
The pollution level ${p_0}$ is varied in order to illustrate both regimes:
${p_0}=40$ for the unstable case and ${p_0}=60$ for the stable case, while the critical
threshold is given by $p_{\mathrm{crit}}=\gamma\lambda/\alpha=50$. 
We will choose for $x$ small initial values.

In the stable regime, all state variables remain nearby the equilibrium values and 
all but $p(t)$ converge towards them. In particular, renewable and non-renewable 
resources do not exhibit any transient overshoot. 
According to the theoretical analysis, when the initial value $p_0$ is perturbed, 
the pollution level $p(t)$ does not necessarily converge to $p_0$. 
This behavior is evident by Figure~\ref{fig:caso2_stable}. The left panel shows 
the asymptotic stability of $x,y,w,z$, while the right panel provides a numerical 
illustration of Lyapunov stability of the equilibrium $P=(0, 0,0,0, p_0)$ with respect 
to the $p$ variable.
The shaded region represents the neighbourhood
$[ p_0-\eta, p_0+\eta]$ with $\eta=0.15$.
For sufficiently small perturbations of size $\epsilon=10^{-3}$ and
$\epsilon=10^{-4}$ (applied to all components of $P$), the trajectory
$p(t)$ remains confined within the prescribed band for all $t\ge0$,
whereas a larger perturbation ($\epsilon=10^{-2}$) may exit it.
This illustrates the dependence of the admissible neighbourhood on the perturbation size.

\begin{figure}[htbp]
 \centering
  \includegraphics[width=0.45\textwidth]{caso1_stabilexywz}
 \quad 
\includegraphics[width=0.45\textwidth]{lyap2_v2}
\caption{Case~1. Stable regime for the equilibrium family $(0,0,0,0,60)$. Left: state variables $x, \, y, \, w, \, z$. Right: trajectories corresponding to perturbations
$\epsilon=10^{-3},10^{-4}$ and $\epsilon=10^{-2}$.}
\label{fig:caso2_stable}
\end{figure}

Figure~\ref{fig:caso2_unstable} shows the corresponding unstable regime.
In this case, the dynamics exhibit a short-lived boom--bust pattern.
Renewable resources $y$ rapidly increase towards their carrying capacity,
which in turn triggers a sharp transient growth of wealth $w$.
However, as the system evolves, the effective growth factor becomes negative
and $w$ subsequently declines towards zero, a trend already clearly visible by $t=8$. 
Also the $x , y, z$ variables then decline to zero. 
On the other hand, figure~\ref{fig:caso2_unstable} (right panel) shows the time 
evolution of the pollution variable $p(t)$ on a logarithmic scale. 
After a short initial transient, $p(t)$ exceeds the critical value $p_0$ and 
grows rapidly over several orders of magnitude,
revealing a runaway behaviour of pollution and confirming the instability of the 
equilibrium configuration.
For visualization purposes, extremely small values of $p$ are truncated
to avoid numerical artefacts associated with the logarithmic scale.

A direct comparison between Figures~\ref{fig:caso2_stable} 
and~\ref{fig:caso2_unstable} highlights the role of the pollution level ${p_0}$: 
keeping all other parameters fixed, crossing the  threshold 
$\alpha p_0 = \gamma \lambda$ qualitatively changes the dynamics of the system.

From a dynamical systems perspective, pollution enters the model as a damping mechanism.
When the parameter condition ensuring stability is satisfied, this damping effect 
dominates the growth terms, yielding bounded trajectories and convergence towards 
the equilibrium manifold.
If the condition is violated, the growth terms prevail, leading to transient expansion 
which, through the influence of $w$ on the equation for $p$, implies a boost of 
pollution, hence a decline of all the other state variables.

It is interesting to note that all these results may be seen as social collapses. 
They mean that, when a society passes close to the values of the half-line we are 
dealing with, then population and wealth will collapse to zero, both in the stable 
and in the unstable case. In the latter case, the collapse is preceded by a growth 
phase, and this dynamic may remind one of the dynamics observed in the well-known models of \cite{meadows1972limits}.


 \begin{figure}[htbp]
 \centering
 \includegraphics[width=0.45\textwidth]{caso1_instabilexywz}
 \quad
 \includegraphics[width=0.45\textwidth]{caso1_instabilelogp}
 \caption{Case 1. Unstable regime for the equilibrium family $(0,0,0,0,40)$.
 Left: state variables $x, \, y, \, w, \, z$.
 Right: pollution $p$ shown on a logarithmic scale. }
 \label{fig:caso2_unstable}
 \end{figure}
 
 \subsection{Simulation case 7}
In this case, consistently with the theoretical analysis, we investigate the 
model in the limit  $k, \alpha \to 0$. For values $k > 10^{-2}$, the equilibrium point 
$P_0$ does not belong to the admissible
set $\overline{\mathcal{C}}$, while for $\alpha > 10^{-4}$ the Jacobian matrix 
exhibits positive eigenvalues.
Accordingly, we fix $k = 10^{-2}$ and $\alpha = 10^{-4}$. The point under 
consideration is \((x_2, P_0)\), with \(x_2\) as in \eqref{eq:x12}.


 Moreover,  to prevent the large pollution levels from driving the population 
toward the attractive critical state $x = 0$, we reduce the pollution-related 
mortality coefficient to  $\gamma_0 = 10^{-2}$. This choice weakens the impact of 
pollution on population dynamics while preserving the analytical structure of 
the model.
 With these parameter values, we obtain the equilibrium point
 \[
 (x_0, y_0, w_0, z_0, p_0) = (60.1229,\, 10,\, 59.1690,\, 0,\, 3.9 \times 10^{5}),
 \]
whose Jacobian matrix evaluated at this point has eigenvalues
 \[
 \rho_1 = -62.5185,\quad
 \rho_2 = -3.2429,\quad
 \rho_3 = -1.6730,\quad
 \rho_4 = -0.0841,\quad
 \rho_5 = -25.4486.
 \]

Because of the different orders of magnitude of the equilibrium coordinates, 
we present the time  evolution of the population $x$ and wealth $w$ in the left panel 
of Figure~\ref{fig:caso14_stable},  while the behavior of the pollution variable $p$ 
is shown in the right panel. The remaining  variables $y$ and $z$ exhibit a very fast 
convergence toward their equilibrium values and are  therefore omitted for clarity.

\begin{figure}[htbp]
 \centering
\includegraphics[width=0.45\textwidth]{caso14xw}
\quad
\includegraphics[width=0.45\textwidth]{caso14p}
 \caption{Case~7. Stable regime for the equilibrium family $(x_2,P_0)$ with \(x_2\) as in \eqref{eq:x12} and \(P_0\) as in case 7.
 Left: state variables $x, \, w$.
 Right: pollution $p$ shown on a logarithmic scale.}
 \label{fig:caso14_stable}
 \end{figure}

\subsection{Simulation case~8}
In this scenario, the following constraints are required in order for 
\(P_0 \in \mathcal{\overline{C}}\):
\begin{gather*}
  \delta \le k \le 50 + \delta \quad  \text{if } \delta \le \frac 1 {100} \\
  \delta +50 - \frac 1 {2 \delta} \le k \le \delta +50 \quad  \text{if } \delta > \frac 1 {100}
 \end{gather*}
In Figure~\ref{fig:caso15_stable}, we report the results obtained for $\delta=10$, 
$k=59.98$, and $\alpha=10^{-4}$. The equilibrium point $(x_0,P_0)$ has coordinates
 \[
 x_0 = 13.0849, \quad
 y_0 = 0.0400, \quad
 w_0 = 4.9980, \quad
 z_0 = 0.0600, \quad
 p_0 = 0.
 \]
 The eigenvalues of the Jacobian matrix are
 \[
 \rho_1 = -0.7454,\quad
 \rho_{2,3} = -0.0044 \pm 6.7067 i,\quad
 \rho_4 = -0.0111,\quad
 \rho_5 = -49.9790.
 \]
 Hence the equilibrium is asymptotically stable. The presence of the complex conjugate 
pair $\rho_{2,3}$ with a small negative real part explains the observed slowly damped 
oscillations, suggesting that the system is operating close to a Hopf threshold.

A qualitatively similar behaviour is obtained for $\delta=10$, $k=59.98$, and 
$\alpha=10^{-1}$: the equilibrium and the transient dynamics remain essentially 
unchanged, while only the fast mode (associated with $\rho_5$) is modified. 
We do not report the graph for this case.

\begin{figure}[htbp]
 \centering
 \includegraphics[width=0.45\textwidth]{caso15k5998xw}
 \quad
 \includegraphics[width=0.45\textwidth]{caso15k5998yz}
 \caption{Case~8: \(\delta=10\), \(k = 59.98\) and \(\alpha = 10^{-4}\).}
 \label{fig:caso15_stable}
 \end{figure}

We observe that, as $\delta$ increases, the admissible interval for $k$ becomes 
progressively narrower and the real part of the complex eigenvalues approaches zero, 
resulting in increasingly weakly damped oscillations. Moreover, varying $\alpha$ from 
$10^{-4}$ to $10^{-1}$ does not significantly affect the qualitative dynamics in the 
considered regimes. We do not report the graphs for these cases.

 In Figure~\ref{fig:esempio15_1}, we show the scenario with 
\(\delta = 0.01\), \(k = 40\) and \(\alpha=10^{-4}\). The point \((x_0,P_0)\) has 
coordinates
\[
 x_0 = 1.0470 \cdot 10^{4}, \quad
 y_0 = 20.0200, \quad
 w_0 = 3999, \quad
 z_0 = 79.9800, \quad
 p_0 = 0.
\]
The eigenvalues of the Jacobian matrix are
 \[ \rho_1 = -0.7454, \quad
 \rho_{2,3} = -0.3129 \pm 5.8319i, \quad
 \rho_4 = -9.3841, \quad
 \rho_5 = -39.1942.
 \]
Also in this case the equilibrium is asymptotically stable, and the complex pair 
$\rho_{2,3}$ accounts for the damped oscillations observed in the trajectories of 
$y$, $w$, and $z$.

 \begin{figure}[!htbp]
 \centering
 \includegraphics[width=0.9\textwidth]{caso15k40_v2}
 \caption{Case~8: \(\delta = 0.01\), \(k = 40\) and \(\alpha=10^{-4}\).}
 \label{fig:esempio15_1}
 \end{figure}

Stability for large values of $\alpha$ appears to be rather unlikely, as it occurs 
only under restrictive parameter conditions. In particular, it can be shown that 
the condition stated in item~(iii) of Proposition~\ref{prop:19} is never satisfied for 
$\delta > \frac{1}{98}$.

 If $\delta < \frac{1}{100}$, the condition holds for
 \(k \in [\delta, \delta + 50)\).
 If $\frac{1}{100} < \delta < \frac{1}{98}$, the condition holds for
 \(k \in \left(\delta + 50 - \frac{1}{2\delta}, \, \delta + 50\right)\).
 For example, in the simulation obtained with \( \delta = 0.0025 \), 
 \( k = \delta + 10^{-3} \), and \( \alpha = 10^{2} \), we obtained the equilibrium point
 \(P = (1.0472, 99.9980, 0.4000, 300.0020, 0)\),
 with eigenvalues of the Jacobian matrix computed at this point given by
 \[ \rho_1 = -0.7454, \quad
 \rho_{2,3} = -2.5 \cdot 10^{-6} \pm 0.0146i, \quad
 \rho_4 = -49.9990, \quad
 \rho_5 = -1.9920 \cdot 10^{4}.
 \]
Convergence occurs over very long time intervals.
 We do not report the corresponding plots, as the system is strongly stiff.
 In particular, the eigenvalues of the Jacobian matrix exhibit a pronounced 
separation of time scales: one component decays very rapidly (of order $10^{4}$), 
while another evolves on a much longer time scale, with characteristic time of order 
$4 \times 10^{5}$.
As a consequence, after a very fast transient, the solution approaches the equilibrium 
extremely slowly, making the graphical representation scarcely informative. 
As \( \alpha \) increases, the equilibrium point remains essentially unchanged;
however, the magnitude of the largest eigenvalue grows accordingly.
For example, for \( \alpha = 10^{4} \), the largest eigenvalue is of the order of 
\( 10^{6} \), thus making the Jacobian matrix increasingly stiff.

\section{Conclusions}
 
Let us summarize the main results of this paper.
\begin{itemize}
\item In our model a steady growth in wealth is not possible. As we have said in 
 Remark \ref{re:osc}, the variable $w$ either is bounded or, if it is not, it 
 experiences strong oscillations. The quest for a perpetual economic growth seems
 either hopeless or even dangerous. These results agree with what we have obtained in 
 our previous papers.
 
\item The collapse of the society is a real possibility. The results about stability of 
 equilibria with $w_0 =0$ mean that, if a trajectory approaches such a stable equilibrium, then there will be a collapse of wealth and population. A very interesting result in section 7 is the fact that, for these equilibria, collapse seems to occur {\it also in the case of instability}.
 Indeed, Figure 2 in Section 7 shows a behavior reminiscent of the World3 model 
 \cite{meadows1972limits}: an initial phase of strong growth followed by collapse.
 In our model this evolution seems due to the role of pollution $p$, which grows very 
 rapidly and causes the collapse of wealth and resources. We have no theoretical result in this case, so it seems that there is still much work to do here
  about the possible collapse paths.

 \item We are able to identify stable equilibria with positive values of $x,w$ 
 (see in particular cases~7 and ~8 in section 5 and the corresponding numerical 
 simulation in section 7). The existence of these equilibria is good news because 
 they represent safe havens towards which to direct the evolution of human societies 
 and thus avoid collapse. We obtained similar results in our previous papers, 
 but in the present case it seems that equilibria of this kind are obtained only for 
 very particular values of the parameters. This worsening of the situation seems due to 
 the introduction of pollution in the model.
 
\item Roughly speaking, a general indication is that the role of pollution is to make 
 it more difficult for society to evolve towards safe havens, and to make the evolution 
 towards social collapse more likely.
\end{itemize}

\begin{thebibliography}{10}

\bibitem{AkYor20}
N.~Akhavan, A.~Yorke;
\emph{Population collapse in elite-dominated
  societies: a differential model without differential equations}, SIAM J.
  Appl. Dyn. Syst. \textbf{19} (2020), 1736--1757.

\bibitem{BaCra24a}
M.~Badiale, I.~Cravero;
\emph{A handy-type model with non-renewable
  resources}, Nonlinear Anal. Real World Appl. \textbf{77} (2024), 104071.

\bibitem{BaCra24b}
M.~Badiale, I.~Cravero;
\emph{A nonlinear {ODE} model for a consumeristic society},
  Mathematics \textbf{12} (2024), no.~8, 1253.

\bibitem{BaCra25b}
M.~Badiale, I.~Cravero;
\emph{Interactions between wealth and natural resources: a nonlinear
  {ODE} model}, Axioms \textbf{14} (2025), no.~5, 430.

\bibitem{BaCra25a}
M.~Badiale, I.~Cravero;
\emph{A nonlinear ordinary differential equation model for a society
  strongly dependent on non-renewable resources}, AppliedMath \textbf{5}
  (2025), no.~1, 43.

\bibitem{Bardi2024PopulationCollapse}
U.~Bardi;
\emph{The coming population collapse}, The Seneca Effect (Substack
  blog post), Dec 2024, Published 9 December 2024,
 https: //senecaeffect.substack.com/p/the-coming-population-collapse.

\bibitem{Calvin2018}
K.~Calvin, B.~Bond-Lamberty;
\emph{Integrated human-earth system
  modeling---state of the science and future directions}, Environmental
  Research Letters \textbf{13} (2018), 063006,
  DOI: 10.1088/1748-9326/aac642.

\bibitem{GraWilSat20}
B.~Grammaticos, R.~Willox, J.~Satsuma;
 \emph{Revisiting the human and
  nature dynamics model}, Regul. Chaotic Dyn. \textbf{25} (2020), no.~2,
  178--198.

\bibitem{herrington2020update}
G.~Herrington;
\emph{Update to limits to growth: Comparing the world3 model
  with empirical data}, Journal of Industrial Ecology (2020), 1--13.

\bibitem{meadows2004limits}
D.~H. Meadows, D.~L. Meadows,  J.~Randers;
\emph{The limits to growth: The
  30-year update}, Chelsea Green Publishing Co., White River Junction, VT,
  2004.

\bibitem{meadows1972limits}
D.~H. Meadows, D.~L. Meadows, J.~Randers,  W.~W. Behrens;
 \emph{The limits   to growth: A report for the club of rome’s project on the predicament of
  mankind}, Universe Books, New York, 1972.

\bibitem{MotRivKal14}
S.~Motesharrei, J.~Rivas, E.~Kalnay;
\emph{Human and nature dynamics
  ({HANDY}): modeling inequality and use of resources in the collapse or
  sustainability of societies}, Ecol. Econ. \textbf{101} (2014), 90--102.

\bibitem{Nebel2024}
A.~Nebel, A.~Kling, R.~Willamowski, T.~Schell;
\emph{Recalibration of   limits to growth: An update of the World3 model}, 
  Journal of Industrial
  Ecology \textbf{28} (2024), 87--99,  DOI: 10.1111/jiec.13442.

\bibitem{Nebel2025Corrigendum}
A.~Nebel, A.~Kling, R.~Willamowski, T.~Schell;
\emph{Correction to ``recalibration of limits to growth: An update of
  the World3 model''}, Journal of Industrial Ecology \textbf{29} (2025),
  1443--1444, DOI: 10.1111/jiec.70042

\bibitem{Nitzbon}
J.~Nitzbon, J.~Heitzig, U.~Parlitz;
\emph{Sustainability, collapse and
  oscillations in a simple World-Earth model}, Environmental Research Letters
  \textbf{12} (2017), 074020.

\bibitem{Perissi2019}
I.~Perissi;
\emph{Highlighting the archetypes of sustainability management by
  means of simple dynamics models}, Journal of Simulation (2019),
  DOI: 10.1080/17477778.2019.1679612.

\bibitem{ShiKad23}
M.~Shillor, T.~A. Kadhim;
\emph{Analysis and simulation of the handy model
  with social mobility, renewables, and nonrenewables}, Electron. J. Differ.
  Equ. \textbf{2023} (2023), 1--22.

\bibitem{Tan2023}
J.~Tan, Q.~Duan, C.~Xiao, C.~He, X.~Yan;
\emph{A brief review of the
  coupled human-{Earth} system modeling: Current state and challenges}, The
  Anthropocene Review (2023), 1--21.

\bibitem{Thissen1978}
W.~A.~H. Thissen;
\emph{Investigations into the club of Rome's World3 model:
  Lessons for understanding complicated models}, PhD thesis, Technische
  Hogeschool Eindhoven, Eindhoven, The Netherlands, 1978.

\bibitem{Ton23}
A.~Tonnelier;
\emph{Sustainability or societal collapse: dynamics and
  bifurcations of the handy model}, SIAM J. Appl. Dyn. Syst. \textbf{22}
  (2023), no.~3, 1877--1905.

\bibitem{Vermeulen1976}
P.~J. Vermeulen, D.~C.~J. de~Jongh;
\emph{Parameter sensitivity of the
  `limits to growth' world model}, Applied Mathematical Modelling \textbf{1}
  (1976), 29--32.

\end{thebibliography}

\end{document} 