Trying to model a system of differential equations for outbreaks

Question

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

Answer 1

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.