Disease Network Model Using NDSolve

Question

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?