I have a simple SIS model with a non-constant population size due to disease-related mortality:

\begin{aligned} \frac{dS}{dt} &= a I - b S I\\\\ \frac{dI}{dt} &= b S I - a I - c I \end{aligned}

where $S$ and $I$ are the susceptible and infected classes, $a$ is the rate of transition back to the susceptible class after being infected, $b$ is the per-capita transmission rate, and $c$ is the disease-related mortality for infecteds.

Finding the equilibria of this model is stumping me for some reason. First, I factor the second equation to get:

\begin{equation} I ( bS - a - c ) = 0 \end{equation}

so $\hat{I} = 0$ or $\hat{S} = \frac{a+c}{b}$ must be accord with $\dot{I} = 0$. The first means that:

\begin{equation} \dot{S} = a(0) - b S (0) = 0 \end{equation}

and so $\hat{S}$ must be its initial condition, $S_0$.

However, there is another equilibrium where the disease spreads and dies out, meaning $\hat{S} < S_0$, which I assumed is equal to $\hat{S} = \frac{a+c}{b}$.

However, using simultation in $R$, this equilibrium isn't correct - $\frac{a+c}{b}$ is still off by a fair amount.

Can anyone see where I'm going wrong? Thanks!


For $I=0$ you can have any value of $S$ to get an equilibrium position, the equilibrium equations do not impose further conditions on $S$. If you add the stability at $I=0$, then you have to add the condition $bS−a−c<0$. Outside that region the values for $I$ will increase.

If you divide both equations, you find that $$ \frac{dS}{dI}=\frac{a-bS}{bS-a-c} $$ is a separable equation, so that $$ I+S-\frac{c}{b}\ln|bS-a| $$ is a conserved quantity. That is, for the limit with $\hat I=0$ you get the equation $$ \hat S-\frac{c}{b}\ln|b\hat S-a|=I_0+S_0-\frac{c}{b}\ln|bS_0-a| \\~\\ e^{b\hat S-a}(b\hat S-a)=e^{b(I_0+S_0)-a}(bS_0-a) \\~\\ \hat S=\frac1b\left(a+W_0(e^{b(I_0+S_0)-a}(bS_0-a))\right) $$ for the limit value $\hat S$. The sign of $bS-a$ has to stay constant, $W_0$ is the Lambert-W function.

This only works for the simplest of models, more complex models will usually not give rise to conserved quantities/constants of motion.

  • $\begingroup$ Okay, thanks. I see that at for $\hat{I}=0$ and $\hat{S}$ equal to any value, stability depends on the ratio $\frac{a+c}{b \hat{S}} > 1$, which means infections will increase before dying out. But how do I find the equilibrium value of $\hat{S}$? My simulations show that when that ratio is greater than 1, $S$ settles down to an equilibrium, but I am struggling to compute what that is. Thanks. $\endgroup$ – user_15 Nov 3 '20 at 11:10
  • $\begingroup$ Thanks for your edited response. However, I'm still not understanding how we can derive the equilibrium value of $\hat{S}$ from your answer. Is it possible? Or can we only derive the conditions for when disease will grow in the population? $\endgroup$ – user_15 Nov 3 '20 at 13:42
  • $\begingroup$ What you now have is an implicit equation. You can solve it numerically or via the Lambert-W function, which is only slightly better than giving the numerical method a name. But the value of $\hat S$ is uniquely given by that equation. $\endgroup$ – Lutz Lehmann Nov 3 '20 at 13:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.