# Selecting a Lyapunov function for a SEI model

Consider the system:

$$\frac{dS}{dt} =\nu N -\frac{\beta S I}{N} -\nu S$$ $$\frac{dE}{dt} =\frac{\beta S I}{N} -(\sigma+\nu)E$$ $$\frac{dI}{dt} = \sigma E -\nu I$$

where $$N=S+E+I$$.

Removing degeneracy we have the following reduced system:

$$\frac{dS}{dt} = -\beta S I +\nu -\nu S$$ $$\frac{dI}{dt} = \sigma -\sigma S -(\sigma+\nu) I$$

where $$\mathcal{R}_0 = \frac{\beta \sigma}{\nu(\nu+\sigma)}$$.

The equilibrium point:

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

I've been looking for a (global) Lyapunov function to prove global stability for the equilibrium point $$e_2$$ but I haven't found one that works with this specific system, Any ideas?

A similar question was asked here: Searching for a Lyapunov function for a SIRS model

• I am looking for a canonical answer.
– user644376
Commented Oct 28, 2021 at 15:00
• I don't know exactly what you mean by canonical, but please see my answer below. I only saw your post recently. There a few minor details to fix, but the time is running short. I hope this helps. Commented Nov 2, 2021 at 12:56

I believe there is no Lyapunov functional for this ODE system. Nevertheless, one can prove the global stability of this fixed point by using information one gathers from nullclines of this system, which are the curves where $$\dot S=0$$ and $$\dot I=0$$. From the reduced system, we have the $$\frac{dS}{dt}=0$$ nullcline given by $$-\beta S I +\nu -\nu S=0,$$ this corresponds to the curve $$I=\frac{\nu}{\beta}\left(\frac1S-1\right).$$ Above it $$\frac{dS}{dt}<0$$ while below it $$\frac{dS}{dt}>0$$.

The nullcline with $$\frac{dI}{dt}=0$$ is given by $$\sigma -\sigma S -(\sigma+\nu) I=0.$$ The straight line is $$S+\left(1+\frac\nu\sigma\right)I=1.$$ Below it: $$\frac{dI}{dt}>0$$, while above it: $$\frac{dI}{dt}<0$$. Due to the normalization $$S+E+I=1$$, hence, all states start in the region $$\Omega$$ (greenish shaded area), which we define as below the line $$S+I=1$$ and in the first quadrant.

Defining the velocity of the flow given in the reduced system as $$\vec v=(-\beta S I +\nu -\nu S, \sigma -\sigma S -(\sigma+\nu) I),$$ On the boundary $$I=0$$ and $$0\leq S\leq1$$, one finds that $$\frac{dS}{dt}=\nu(1-S)\geq 0$$ and $$\frac{dI}{dt}=\sigma(1-S)\geq0$$, hence the trajectories cannot leave the region $$\Omega$$ there. On the boundary $$S=0$$ and $$0\leq I\leq1$$, one finds that $$\vec v\cdot(1,0)=\nu\geq0$$, hence the trajectories cannot leave the region $$\Omega$$ there either. Finally, on the boundary $$S+I=1$$, whose normal vector is $$\hat n=(\frac1{\sqrt{2}}, \frac1{\sqrt{2}})$$, $$\vec v\cdot \hat n=-\frac1{\sqrt{2}}\beta SI\leq0$$. Hence, again no trajectories can leave region $$\Omega$$ across this boundary. Therefore, we conclude that $$\Omega$$ is an invariant region under the flow defined by the reduced ODE system.

One can show that the shaded blue area between the nullclines is an attractor. We call it attractor $$A$$. Furthermore, there are no other fixed points in $$A$$ besides the unstable fixed point $$e_1=(1,0)$$ and the stable fixed point $$e_2=\left(\frac{\nu\left(\sigma + \nu\right)}{\beta \sigma}, \frac\sigma{\sigma+\nu} -\frac\nu\beta \right)$$. Also, since the trajectories cannot leave region $$A$$ since on the nullcline 1 (lower boundary of $$A$$) the flow points upwards and on nullcline 2 (upper boundary of $$A$$) the flow points leftwards. In addition, in region $$A$$ $$\frac{dS}{dt}\leq0$$ and $$\frac{dI}{dt}\geq0$$ for all time, then periodic solutions (limit cycles) are ruled out. Furthermore, due to the Poincaré-Bendixon theorem (PBT), all trajectories in that region with converge to $$e_2$$. We don't even need to use the PBT inside region $$A$$, if we notice that in each vertical cut of it, the flow points from left to right towards fixed-point 2. Take a cut at $$S=S_0$$ for example, then we obtain $$\vec v\cdot (-1,0)=\beta S_0 I -\nu +\nu S_0\geq0$$. See the trajectories in Fig. 2, starting at one such cut, illustrating this point.

Below nullclines 1 and 2, we have $$\frac{dS}{dt}>0$$ and $$\frac{dI}{dt}>0$$, hence no fixed-points and no limit cycles exist there. In region $$\Omega$$, below nullcline 1 and above nullcline 2, we have $$\frac{dS}{dt}>0$$ and $$\frac{dI}{dt}<0$$, hence, there are no fixed-points nor limit cycles therein. Finally, in region $$\Omega$$ above both nullclines, we have $$\frac{dS}{dt}<0$$ and $$\frac{dI}{dt}<0$$, hence there are no fixed points nor any periodic orbits there.

• Firstly, thank you for your answer! How do you know there is no Lyapunov function(s) for this system? Also, if you don't mind, can you post your code for the figure? I am also stuck proving global stability for another but similar model, maybe you can have a look when I post it?
– user644376
Commented Nov 3, 2021 at 14:25
• I can't prove there isn't a LF. I suspect there isn't one. But if your main objective is to prove global stability, then you don't need it. At least, for the present case. Also, one can still improve on my answer. For example, you can take vertical cuts in the region A (between the nullclines) and you can easily find that all trajectories crossing the cut flow from right to left. Then, inside the region A, you wouldn't even need the Poincaré-Bendixon theorem. Commented Nov 3, 2021 at 16:32
• Before I post the code, let me know if the answer I provided is acceptable or not to you. Commented Nov 3, 2021 at 16:35
• As long as it proves global stability, which it does(I need to do a deeper read tomorrow) then, yes it is acceptable. I will respond back to your response tomorrow.
– user644376
Commented Nov 3, 2021 at 16:56
• I just added a few comments and a new figure displaying trajectories to better illustrate my point. Notice also, that with this construction of nullclines, when $R_0\rightarrow1$, $e_2\rightarrow e_1$ and you have a bifurcation. Commented Nov 3, 2021 at 19:25