Dynamical Systems Analysis of Epidemiological Models

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.

SIR Model

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

Old Equations
\[ \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. \]
New Equations
\[ \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

Finding Equilibrium Points

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

Disease-Free Equilibrium Point

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.

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

\[ \left\lbrace \begin{aligned} R' &= 0 \\ S' &= 0 \\ I' &= 0 \\ I &\neq 0 \\ \end{aligned} \right. \implies \left\lbrace \begin{aligned} R &= \dfrac{r}{d}I \\[1em] 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. \]

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.

\[ (S, I, R) = (780, 3, 197) \] Parameters
\[ \begin{aligned} \sigma &= 10 \\ d &= 0.01 \\ c &= 1\\ \beta &= 0.9\\ \delta &= 0.02 \\ v &= 0.7 \end{aligned} \]

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

Stable or Unstable?

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.

\[ (S, I, R) = (1000, 0, 0) \] Parameters
\[ \begin{aligned} \sigma &= 10 \\ d &= 0.01 \\ c &= 1\\ \color{red}\beta &\color{red}= 0.5\\ \delta &= 0.02 \\ v &= 0.7 \end{aligned} \]

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.

\[ (S, I, R) = (780, 3, 197) \] Parameters
\[ \begin{aligned} \sigma &= 10 \\ d &= 0.01 \\ c &= 1\\ \color{red}\beta &\color{red}= 0.9\\ \delta &= 0.02 \\ v &= 0.7 \end{aligned} \]
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.

“R-zero” or “R-naught”

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

(#1) Derive the condition for the existence of an Endemic Equilibrium (by definition with \(I > 0\)).
Recall that
\[ \left\lbrace \begin{aligned} I' &= c\beta\dfrac{I}{S+I+R}S - (\delta+v)I \\[1em] R' &= vI - dR \end{aligned} \right. \]

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 R0, and so we’ve shown that if there’s an Endemic Equilibrium, then R0 > 1.

If we calculate R0 for our example,

R0 = \(\dfrac{c\beta}{\delta +v} = \dfrac{1\times 0.9}{0.02+0.7} = 1.25 \)

(#2) The second method, which is widely used, makes explicit that R0 is a linear approximation.

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

Since R0 = 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\)S0, 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 S0, and once again arrive at
R0 = \(\dfrac{c\beta}{\delta +v}\)
This time as the meaning of new cases/existing case.

(#3) The third way to calculate R0 is to say that we are looking for cases where the Disease-Free Equilibrium is unstable.

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 R0 can be put very simply in nonlinear dynamics:

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

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

We can therefore derive R0 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 R0, 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.

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 R0 = \(\frac{c\beta}{\delta+v}\).)

SI model

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

Finding Equilibrium Points

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

\[ \begin{aligned} S' = 0 \qquad &\implies& (\sigma - dS) &= c\beta\dfrac{SI}{S+I} \\[1.5em] &\implies& (S+I)(\sigma – dS) &= c\beta SI \\[1.5em] &\implies& S(\sigma – dS) + I (\sigma- dS) &= c\beta SI \\[1.5em] &\implies& S(\sigma – dS) &= c\beta SI - I (\sigma- dS) \\[1.5em] &\implies& S(\sigma – dS) &= I\bigg(c\beta S-(\sigma-dS)\bigg) \\[1.5em] \end{aligned} \]
Therefore, the \(S'=0\) nullcline is
\[ I = \dfrac{S(\sigma – dS)}{c\beta S-(\sigma-dS)} \]
\[ 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 R0 for this model.)

R0 = \(\dfrac{c\beta}{\delta}\)

Finding the Endemic Equilibrium Point

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

\[ (S, I) = (20, 490) \] Parameters
\[ \begin{aligned} \sigma &= 10 \\ d &= 0.01 \\ c &= 1\\ \beta &= 0.5\\ \delta &= 0.02 \end{aligned} \]

Stable or Unstable?

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.

Bifurcation Diagram

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 R0 passes 1 is a transcritical bifurcation(see Modeling Life p 156-158), with the Disease-Free EP and the Endemic EP exchanging stability as R0 passes 1.

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.

SI model - density dependent \(c\) (de Boer)

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

Finding Equilibrium Points

Disease-Free Equilibrium (DFE)
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.

Endemic Equilibrium (EE)

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

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

Linear Stability Analysis for the DFE

We can also approach R0 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 R0.

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