This page is more advanced than the previous, and is intended to support students and teachers working with the text Modeling Life (Springer Nature)

It focuses on some simpler epidemiologic models, and studies them with the techniques of nonlinear dynamics: the existence of Equilibrium Points and the analysis of their stability and instability by means of simulations, nullclines, and Linear Stability Analysis

Here we will study 3 Models

• SIR Model

  A 3 variable SIR model obtained from the previous 5 variable model by collapsing all Infecteds (\(E, I , M\)) into one category (\(I\)).

• SI Model

  A 2 variable model obtained from the 3 variable model by ignoring recovery (\(R\)) from the infection, and

• SI-density Model

  A 2 variable model with a density dependent probability of encounter.

Let’s start with a model with fewer variables, so we can study the dynamics more easily.

• We will ignore the exposed/infected/infectious distinctions, and call them all “infected” (\(I\)).
• We will also assume that people are injected into the susceptible population at a rate \(\sigma\), that the overall population death rate is \(d\), and that infected people have a higher death rate \( \delta \) (\(\delta > d\)).

\[ \left\lbrace \begin{aligned} S' &= -c\beta\dfrac{I}{S+E+I+M+R}S \\[1em] \color{red}E' &= c\beta\dfrac{I}{S+E+I+M+R}S - \alpha E \\[1em] \color{red}I' &= \alpha E - \gamma I \\[1em] \color{red}M' &= \gamma I - w M \\[1em] R'&= w M \end{aligned} \right. \]

\[ \left\lbrace \begin{aligned} S' &= {\color{blue}\sigma} - dS - c\beta\dfrac{SI}{S+I+R} \\[1em] \color{red}I' &= c\beta\dfrac{I}{S+I+R}S - (\delta+v)I \\[1em] R' &= vI - dR \end{aligned} \right. \]

Parameter	Short Description	Default Value	Unit
\( \sigma \)	population growth rate	10	person(s)/day
\( d \)	per capita normal death rate	0.01	1/person/day
\( c \)	encounter rate	1	encounter(s)/day
\( \beta \)	probability of infection per encounter per day	0.5	1/encounter/day
\( \delta \)	per capita infected death rate	0.02	1/person/day
\( v \)	per capita infection recovery rate per day	0.7	1/person/day

1-6 Parameters

The first step in understanding a differential equation is, as always, to find its equilibrium points. In epidemiology, we are interested in two kinds of equilibrium points:

• a Disease-Free Equilibrium Point, at which the infected population \(I=0\),
• and possibly other equilibrium points, at which \(I \neq 0\). Such points are called endemic equilibrium points (or outbreak equilibrium points).

Clearly, there is a Disease-Free Equilibrium (DFE) Point at \[ (S, I, R) = (\dfrac{\sigma}{d}, 0 ,0) = (\dfrac{10}{0.01}, 0 ,0) = (1000, 0 , 0) \] We can see this by realizing that \(I = 0\) makes \(I' = 0\). Then if \(R' =0\), \(R\) must equal 0 also. Plugging in \(I = 0\), \(R =0\) into the \(S' = 0\) equation leaves \(S = \frac{\sigma}{d}\). So \( (S, I, R) = (\frac{\sigma}{d}, 0 ,0) \) is indeed an Equilibrium Point.

For those parameter values, the only equilibrium point is the Disease-Free. But for other parameter values, there can be Endemic Equilibrium Points. We will experiment here with a highly infectious agent, for which beta, the probability of infection per encounter, could be as high as 0.9.

With new parameter values, let's look for an endemic equilibrium, for which the infected population \(I\) does not equal to \(0\).

One way to find the Outbreak Equilibrium Point is to calculate the nullclines, plot them and see whether they intersect. The two nullclines for \(S' = 0\) and \(I' = 0\), calculated above, are: \[ \left\lbrace \begin{aligned} S &= \dfrac{1}{d}\left(\sigma - \dfrac{\beta SI}{S+I+\frac{r}{d}I}\right) \\[1em] \dfrac{\beta S\bcancel{I}}{I+S+\frac{r}{d}} &= (\delta + v)\bcancel{I} \\ \end{aligned} \right. \]

We would like to know if these two nullclines meet, for if they do, that would be a non-trivial Outbreak Equilibrium. Let's plot these two curves in \((S, I)\) space.

Here we used Python to plot the two curves, and we see that they intersect.

We can also use Python (or any other programming language) to directly calculate the equilibrium points to be at roughly \( S =780, I = 3\), which implies that \(R = \frac{vI}{d} = 197\).

As always, the next question is, are the equilibrium points stable or unstable?

We have several ways of determining the stability of equilibrium points; the two main ones are simulations and linear stability analysis.

Simulation tells us that either of these EPs can be stable, for different parameter values.

If \(c\) or \(\beta\) is very low, the DFE is stable: if \(\beta = 0.5\), the simulation tells us that there is indeed a stable disease free EP.

Note that \(R\) and \(I\) are indeed going to zero, but \(R\) is changing very slowly. With these parameters, it takes 1 year for \(R\) to recover completely, although \(I\) goes to zero within days.

Trying several different initial conditions, in a state space trajectory plot, demonstrates that the point \((1000, 0, 0)\) is indeed a stable equilibrium point for the system.

If \(\beta\) is high, the Endemic Equilibrium Point is stable: We saw that a higher value of \(\beta\), such as \(\beta = 0.9\), can create an outbreak equilibrium. A simulation shows a damped oscillation in an approach to the Endemic Equilibrium point, roughly \[S = 780, I = 3, R = 197\] Compare this simulation result with the numerically calculated values above. They are very close.

The time series above shows an oscillatory approach to the equilibrium point, which is confirmed by the trajectory in the state space below.

The simulations confirm that the Endemic Equilibrium is stable, with a damped oscillatory approach.

Note that the model predicts oscillations or waves of infection simply in response to the dynamics. These waves of infection are a source of intense public interest, because policy makers are struggling to deal with a second wave of coronavirus outbreak.

And note that these waves are inherent in the rebound dynamics of this SIR model. They do not include any policy changes made during the course of the epidemic. Current interest is focusing on secondary waves of infection as a result of changes in social contact rules (altering \(c\) or \(\beta\)), but those waves will be coupled to any waves inherent in the rebound dynamics.

The news these days frequently mentions the concept of R₀ (pronounced “R-zero”(mostly American) or “R-naught” (mostly British)). R₀ attempts to quantify, in a single number, whether an infection will spread (R₀ > 1) or die out (R₀ < 1).

The concept is explained very nicely by Kate Winslet in the movie Contagion

How would we calculate such a number? Rob J de Boer’s text Biological Modeling of Populations gives 3 ways to calculate R₀.

Recall that

So when \(R' = 0\), that implies that \(R = \dfrac{v}{d}I\). we will use this to substitute for \(R\) in the following.

Then let’s ask what it means for there to be an EE (Endemic Equilibrium).

If \(I' = 0\), then \[ c\beta\dfrac{I}{S+I+R}S =(\delta+v)I \] assuming \(I > 0\), we can divide both sides by \(I\), to get \[ c\beta\dfrac{\bcancel{I}}{S+I+R}S =(\delta+v)\bcancel{I} \implies c\beta\dfrac{S}{S+I+R} = \delta+v \] now we will substitute \(\dfrac{v}{d}I\) for \(R\), to get \[ c\beta \dfrac{S}{S + I + \frac{v}{d}I} = \delta + v \] rearranging terms \[ I = \dfrac{1}{1+\frac{v}{d}}({\color{red}c\beta \dfrac{S}{\delta + v}} - {\color{blue}S}) \] since \(I\) must be > 0, that says that \({\color{red}c\beta \dfrac{S}{\delta +v}} > {\color{blue} S} \) , which dividing by \(S\) gives \[ {\color{red}c\beta \dfrac{\bcancel{S}}{\delta +v}} > {\color{blue} \bcancel{S}} \implies c\beta\dfrac{1}{\delta +v} > 1 \] \(\frac{c\beta}{\delta +v}\) is the quantity called R₀, and so we’ve shown that if there’s an Endemic Equilibrium, then R₀ > 1.

If we calculate R₀ for our example,

The best way to look at R₀ is to see that R₀ is the ratio of new cases/existing case. Then, clearly, if R₀ > 1, the epidemic will grow, and if R₀ < 1, it will decline.

Since R₀ = new cases/existing case, the number of new cases is equal to new case rate\(\times\) \(S_0\) \(\times\) infectivity duration That’s how many people will be infected during the infectivity lifetime, calculated at the initial population \(S = S_0\).

The infectivity duration is the inverse of the infection clearance rates. There are 2 infection clearance rates, \(\delta\) and \(v\), so the inverse of the sum of these \(\frac{1}{\delta + v}\) is the infectivity time.

So the ‘new case rate’ at \(S = S_0\) is equal to \(\beta\times\)S₀, and over \(\frac{1}{\delta + v}\) days, it amounts to a total infected of

\[ \dfrac{c\beta\times S_0}{\delta +V} \] We then calculate the per capita “new cases per old case” by dividing by S₀, and once again arrive at This time as the meaning of new cases/existing case.

---

Disease-Free Equilibrium

If the Disease-Free Equilibrium is unstable, it follows that the population must go to a non-zero equilibrium, that is, the EE.

Now, it should be noted that “DFE unstable” is not the only way for a deadly epidemic to spread. There are 2 kinds of cases where the DFE is stable, and yet epidemics occur.

• The first is when the DFE is stable, but for certain initial conditions, there is a large excursion away from the equilibrium point, followed by a gradual excursion back to the EP. That large transient excursion could kill millions of people.
• The second is when the DFE is stable, but there is another endemic equilibrium that is also stable. Then a sufficient perturbation could kick the system out of the basin of the DFE and into the basin of the EE.
• We are interested here in the third way an Epidemic can occur, which is by the DFE becoming unstable.

The definition of R₀ can be put very simply in nonlinear dynamics:

R₀ is the linear approximation to the nonlinear SIR equations at the disease-free equilibrium point \((S = S_0, I = 0, R = 0)\).

The criterion (R₀ > 1) therefore guarantees that the disease-free state is unstable, and the infection will grow.

We can therefore derive R₀ for a particular model by

• finding the Equilibrium Points of the model
• using the Jacobian of the vector field to do the linear stability analysis at those EPs.

Let's go back to our 'high-\(\beta\)' model. We saw by simulation that the DFE was unstable, and that there was an Endemic Equilibrium that appeared to be a stable EP with an oscillatory approach.

If we use Python to directly calculate the eigenvalues of these two EPs, the results exactly confirm the simulations. For the DFE, (\(\frac{\sigma}{d}\), 0 , 0), Python calculates the Jacobian matrix to be

\[ \displaystyle \mathbf{J} = \left[\begin{matrix}- d & - \beta c & 0\\0 & \beta c - \delta - v & 0\\0 & v & - d\end{matrix}\right] \]

We see that only the main diagonal contributes to the eigenvalue calculation, since all the other products are 0. Since the product of three numbers can only be zero when at least one of them is zero, we conclude that there are 3 eigenvalues: 2 of them are negative real (\(-d\)) and the third is equal to \(c\beta - \delta- v\). \[ \begin{aligned} \lambda_1 &= -d \\ \lambda_2 &= -d \\ \lambda_3 &= c\beta - \delta - v \end{aligned} \]

Therefore, the Disease Free Equilibrium is unstable if and only if \[ c\beta-\delta-v > 0 \] Rearranging terms, this gives us that the DFE is unstable if and only if that eigenvalue is positive, that is \[ \dfrac{c\beta}{\delta + v} > 1 \] which confirms our calculation of R₀, now as the criterion for the stability of the Disease Free EP.

---

Endemic Equilibrium

If we try calculate the eigenvalues of the Outbreak (Endemic) Equilibrium, we would see that the symbolic calculation is huge and messy algebra, so we will only calculate it numerically, for these parameter values. They are:

\[ \begin{aligned} \lambda_1 &= -0.01 \\ \lambda_2 &= -0.0063+0.042i \\ \lambda_3 &= -0.0063-0.042i \end{aligned} \]

The eigenvalues are therefore one negative real, and the other two are a pair of complex conjugate eigenvalues with negative real part, indicating that the Endemic Equilibrium is indeed a stable EP with a damped oscillatory approach, exactly what we saw in the simulation.

Question
Note that the real eigenvalue has a very small value, and so do the real parts of the imaginary eigenvalues. What parameters cause this? and what are the consequences of this for the behavior of the model?

---

Bifurcation Diagram

Finally, we can summarize our results in a bifurcation diagram. By running hundreds of simulations for varying values of \(\beta\), we arrive at the bifurcation diagram below. (Recall R₀ = \(\frac{c\beta}{\delta+v}\).)

It helps to have a 2-variable SI model to study qualitative dynamics. We can reduce our model above to a 2-variable model by assuming that people do not recover from the infection; they stay infected and infectious for ever. This amounts to assuming that the recovery period is long with respect to the time scale we are interested in.

With that assumption, \(v = 0\) and the new equations become

\[ \begin{aligned} S' &= \sigma - dS - c\beta\dfrac{SI}{S+I} \\[1em] I' &= c\beta\dfrac{SI}{S+I} - \delta I \end{aligned} \]

Parameter	Short Description	Default Value	Unit
\( \sigma \)	population growth rate	10	person(s)/day
\( d \)	per capita normal death rate	0.01	1/person/day
\( c \)	encounter rate	1	encounter(s)/day
\( \beta \)	probability of infection per encounter per day	0.5	1/encounter/day
\( \delta \)	per capita infected death rate	0.02	1/person/day

1-6 Parameters

Our first task is, as ever, to find the Equilibrium Points

Disease-Free Equilibrium Clearly, there is a Disease-Free Equilibrium (DFE) at \(S = \frac{\sigma}{d}, I = 0\).

Endemic Equilibrium To find a second equilibrium point at which \(I \neq 0\) (an Endemic Equilibrium), we must look at the nullclines \(S' = 0\) and \(I'= 0\) and see whether they intersect.

Let's calculate the nullclines for \(S'=0\) and \(I'=0\).

Therefore, the \(S'=0\) nullcline is

\[ I' = 0 \qquad \implies c\beta \dfrac{SI}{S+I} = \delta I \\[1.5em] \]

Since \(I \neq 0\), we can divide both sides by \(I\),

\[ \begin{aligned} &&c\beta \dfrac{S\bcancel{\color{red}I}}{S+I} &= \delta \bcancel{\color{red}I} \\[1.5em] &\implies& c\beta \dfrac{S}{S+I} &= \delta \\[1.5em] &\implies& c\beta S &= \delta (S+I) = \delta S + \delta I\\[1.5em] &\implies& I &= \dfrac{c\beta S - \delta S}{\delta}\\[1.5em] \end{aligned} \]

(Note that since \(I > 0\) then \(c\beta > \delta\) which means \(\frac{c\beta}{\delta} > 1\). This is R₀ for this model.)

We now find the Endemic Equilibrium Point by plotting the two nullclines and seeing where they cross.

The next task is to determine the stability of the equilibrium points. Our first method will be to use simulation.

Let's run a simulation of the model with a low value of beta, say \(\beta = 0.015\).

We see that the Disease-Free Equilibrium, \(S = \frac{\sigma}{d} = 1000\), \(I = 0\) is stable.

But if we raise \(\beta\) to 0.03, now the simulation shows a stable Endemic Equilibrium Point at as predicted.

By making a large number of simulations, varying \(\beta\), we can form the Bifurcation Diagram for the 2 variable model:

We see that, as R0 passes 1, the Disease-Free equilibrium point becomes unstable, and the Endemic Equilibrium becomes stable.

But we are not yet sure what type of bifurcation this is. For that, we have to extend the analysis of the model to include negative solutions for I. Of course all negative values for S and I are impossible; populations can't be negative. But including the impossible negative solutions give us the answer to our question:

Clearly, the bifurcation that occurs as R₀ passes 1 is a transcritical bifurcation(see Modeling Life p 156-158), with the Disease-Free EP and the Endemic EP exchanging stability as R₀ passes 1.

Exercise
Verify all of these statements by carrying out the eigenvalue analysis of the stability of the Disease-Free Equilibrium Point. This can be done by hand since there are only 2 variables in this model.

Finally, we want to consider a slightly simpler model. In this model, we will make the susceptible meets infected term to be, not \(c\beta\frac{SI}{S+I}\) as previously, but rather, just \(c\beta SI\).

What could this be modeling? We derived the original model by assuming that an \(S\)(susceptible) bumps into a random person. The probability of this random person being an infected is \(\frac{I}{S+I}\), the probability of that infected person infecting the susceptible is \(\beta\), and the number of encounters/day is \(c\), so that gave us the key infection term as \(c\beta\frac{SI}{S+I}\).

But notice that we are assuming that \(c\) is constant. Is that correct? de Boer points out (personal communication) that ‘c = constant’ is a good model for sexually transmitted diseases, where ‘contact’ means ‘sexual contact’. In that case, \(c\) is a choice, how many sexual contacts per day, and does not go up and down with population density.

On the other hand, for a flu-like disease like COVID-19, the number of contacts/day goes up with \(N = S+I\). If \(c = k(S+I)\), then the equations become \[ \begin{aligned} S' &= \sigma – dS – k \beta SI \\[1em] I' &= k\beta SI – \delta I \end{aligned} \]

Parameter	Short Description	Default Value	Unit
\( \sigma \)	population growth rate	10	person(s)/day
\( d \)	per capita normal death rate	0.01	1/person/day
\( k \)	density dependence	1
\( \beta \)	probability of infection per encounter per day	0.5	1/encounter/day
\( \delta \)	per capita infected death rate	0.02	1/person/day

1-6 Parameters

The disease free equilibrium point is easily found by setting \(I = 0\), and calculating \(S = \frac{\sigma}{d}\). So \((\frac{\sigma}{d}, 0 )\) is our DFE.

We can find the Endemic Equilibrium by calculating the nullclines and seeing whether they intersect.

The \(S'=0\) nullcline is \(I = \dfrac{\sigma – dS}{k\beta S}\)

The \(I' =0\) nullcline is \(I = \dfrac{\delta}{k\beta S} \)

Exercise
Plot the nullclines for varying values of the parameters to determine the existence and values of any Endemic Equilibrium Points.

We can also approach R₀ from the point of view of the stability of the DFE. The Jacobian of the model, evaluated at the EP \((\frac{\sigma}{d}, 0)\), is

\[ J = \left[\begin{matrix} -d & -\beta S \\ 0 & \beta S - \delta \end{matrix}\right] \]

Since the Jacobian matrix has a 0 in the bottom left entry, the eigenvalue calculation is particularly easy: the product of 2 numbers is zero only when at least one of them is 0, and so

\[ (-d-\lambda)\times(\beta S – \delta-\lambda) = 0 \implies \lambda_1 = -d, \lambda_2 = \beta S - \delta. \] Now, \( \lambda_2 = \beta S – \delta\) is > 0 only when \(\beta S > \delta\), or when \(\frac{\beta S}{\delta} > 1\).

This is R₀.

Exercise
Plot the location and stability of the EPs for this model in a bifurcation diagram.