SIRS Model doesn't depend on initial conditions?

Question

So i have been working (as an undergrad, by working i mean "Redoing a few things my professor does") in a SIRS model for epidemies. SIRS stands here for:

Susceptible -> Infected -> Recovered -> Susceptible.

So, the system of 1st order differential equations that rule the model is: $$ S'(t) = - \beta I(t) S(t) + \mu R(t)$$ $$ I'(t) = \beta I(t) S(t) - \gamma I(t)$$ $$ R'(t) = \gamma I(t) - \mu R(t)$$ Where: $R'(t)$ stands for derivative of $R$ with respect to time $t$ (the same holds for $I$ and $S$). $\gamma, \beta$ and $\mu$ are constants, that depend on the properties of the diesase (if its very contagious, if it has a high chance of killing you, and so on).

So i intended to solve these equations using Euler's Method in Python3, here is my code:

while t<140:
  Sold = S
  Iold = I
  S = S + dt*(u*R - B*I*S)
  I = I + dt*(B*I*S - G*I)
  R = R + dt*(G*I - u*R)
  t = t + dt
  Arq.write("{} {} {} {} {} \n".format(t,S,I,R,R+I+S))

Notice that there is a definition of old variables. I was using them to save the previous value of the function before using it in the other equation. (I mean, instead of updating $R$ using the new value of $I$, use the old value of $I$, before being updated in the line above). I stopped using the old variables and nothing changed graphically speaking. If you want the code to run in using old variables, here it is:

while t<140:
  Sold = S
  Iold = I
  S = S + dt*(u*R - B*I*S)
  I = I + dt*(B*I*Sold - G*I)
  R = R + dt*(G*Iold - u*R)
  t = t + dt
  Arq.write("{} {} {} {} {} \n".format(t,S,I,R,R+I+S))

Okay. So let's now talk about what is happening. I'm using a 'normalized' population, i.e, S + I + R = 1 everytime. So, i choosed arbitrary values between 0 and 1 for the constants $\mu, \gamma$ and $\beta$, and set the following initial values for $I,R$ and $S$: $I(0) = 0.1$, $R(0) = 0.05$ and $S(0) = 0.85$. I got the following result, in a graphic population $\times$ time:

(I'm sorry, i don't know how to reduce the image size here). So now, i'll change the initial values to very different values: $S(0) = 0.20$, $R(0) = 0.55$, and $I(0) = 0.25$. Here is the result: And what we can see here is that, despite the initial values being very different, the two cases seem to converge to the same values of the three populations. I won't put other images here otherwise this question would be a kilometer long, but all the initial conditions i've put didn't changed the result. (or didn't seem to).

Why is this? Is my model correct?If not, what is the problem with it? If it is correct, can this be explained mathematically? This is a model that tries to simulate diseases. Are there any example of a non-killing disease that had (or have) this property?

Thanks in advance. (I would like to add that an edit to add colors to the code and reduce the image's size would be very appreciated (especially if the one who did it explained how! :))

Answer 1

Your model seems to have a single stable critical point that attracts all orbits. You're plotting $(S(t), I(t), R(t))$ as functions of time for specific initial conditions, but actually the easier way to see this directly is to plot the vector field $(\dot S, \dot I)$ in the $S$-$I$ plane like so:

Immediately you can see that for this particular choice of parameters (I don't remember which values I used), every solution of the ODE will converge to the same critical point.

To analyze this behaviour mathematically, you first identify the location of the critical point by solving $\mathrm{rhs}(S_0,I_0,R_0) = 0$ ($R=1-S-I$), and then linearize the ODE around that critical point and look at the eigenvalues of the Jacobian of the r.h.s. This doesn't quite show that all orbits converge to the critical point, only those that start close enough to it. To show that all orbits converge to the same critical point you could look for a Lyapunov function for your system.

These are all very standard topics discussed in any ODE or dynamical systems textbook. One recent book that I liked was Exploring ODEs by Trefethen, Birkisson, Driscoll (free online), look at Chapters 15–16.

I stopped using the old variables and nothing changed graphically speaking.

Imagine you're solving $\dot y = f(y)$ and you make a small error and you accidentally end up solving $\dot y = f(y) + \delta(t)$. In the case of $\dot I$, $\delta$ would be $\beta I (S(t_{n+1}) - S(t_n))\,\delta t$, although this expression isn't very well defined mathematically. This $\delta$ is merely a perturbation of size $|\dot S|\delta t^2$. If you used a higher-order method that forward Euler you would identify this by noticing that your method doesn't converge with the right order. Instead of accurately solving the right ODE, you (accurately) solve a slightly wrong ODE, and the difference for your model cannot be seen because it just converges to the same critical point. Other ODEs with higher-order methods wouldn't be this forgiving.