# Selecting a suitable Lyapunov function for the following systems?

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}

\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).