General remarks
Usually, when I tackle this kind of problem, I'm looking at what part of the space is forward invariant first. It allows me to better understand what I have to consider in my analysis. Here, the SIS-like models are usually confined in $\mathbb{R}^N_+$ since the quantities are the number of individuals in each category.
An attracting region of interest:
Let $$\mathcal{S} = \left\{ (S, \mathbf{x}) \in \mathbb{R}_+^5 \mid S(t) \le S^* \right\} = \left[ 0 , S^* \right] \times \mathbb{R}_+^4.$$
From the first equation of your system, $$\dot{S} = \mu - (\text{infection terms}) - \nu S.$$
Since the infection terms are non-negative, we have the inequality $$\dot{S} \le \mu - \nu S.$$
Using Gronwall implies $S(t) \leq S^* = \mu/\nu$ for large $t$. i.e., the region $\mathcal{S}$ is attracting.
Moreover, this region is forward invariant. To prove that, we can use Nagumo's theorem, which states that it is sufficient to prove that the vector field is pointing inward into your set on its border. We can look at the $1$D-dynamics of $S$ over the interval $[0, S^*]$. At the edges, one has
\begin{align*} \dot{S}_{|S = S^*} &= \mu - \beta_1 S^* I_2 - \beta_2 S^* J - \beta_3 S^* A - \nu S^*\\ &\leq \mu - \nu S^* = \mu - \nu \frac{\mu}{\nu} = 0. \end{align*} and \begin{align*} \dot{S}_{|S = 0} &= \mu - \beta_1 S^* I_2 - \beta_2 S^* J - \beta_3 S^* A - \nu S^*\\ &= \mu \geq 0. \end{align*}
Then, $\mathcal{S}$ is attracting and is forward invariant.
This fact allows us to simplify our life and only study the dynamics on $\mathcal{S}$. We will, in the following, even remove $S$ explicitly from our study since we will restrict ourselves to $\mathcal{S}$.
Cleaner formulation (vectorial notation)
The primary objective is to distinguish the dynamics of the infectious compartments from those of the susceptible compartment. This allows us to represent linear transitions and non-linear new infections using matrices, which are easier to analyze. Moreover, the vectorial matrix is much more natural to think about in terms of Lyapunov functions, in my opinion.
Let the vector of infectious compartments be $\mathbf{x} = (I_1, I_2, J, A)^\top$.
The equations for these infectious compartments can be written in the general form, for any $(S,\mathbf{x}) \in \mathcal{S}$,
$$ \dot{\mathbf{x}} = \mathbf{f}(S, \mathbf{x}) - \mathbf{V}\mathbf{x} \tag{1}$$
where,
$\mathbf{f}(S, \mathbf{x})$: the nonlinear coupling between $S$ and the other quantities. For a given number of susceptible $S$, one can rewrite the function as, $$ \mathbf{f}(S, \mathbf{x}) = %\begin{bmatrix} %p \beta_1 S I_2 + q \beta_2 S J + r \beta_3 S A \\ %(1-p) \beta_1 S I_2 + (1 - q) \beta_2 S J + (1 - r) \beta_3 S A \\ %0 \\ %0 %\end{bmatrix} = S \begin{bmatrix} p \beta_1 I_2 + q \beta_2 J + r \beta_3 A \\ (1-p) \beta_1 I_2 + (1 - q) \beta_2 J + (1 - r) \beta_3 A \\ 0 \\ 0 \end{bmatrix} = S \begin{bmatrix} 0 & p\beta_1 & q\beta_2 & r\beta_3 \\ 0 & (1-p)\beta_1 & (1-q)\beta_2 & (1-r)\beta_3 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix}\begin{bmatrix} I_1 \\ I_2 \\ J \\ A \end{bmatrix} := S\, \mathbf{F} \mathbf{x}.$$
$\mathbf{V}\mathbf{x}$: all other linear dynamics, $$ \begin{bmatrix} b_1 I_1 - \xi_1 J \\ b_2 I_2 - \epsilon I_1 - \xi_2 J \\ b_3 J - p_1 I_2 \\ b_4 A - p_2 J \end{bmatrix} = \begin{bmatrix} b_1 & 0 & -\xi_1 & 0 \\ -\epsilon & b_2 & -\xi_2 & 0 \\ 0 & -p_1 & b_3 & 0 \\ 0 & 0 & -p_2 & b_4 \end{bmatrix} \begin{bmatrix} I_1 \\ I_2 \\ J \\ A \end{bmatrix} := \mathbf{V} \mathbf{x} $$
This matrix $\mathbf{V}$ has the following properties:
- its off-diagonal entries are non-positive (it's a $Z$-matrix).
- its inverse, $\mathbf{V}^{-1}$, exists and all its entries are non-negative ($\mathbf{V}^{-1} \ge 0$).
Remark: I didn't prove it theoretically, but I used symbolic calculus to confirm that. You will need to prove that theoretically, but it's fastidious (due to the notation), so I skipped that. I used Python's library Sympy. P.S.: I can think about an efficient proof.
import sympy mu, nu = sympy.symbols('mu nu', positive=True) epsilon, p1, p2 = sympy.symbols('varepsilon p_1 p_2', positive=True) alpha, xi1, xi2 = sympy.symbols('alpha xi_1 xi_2', positive=True) b1 = epsilon + mu b2 = p1 + nu b3 = xi1 + xi2 + p2 + nu b4 = alpha + nu V = sympy.Matrix([ [b1, 0, -xi1, 0], [-epsilon, b2, -xi2, 0], [0, -p1, b3, 0], [0, 0, -p2, b4] ]) adj_V = V.adjugate() all_elements_non_negative = True for i in range(adj_V.rows): for j in range(adj_V.cols): element = sympy.expand(adj_V[i, j]) is_non_neg = element.is_nonnegative print(f"Element ({i},{j}) is non-negative? {is_non_neg}") if not is_non_neg: all_elements_non_negative = False print(f" --> Element expression: {element}") print(f"Is the entire adjugate matrix non-negative? {all_elements_non_negative}")
Since $\mathbf{V}^{-1}$ is a non-negative matrix, we know some things about it. By Perron-Frobenius, it has a leading eigenvalue for which the associated eigenvectors are non-negative. More on that later.
Local stability of Disease Free Equilibrium (DFE)
Note that at an equilibria $\mathbf{x}^*$, $\mathbf{f}(S^*, \mathbf{x}^*) - \mathbf{V} \mathbf{x}^* = 0$, then,
$$ \mathbf{f}(S^*, \mathbf{x}^*) - \mathbf{V} \mathbf{x}^* = 0 \Leftrightarrow S^* \mathbf{F} \mathbf{x}^* = \mathbf{V} \mathbf{x}^* \Leftrightarrow \mathbf{x}^* = S^* \mathbf{V}^{-1} \mathbf{F} \mathbf{x}^* \Leftrightarrow 0 = (S^* \mathbf{V}^{-1} \mathbf{F} - \mathbf{I}) \mathbf{x}^* $$
Let's look at the linearization of the dynamics around $\mathbf{x}^*$. For that, we look at the Jacobian matrix of the dynamics (1).
Still with Sympy, I have the expression of $\mathbf{V}^{-1}$
$$ \mathbf{V}^{-1} = \frac{1}{b_1b_2b_3b_4 - b_4 p_1(b_1\xi_2 + \epsilon\xi_1)} \begin{bmatrix} -b_4(p_1\xi_2 - b_2b_3) & b_4 p_1 \xi_1 & b_4 b_2 \xi_1 & 0 \\ b_4 b_3 \epsilon & b_4 b_1 b_3 & b_4(b_1\xi_2 + \epsilon\xi_1) & 0 \\ b_4 p_1 \epsilon & b_4 b_1 p_1 & b_4 b_1 b_2 & 0 \\ p_1 p_2 \epsilon & b_1 p_1 p_2 & b_1 b_2 p_2 & -p_1(\epsilon\xi_1 + b_1\xi_2) - b_1b_2b_3 \end{bmatrix} $$
Then, the dominant non-zero eigenvalue of the matrix $\mathbf{J} = \mathbf{V}^{-1} \mathbf{F} $ is given by: $$ \lambda_{\max} = \frac{S^* \left( \beta_1 b_3 (b_1(1-p) + p\epsilon) + \beta_2 p_1 (b_1(1-q) + q\epsilon) + \beta_3 p_1 p_2 (b_1(1-r) + r\epsilon)/b_4 \right)}{b_1b_2b_3 - p_1(b_1\xi_2 + \epsilon\xi_1)} = \mathcal{R}_0$$
Then, the DFE is locally stable if the maximal eigenvalue of $\mathbf{V}^{-1}\mathbf{F} <1 \Leftrightarrow \mathcal{R}_0 <1$ by the Linearization principle.
Global Asymptotic Stability of DFE via Lyapunov argument
Your attempt with $V = I_1+I_2+J+A$ is a good start, but fails because it weights each infectious compartment equally. The key insight of the correct method is to weight each compartment by its relative contribution to generating future infections.
The Lyapunov Function
We use the linear Lyapunov function: $$ L(\mathbf{x}) = \mathbf{w}^\top \mathbf{V}^{-1} \mathbf{x} $$ where $\mathbf{w}^\top$ is the left eigenvector of the matrix $\mathbf{V}^{-1}\mathbf{F}$ corresponding to its largest eigenvalue. i.e.,$$ \mathbf{w}^\top (\mathbf{V}^{-1}\mathbf{F}) = \mathcal{R}_0 \mathbf{w}^\top $$
Indeed,
- since $\mathbf{F} \ge 0$ and $\mathbf{V}^{-1} \ge 0$ (at least numerically), their product $\mathbf{V}^{-1}\mathbf{F}$ is a non-negative matrix. The Perron-Frobenius theorem for non-negative matrices guarantees the existence of a real, non-negative eigenvalue equal to the spectral radius ($\mathcal{R}_0$) and a corresponding non-negative left eigenvector $\mathbf{w}^\top$.
- since $\mathbf{w} \ge 0$, $\mathbf{V}^{-1} \ge 0$, and $\mathbf{x} \ge 0$ in the feasible region, it follows that $L \ge 0$. Furthermore, because $\mathbf{w}$ is non-zero and all entries of $\mathbf{V}^{-1}$ are non-negative (I think we should be able to prove there are positive not only non-negative), $L=0$ if and only if $\mathbf{x}=\mathbf{0}$. This confirms $L$ is a valid Lyapunov function candidate.
Remark: In your statement, you didn't specify, but I think you must show that the cone $\mathbb{R}^5_+$ is forward invariant. This is the feasible region I was referring to above. It makes sense, but you didn't specify it, so I let you prove it.
Theorem: Under the assumption $\mathcal{R}_0 \leq 1$. The DFE is globally asymptotically stable for the dynamics given by (1).
Proof: The proof is divided into 3 steps.
Differentiate the Lyapunov function: $$ \dot{L}(\mathbf{x}) = \mathbf{w}^\top \mathbf{V}^{-1} \dot{\mathbf{x}} $$
Now, substitute the system dynamics by writing the new infection term as $\mathbf{f}(S, \mathbf{x}) = \frac{S(t)}{S^*} \mathbf{F}\mathbf{x}$. The full system for the compartments is $\dot{\mathbf{x}} = \frac{S(t)}{S^*} \mathbf{F}\mathbf{x} - \mathbf{V}\mathbf{x}$. Then, \begin{align*} \dot{L}(\mathbf{x}) &= \mathbf{w}^\top \mathbf{V}^{-1} \left( \frac{S(t)}{S^*} \mathbf{F}\mathbf{x} - \mathbf{V}\mathbf{x} \right) = \frac{S(t)}{S^*} \mathbf{w}^\top \mathbf{V}^{-1} \mathbf{F}\mathbf{x} - \mathbf{w}^\top (\mathbf{V}^{-1} \mathbf{V})\mathbf{x} \\ &= \frac{S(t)}{S^*} \mathbf{w}^\top (\mathbf{V}^{-1} \mathbf{F})\mathbf{x} - \mathbf{w}^\top \mathbf{x} \end{align*}
Now we substitute the property of our weight vector, $\mathbf{w}^\top (\mathbf{V}^{-1}\mathbf{F}) = \mathcal{R}_0 \mathbf{w}^\top$, into the equation. $$ \dot{L}(\mathbf{x}) = \frac{S(t)}{S^*} (\mathcal{R}_0 \mathbf{w}^\top) \mathbf{x} - \mathbf{w}^\top \mathbf{x} = \left( \frac{S(t)}{S^*} \mathcal{R}_0 - 1 \right) \mathbf{w}^\top \mathbf{x} $$
An attracting region where the derivative is always negative:
We know from the previous discussion that $$\mathcal{S} = \left\{ (S, \mathbf{x}) \in \mathbb{R}_+^5 \mid S(t) \le S^* \right\},$$
is 'globally' attracting and forward invariant. Then, there exists some time $t\geq 0$ such that $S(t) \in \mathcal{S}$. (We have to be cautious here, need to think more.) We then only need to look at the behaviors of our Lyapunov function on $\mathcal{S}$. For any $S \in \mathcal{S}$, we have $\frac{S}{S^*} \le 1$.
If we assume $\mathcal{R}_0 \le 1$, then the product $\frac{S}{S^*} \mathcal{R}_0 \leq \mathcal{R}_0 \leq 1$. Therefore, the term $\left( \frac{S(t)}{S^*} \mathcal{R}_0 - 1 \right)$ must be less than or equal to zero. Since we already established that $\mathbf{w}^\top \mathbf{x} \ge 0$, we have proven that, for any $(S, \mathbf{x}) \in \mathcal{S}$, $$ \dot{L}(\mathbf{x}) \leq 0.$$
Apply LaSalle's Invariance Principle for global stability: To prove asymptotic stability, we must identify the largest invariant set in $\mathcal{S}$ where $\dot{L}(\mathbf{x}) = 0$.
$\dot{L}(\mathbf{x}) = 0$ requires that $\left( \frac{S(t)}{S^*} \mathcal{R}_0 - 1 \right) \mathbf{w}^\top \mathbf{x} = 0$. This means that on the invariant set, for all time, either $\mathbf{x}(t)=\mathbf{0}$ or $\frac{S(t)}{S^*} \mathcal{R}_0 = 1$.
Case 1: $\mathcal{R}_0 < 1$. Since $S(t) \le S^*$, we have $\frac{S(t)}{S^*} \mathcal{R}_0 < 1 $. The condition $\frac{S(t)}{S^*} \mathcal{R}_0 = 1$ is impossible. Therefore, the only possibility for $\dot{L}(\mathbf{x})=0$ is that $\mathbf{w}^\top \mathbf{x} = 0$. Since $\mathbf{w}$ and $\mathbf{x}$ are non-negative vectors and $\mathbf{w}$ is non-zero, this forces $\mathbf{x}=\mathbf{0}$ (i.e., $I_1=I_2=J=A=0$). If the system is in a state where all infectious compartments are zero, the first ODE becomes $\dot{S} = \mu - \nu S$, whose only equilibrium is $S(t) = \mu/\nu = S^*$. Thus, the only invariant set is the DFE, $e_1$.
Case 2: $\mathcal{R}_0 = 1$. Now, $\dot{L} = 0$ implies either $\mathbf{x}(t)=\mathbf{0}$ (which leads to the DFE as above) or $S(t)=S^*$. Let's analyze the second possibility. If a trajectory stays in the set where $S(t)=S^*$ for all $t$, then we must have $\dot{S}=0$. Plugging $S=S^*$ into the first ODE gives $$0 = \mu - (\beta_1 S^* I_2 + \beta_2 S^* J + \beta_3 S^* A) - \nu S^*.$$
Since $\mu = \nu S^*$, this simplifies to $$\beta_1 S^* I_2 + \beta_2 S^* J + \beta_3 S^* A = 0.$$
As all parameters and variables are non-negative, this forces $I_2 = J = A = 0$. With these compartments being zero, the system for $I_1$ becomes $\dot{I_1} = -b_1 I_1$, which implies $I_1(t) \to 0$. Therefore, any trajectory must approach a state where all infectious compartments are zero. Again, the only invariant set is the DFE, $e_1$.
By LaSalle's Invariance Principle, since the only invariant set contained in $\mathcal{S}$ where $\{\dot{L}=0\}$ is the DFE, all solutions in the feasible region converge to the DFE, $e_1$, when $\mathcal{R}_0 \le 1$. Thus, the DFE is globally asymptotically stable.
P.S.: This is not perfect (I see at least 2-3 improvements possible) and there must be typos, etc, but to write it took me 2 and a half hours straight, so I will come back to it later. I need to get some fresh air.
