3
$\begingroup$

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: enter image description here

(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: enter image description here 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! :))

$\endgroup$
  • 3
    $\begingroup$ Your model seems to have a single stable critical point that attracts all orbits (example vector field plot of $(\dot s,\dot i)$). Do you know how to analyze autonomous ODEs by linearizing them around a critical point? Also, I suspect the reason it doesn't matter if you use the right values during integration is that it only introduces a small perturbation, and that perturbation decays anyway. $\endgroup$ – Kirill Oct 22 '17 at 21:55
  • $\begingroup$ I don't know what you mean by "Linearize ODEs around a critical point". I think i understand the concept of a stable critical point, but don't know what is linearize ODEs around a critical point. $\endgroup$ – Vitor C Goergen Oct 22 '17 at 23:39
  • $\begingroup$ @kirill -- I think you should make this a regular answer since that's the correct answer to the question. $\endgroup$ – Wolfgang Bangerth Oct 23 '17 at 14:01
3
$\begingroup$

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:

enter image description here

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.

$\endgroup$
  • $\begingroup$ Interesting. So it seems like we are dealing with a problem on the model itself. Is there any way of changing the equation in order to avoid this? I also thank you for the book recommendation. I'll take a look. Thanks! $\endgroup$ – Vitor C Goergen Oct 23 '17 at 21:31
  • 1
    $\begingroup$ It's not a problem with the equation, just a property of the equation in conjunction with the values you chose for the coefficients. $\endgroup$ – Wolfgang Bangerth Oct 24 '17 at 3:49
  • $\begingroup$ @VitorCGoergen Yeah, it's not a problem per se, at some point you have to ask what it is that you're modelling, and whether this particular behaviour represents that. $\endgroup$ – Kirill Oct 24 '17 at 20:04

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.