# Trying to model a system of differential equations for outbreaks

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

• in your initial conditions, you should have == instead of =. Make sure you quit your Mathematica session before trying to run it again. – Eric Brown Nov 7 '15 at 22:44
• @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 – xCanaan Nov 7 '15 at 22:53
• I think you're closer, not sure if your Initial Value Problem (IVP), e.g. s[0]==p needs the derivatives specified too. – Eric Brown Nov 7 '15 at 22:55
• I posted my new code that seems to be closer, but still throws up the fewer dependent variables so the system is over determined. – xCanaan Nov 7 '15 at 23:00

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}]


Plot[
s[t] /. soln
, {t, 0, 60}, PlotRange -> All, PlotLabel -> "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.

• 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! – xCanaan Nov 7 '15 at 23:08
• @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. – Eric Brown Nov 8 '15 at 1:48
• 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. – Eric Towers Nov 8 '15 at 6:07