**i) SI MODEL**

Consider

\begin{align} \frac{dS}{dt} &= \mu N -\frac{\beta S I}{N} - \nu S\\[2ex] \frac{dI}{dt} &= \frac{\beta S I}{N} -\nu I \end{align} Where $N=S+I$ is the total population.

If we assume $\mu =\nu$, the above reduces to:

\begin{align} \frac{dS}{dt} &= -\beta S I + \nu I\\[2ex] \frac{dI}{dt} &= \beta S I -\nu I \end{align}

The equilibrium points:

\begin{align*} e_1 : \left( S_1^*, I_1^*\right)&= \left(1, 0\right), \\[2ex] e_2 : \left( S_2^*, I_2^*\right)&= \left(\frac{\nu }{\beta}, \frac{\nu}{\beta}\left(\frac{\beta}{\nu}-1 \right)\right) \end{align*}

where $\mathcal{R}_0 = \beta/\nu$.

The set $\Omega = \left\lbrace \left(S,I\right)\in \mathbb{R}_+^2 : S\geq 0, I \geq 0, S+I \leq 1 \right\rbrace$ is our domain of definition. This set is a positively invariant set for our system.

Now to prove the global (asymptotically) stability of $e_1$ is straightforward, we use the function $V = I$ and the result follows.

Now my question is, how would I construct a Lyapunov function for $e_2$? I know it should be in the form of a Volterra function but, how can I choose a particular one?

**ii) SIS MODEL**

Consider

\begin{align} \frac{dS}{dt} &= \mu N -\frac{\beta S I}{N}+ \gamma I - \nu S\\[2ex] \frac{dI}{dt} &= \frac{\beta S I}{N} -(\gamma+\nu) I \end{align} Where $N=S+I$ is the total population.

If we assume $\mu =\nu$, the above reduces to:

\begin{align} \frac{dS}{dt} &= -\beta S I + (\gamma +\nu) I\\[2ex] \frac{dI}{dt} &= \beta S I -(\gamma+\nu) I \end{align}

The equilibrium points:

\begin{align*} e_1 : \left( S_1^*, I_1^*\right)&= \left(1, 0\right), \\[2ex] e_2 : \left( S_2^*, I_2^*\right)&= \left(\frac{\left(\gamma+\nu\right)}{\beta}, \frac{\left(\gamma+\nu\right)}{\beta}\left(\frac{\beta}{\gamma+\nu} -1\right)\right), \end{align*}

where $\mathcal{R}_0 = \beta/(\gamma+\nu)$.

Again, as in case (i) we have the set $\Omega = \left\lbrace \left(S,I\right)\in \mathbb{R}_+^2 : S\geq 0, I \geq 0, S+I \leq 1 \right\rbrace$ as our domain of definition. This set is a positively invariant set for our system.

Analogously as earlier, to prove the global (asymptotically) stability of $e_1$ is straightforward, we use the function $V = I$ and the result follows.

How would I chose a suitable Lyapunov function for this model to analyse the endemic equilibrium $e_2$? Well, if we find a suitable Lyapunov function for case (i) then we can just replace $\nu$ with $\gamma+\nu$. However I'd like to see both computations.

**iii) SIR MODEL**

Consider

\begin{align} \frac{dS}{dt} &= \mu N -\frac{\beta S I}{N} - \nu S\\[2ex] \frac{dI}{dt} &= \frac{\beta S I}{N} -(\gamma+\nu) I\\[2ex] \frac{dR}{dt} &= \gamma I -\nu R \end{align} Where $N=S+I+R$ is the total population.

If we assume $\mu =\nu$, the above reduces to:

\begin{align} \frac{dS}{dt} &= -\beta S I + \nu I +\nu R \\[2ex] \frac{dI}{dt} &= \beta S I -(\gamma +\nu) I\\[2ex] \frac{dR}{dt} &= \gamma I -\nu R -\xi R \end{align}

Using $R=N-S-I$ to kill degeneracy, this further reduces to

\begin{align} \frac{dS}{dt} &= -\beta S I + \nu -\nu S \\[2ex] \frac{dI}{dt} &= \beta S I -(\gamma +\nu) I\\[2ex] \end{align}

The equilibrium points read

\begin{align*} e_1 : \left( S_1^*, I_1^*, R_1^*\right)&= \left(1, 0, 0\right), \\[2ex] e_2 : \left( S_2^*, I_2^*, R_2^*\right)&= \left(\frac{\left(\gamma+\nu\right)}{\beta}, \frac{\nu}{\beta}\left(\frac{\beta}{\gamma +\nu}-1\right), \frac{\gamma}{\beta}\left(\frac{\beta}{\gamma +\nu}-1\right)\right)\\[1ex] \end{align*}

where $\mathcal{R}_0 = \beta/(\gamma+\nu)$.

Again, as in case (i) and (ii) we have the set $\Omega = \left\lbrace \left(S,I\right)\in \mathbb{R}_+^2 : S\geq 0, I \geq 0, S+I \leq 1 \right\rbrace$ as our domain of definition. This set is a positively invariant set for our system.

**Theorem**

If $\mathcal{R}_0 \leq 1$, then the disease-free equilibrium $e_1$ is globally asymptotically stable in $\Omega$.

**Proof**

We use the function $V = I$ and the result follows.

**Theorem**

If $\mathcal{R}_0 > 1$, then the endemic equilibrium $e_2$ is globally asymptotically stable in the interior of $\Omega$.

**Proof**

To prove the global asymptotically stability for $e_2$, consider the Lyapunov function (in the form of a Volterra function):

$$V(S,I) = \left(S-S_2^*\right)+ \left( I-I_2^*\right) -S_2^* \ln \frac{S}{S_2^*} - I_2^* \ln \frac{I}{I_2^*}. $$

working out the time derivatives along our reduced system, we arrive to:

$$\dot V = -\beta (S_2^* -S)^2 [\frac{1}{S}] \leq 0 $$

We see $\dot V$ is negative except when $\dot V$ takes on the equilibrium values so that $\dot V =0$, so we have semi-definiteness. The largest compact invariant set in $\Omega$ so that $\dot V$ is $0$ is the singleton $\lbrace{ e_2\rbrace}$. Hence, concluding from LaSalle's invariance principle, $e_2$ is globally asymptotically stable in $\Omega$.

**iv) SIRS MODEL**

Consider

\begin{align} \frac{dS}{dt} &= \mu N -\frac{\beta S I}{N} +\xi R - \nu S\\[2ex] \frac{dI}{dt} &= \frac{\beta S I}{N} -(\gamma+\nu) I\\[2ex] \frac{dR}{dt} &= \gamma I-\xi R -\nu R \end{align} Where $N=S+I+R$ is the total population.

If we assume $\mu =\nu$, the above reduces to:

\begin{align} \frac{dS}{dt} &= -\beta S I + \nu I +(\xi+\nu) R \\[2ex] \frac{dI}{dt} &= \beta S I -(\gamma +\nu) I \\[2ex] \frac{dR}{dt} &= \gamma I -(\xi +\nu) R \end{align}

with equilibrium points:

\begin{align*} e_1 : \left( S_1^*, I_1^*, R_1^*\right)&= \left(1, 0, 0\right), \\[2ex] e_2 : \left( S_2^*, I_2^*, R_2^*\right)&= \left(\frac{\left(\gamma + \nu\right)}{\beta}, \frac{\left(\gamma+\nu \right) \left( \xi + \nu \right) \left( \frac{\beta}{\gamma+\nu} -1 \right) }{\beta\left(\gamma + \xi+\nu \right)}, \frac{\gamma \left(\gamma+\nu \right)\left( \frac{\beta}{\gamma+\nu} -1 \right) }{\beta\left(\gamma + \xi+\nu \right)}\right) \end{align*}

So my question(s) are; how do I find suitable Lyapunov functions for systems (i), (ii) and (iv) similar to my example in (iii)? Can you show the solution to systems (i), (ii) and (iv) like how I presented in (iii)?

**EDIT**:

For systems (i) and (ii), due to constant population, we can reduce it to a 1-dimensional system then perform the computations.

(i):

By substituting $S=N-I$ into (2.4) we have

\begin{align} \frac{dI}{dt} &= (\beta - \nu)I - \beta I^2 \end{align}

Solving (2.5) with initial condition $I(0)= I_0$ analytically, we have the solution to the system:

\begin{align} I(t) &= \frac{I_0 (\beta - \nu)}{\beta I_0 - e^{-\left(\beta - \nu\right)t} \left[ \beta I_0 - \left(\beta - \nu\right)\right]}.\\[2ex] S(t) &= 1-I(t). \end{align}

We can make inferences of the long term behaviour of this model by examining the possible values of $\left(\beta -\nu\right)$, that is, of course when (2.6) is feasible. We have two cases \begin{align*} \beta - \nu & < 0 \\[2ex] \beta - \nu & > 0 \end{align*}

If $\beta - \nu < 0$, then $e^{-\left(\beta - \nu\right)t} \rightarrow \infty \text{ as } t \rightarrow \infty$ hence

\begin{align} \lim_{t \to \infty} I(t) = 0. \end{align}

If $\beta - \nu > 0$, then $e^{-\left(\beta - \nu\right)t} \rightarrow 0 \text{ as } t \rightarrow \infty$ hence we have the limit

\begin{align} \lim_{t \to \infty} I(t) = \frac{I_0\left(\beta-\nu\right)}{\beta I_0} = 1-\frac{\nu}{\beta}. \end{align}

(ii):

Analogously as in section 2.1, we obtain the complete solution to our system:

\begin{align} I(t) &= \frac{I_0 (\beta - \gamma-\nu)}{\beta I_0 - e^{-\left(\beta - \gamma -\nu\right)t} \left[ \beta I_0 - \left(\beta -\gamma - \nu\right)\right]}.\\[2ex] S(t) &= 1-I(t). \end{align}

As before, we make inferences of the long term behaviour of this model by examining the possible values of $\left(\beta -\gamma -\nu\right)$, that is, of course when (2.14) is feasible. We have two cases \begin{align*} \beta - \gamma-\nu & < 0 \\[2ex] \beta - \gamma -\nu & > 0 \end{align*}

If $\beta - \gamma-\nu < 0$, then $e^{-\left(\beta - \gamma -\nu\right)t} \rightarrow \infty \text{ as } t \rightarrow \infty$ hence

\begin{align} \lim_{t \to \infty} I(t) = 0. \end{align}

If $\beta - \gamma -\nu > 0$, then $e^{-\left(\beta - \gamma -\nu\right)t} \rightarrow 0 \text{ as } t \rightarrow \infty$ hence we have the limit

\begin{align} \lim_{t \to \infty} I(t) = \frac{I_0\left(\beta-\gamma -\nu\right)}{\beta I_0} = 1-\frac{\gamma+\nu}{\beta}. \end{align}

However the above isn't proving global stability..

For system (iv), we need to follow something similar to (iii).