6
$\begingroup$

I have this SEIR system of differential equations for modeling pandemics

p=10 (*total population*)
b=0.9 (*infection rate*)
d=0.5 (*rate at which infected become contagious*)
y=0.2 (*death/remove rate*)
outbreak = {s'[t] == -b*s[t]*i[t]/p, 
                e'[t] == b*s[t]*i[t]/p - d*e[t], 
                i'[t] == d*e[t] - y*i[t], 
                r'[t] == y*i[t],
                p == s[t] + e[t] + i[t] + r[t]}
initialconditions = {s[0] = p, e[0] = 1, i[0] = 0, r[0] = 0, s'[0] = 0, e'[0] = 0, i'[0] = 0, r'[0] = 0}

where p(population size), b(infection rate), d(contagious rate), y(removal rate) are constants

An example of a similar can be found here, that I'm trying to reproduce.

http://www.public.asu.edu/~hnesse/classes/seir.html

I just can't figure out how to use the NDSolve command as I keep throwing up errors.

Edit: I have this now

ClearAll["Global`*"]
outbreak = {s'[t] == -b*s[t]*i[t]/p, e'[t] == b*s[t]*i[t]/p - d*e[t], 
   i'[t] == d*e[t] - y*i[t], r'[t] == y*i[t], 
   p == s[t] + e[t] + i[t] + r[t]};
initialconditions = {s[0] == p, e[0] == 1, i[0] == 0, r[0] == 0, 
   s'[0] == 0, e'[0] == 0, i'[0] == 0, r'[0] == 0};
p = 10; (*total population*)
b = 0.9 ;(*infection rate*)
d = 0.5; (*rate at which infected become contagious*)
y = 0.2 ;(*death/remove rate*)
NDSolve[Join[outbreak, initialconditions], {s[t], e[t], i[t], 
  r[t]}, {t, 0, 60}]

But it gives the error There are fewer dependent variables, {e[t],i[t],r[t],s[t]}, than equations, so the system is overdetermined. >>

$\endgroup$
  • $\begingroup$ in your initial conditions, you should have == instead of =. Make sure you quit your Mathematica session before trying to run it again. $\endgroup$ – Eric Brown Nov 7 '15 at 22:44
  • $\begingroup$ @EricBrown I end up with this NDSolve[Join[outbreak, initialconditions], {s[t], e[t], i[t], r[t]}, {t, 0, 60}] and it says the error comes up that there are fewer dependant variables than equations so the system is overdetermined $\endgroup$ – xCanaan Nov 7 '15 at 22:53
  • $\begingroup$ I think you're closer, not sure if your Initial Value Problem (IVP), e.g. s[0]==p needs the derivatives specified too. $\endgroup$ – Eric Brown Nov 7 '15 at 22:55
  • $\begingroup$ I posted my new code that seems to be closer, but still throws up the fewer dependent variables so the system is over determined. $\endgroup$ – xCanaan Nov 7 '15 at 23:00
7
$\begingroup$

I think you're very close. I reformatted the code a bit, and removed the conservation equation for p[t].

p = 10 (*total population*)
b = 0.9 (*infection rate*)
d = 0.5 (*rate at which infected become contagious*)
y = 0.2 (*death/remove rate*)

outbreak = {
  s'[t] == -b*s[t]*i[t]/p,
  e'[t] == b*s[t]*i[t]/p - d*e[t],
  i'[t] == d*e[t] - y*i[t],
  r'[t] == y*i[t]
 }
initialconditions = {
  s[0] == p,
  e[0] == 1,
  i[0] == 0,
  r[0] == 0
}

soln = Last@NDSolve[
   {
    outbreak,
    initialconditions
   },
   {s[t], e[t], i[t], r[t]}
   , {t, 0, 60}]

soln

Plot[
 s[t] /. soln
 , {t, 0, 60}, PlotRange -> All, PlotLabel -> "s(t)"]

s(t)

I would have to check this carefully, but I think that Mathematica's ODE integrator is capable of maintaining constant "density" without resorting to an additional variable for p[t]. Specifying this additional equation will lead to an overdetermined system, unless you can put it in the variable list to be solved.

Also, for IVP problems of this type, only s[0] or s'[0] need to be specified, and for SEIR I suspect that we know s[0] much more often than the derivatives. (Furthermore, at $t=0$, the derivatives are not 0! It will change straight away.)

Good luck! Mathematica is great for solving ODEs found in disease modeling.

$\endgroup$
  • 1
    $\begingroup$ Oh this is just perfect! Yea, I have a project in which my group has to model an Ebola outbreak, and have to consider various factors to control an outbreak. I just couldn't get the code to work! Thanks! $\endgroup$ – xCanaan Nov 7 '15 at 23:08
  • $\begingroup$ @xCanaan You're welcome. For a lark, you might want to plot p[t] (which you can get by summing the solutions $s(t) + e(t) + i(t) + r(t)$ just to make sure that the quantity is conserved. It should (numerically) be a horizontal line. Significant deviations would indicate that you need to change methods, e.g. use Method -> "BDF" For some particular parameters and infectious disease models, the equations are remarkably stiff. $\endgroup$ – Eric Brown Nov 8 '15 at 1:48
  • $\begingroup$ I seem to get p[t]=11 using this method. (It's constant up to machine precision.) This is caused by the ICs s[0]==p, e[0]==1. Taking s[0]==p-1 works. $\endgroup$ – Eric Towers Nov 8 '15 at 6:07

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.