Disease Network Model Using NDSolve

I have a basic disease model using NDSolve:

b0 = .4;
b1 = .02;
q = .15;
k = .61;
rr = .135;
dr = .82;
tau = 75;
init = {0.99987, 0, 0.00013, 0, 0};

seir[B0_, B1_, Q_, K_, γ_, f_, τ_] := {s[t], e[t], i[t], r[t], d[t], β[t]} /.
First@NDSolve[
{s'[t] == -β[t]*s[t]*i[t],
e'[t] == β[t]*s[t]*i[t] - (K)*e[t],
i'[t] == (K)*e[t] - γ*i[t],
r'[t] == γ *(1 - f)*i[t],
d'[t] == γ *f*i[t],

s[0] == init[[1]], e[0] == init[[2]], i[0] == init[[3]],
r[0] == init[[4]], d[0] == init[[5]],

β[t] == B0 + ((B1 - B0)/(1 + Exp[-Q*(t - τ)]))},
{s[t], e[t], i[t], r[t], d[t], β[t],

}, {t, 0, 365}]

Plot[Evaluate[seir[b0, b1, q, k, rr, dr, tau]], {t, 0, 160},
PlotLegends -> {"Susceptible", "Exposed", "Infected", "Removed",
"Dead", "Contact Rate"}, Frame -> True, PlotRange -> {0, 1},
ImageSize -> {500, 500}]


How can I apply idea this to a network of n nodes to account for exposed individuals moving from one place to another. For example:

Where P(t) is some function I'm using based on distance between nodes.

Here is an attempt:

b0 = .3;
b1 = .02;
q = .15;
k = 1.6;
rr = .135;
dr = .82;
tau = 75;
w = .5;
init = {0.99987, 0, 0.00013, 0, 0, 0.99987, 0, 0.00013, 0, 0, 0.999,
0, 0.001, 0, 0, 0.9993, 0, 0.0007, 0, 0};
coords = {{40.6642738, -73.9385004}, {34.0193936, -118.4108248}, \
{41.8375511, -87.6818441}, {29.7804724, -95.3863425}};
P[u_, v_] := 1/EuclideanDistance[coords[[u]], coords[[v]]];

seir[B0_, B1_, Q_, K_, \[Gamma]_, f_, \[Tau]_, \[Omega]_] := {
s1[t], e1[t], i1[t], r1[t], d1[t],
s2[t], e2[t], i2[t], r2[t], d2[t],
s3[t], e3[t], i3[t], r3[t], d3[t],
s4[t], e4[t], i4[t], r4[t], d4[t],
\[Beta][t]} /. First@NDSolve[{
s1'[t] == -\[Beta][t]*s1[t]*i1[t],
e1'[t] == \[Beta][t]*s1[t]*i1[t] - (K)*
e1[t] + (e2[t]*\[Omega]*P[1, 2] + e3[t]*\[Omega]*P[1, 3] +
e4[t]*\[Omega]*P[1, 4]) - (e1[t]*P[1, 2]*\[Omega] +
e1[t]*P[1, 3]*\[Omega] + e1[t]*P[1, 4]*\[Omega]),
i1'[t] == (K)*e1[t] - \[Gamma]*i1[t],
r1'[t] == \[Gamma]*(1 - f)*i1[t],
d1'[t] == \[Gamma]*f*i1[t],

s2'[t] == -\[Beta][t]*s2[t]*i2[t],
e2'[t] == \[Beta][t]*s2[t]*i2[t] - (K)*
e2[t] + (e1[t]*\[Omega]*P[2, 1] + e3[t]*\[Omega]*P[2, 3] +
e4[t]*\[Omega]*P[2, 4]) - (e2[t]*P[2, 1]*\[Omega] +
e2[t]*P[2, 3]*\[Omega] + e2[t]*P[2, 4]*\[Omega]),
i2'[t] == (K)*e2[t] - \[Gamma]*i2[t],
r2'[t] == \[Gamma]*(1 - f)*i2[t],
d2'[t] == \[Gamma]*f*i2[t],

s3'[t] == -\[Beta][t]*s3[t]*i3[t],
e3'[t] == \[Beta][t]*s3[t]*i3[t] - (K)*
e3[t] + (e1[t]*\[Omega]*P[3, 1] + e2[t]*\[Omega]*P[3, 2] +
e4[t]*\[Omega]*P[3, 4]) - (e3[t]*P[3, 1]*\[Omega] +
e3[t]*P[3, 2]*\[Omega] + e3[t]*P[3, 4]*\[Omega]),
i3'[t] == (K)*e3[t] - \[Gamma]*i3[t],
r3'[t] == \[Gamma]*(1 - f)*i3[t],
d3'[t] == \[Gamma]*f*i3[t],

s4'[t] == -\[Beta][t]*s4[t]*i4[t],
e4'[t] == \[Beta][t]*s4[t]*i4[t] - (K)*
e4[t] + (e1[t]*\[Omega]*P[4, 1] + e2[t]*\[Omega]*P[4, 2] +
e3[t]*\[Omega]*P[4, 3]) - (e4[t]*P[4, 1]*\[Omega] +
e4[t]*P[4, 2]*\[Omega] + e4[t]*P[4, 3]*\[Omega]),
i4'[t] == (K)*e4[t] - \[Gamma]*i4[t],
r4'[t] == \[Gamma]*(1 - f)*i4[t],
d4'[t] == \[Gamma]*f*i4[t],

s1[0] == init[[1]], e1[0] == init[[2]], i1[0] == init[[3]],
r1[0] == init[[4]], d1[0] == init[[5]],
s2[0] == init[[6]], e2[0] == init[[7]], i2[0] == init[[8]],
r2[0] == init[[9]], d2[0] == init[[10]],
s3[0] == init[[11]], e3[0] == init[[12]], i3[0] == init[[13]],
r3[0] == init[[14]], d3[0] == init[[15]],
s4[0] == init[[16]], e4[0] == init[[17]], i4[0] == init[[18]],
r4[0] == init[[19]], d4[0] == init[[20]],

\[Beta][t] == B0 + ((B1 - B0)/(1 + Exp[-Q*(t - \[Tau])]))}, {

s1[t], e1[t], i1[t], r1[t], d1[t],
s2[t], e2[t], i2[t], r2[t], d2[t],
s3[t], e3[t], i3[t], r3[t], d3[t],
s4[t], e4[t], i4[t], r4[t], d4[t],

\[Beta][t]}, {t, 0, 365}]


As you can see this gets ugly fast. How can I simplify?

• I shelved this project for a while, but am working on it again. I appreciate any suggestions. Jun 13, 2015 at 17:11