1
$\begingroup$

Our system:

\begin{align} \dot S &= \mu -\beta_1 S I_2 -\beta_2 S J-\beta_3 S A - \nu S\\[2ex] \dot {I_1} &= p \beta_1 S I_1 +q \beta_2 S J +r \beta_3 S A + \xi_1 J- b_1 I_1 \\[2ex] \dot{I_2} &=(1 - p) \beta_1 S I_2 + (1 - q) \beta_2 S J + (1 - r) \beta_3 S A + \epsilon I_1+ \xi_2 J - b_2 I_2 \\[2ex] \dot{J} &=p_1 I_2 -b_3 J \\[2ex] \dot{A} &= p_2 J -b_4 A \end{align}

With $N= S+I_1+I_2 +J+A$. Implying $$\dot N = \mu - \nu N -\alpha A \leq \mu-\nu N $$

I tried:

tmax = 1000;
\[Beta]1 = 0.2;
\[Beta]2 = 0.1;
\[Beta]3 = 0.1;
p = 0.8;
q = 0.9;
r = 0.7;
\[Xi]1 = 0.1;
\[Xi]2 = 0.3;
\[Epsilon] = 0.2;
b1 = 0.4;
b2 = 0.3;
b3 = 0.5;
b4 = 0.2;
p1 = 0.3;
p2 = 0.2;
\[Mu] = 0.3;
\[Nu] = 0.1;
SIIJA = NDSolveValue[{
    S'[t] == \[Mu] - \[Beta]1*S[t]*I2[t] - \[Beta]2*S[t]*
       J[t] - \[Beta]3*S[t]*A[t] - \[Nu]*S[t],
    I1'[t] == 
     p*\[Beta]1*S[t]*I2[t] + q*\[Beta]2*S[t]*J[t] + 
      r*\[Beta]3*S[t]*A[t] + \[Xi]1 *J[t] - b1*I1[t],
    I2'[t] == (1 - p) \[Beta]1*S[t]*I2[t] + (1 - q)*\[Beta]2*S[t]*
       J[t] + (1 - r)*\[Beta]3*S[t]*A[t] + \[Epsilon]*I1[t] + \[Xi]2*
       J[t] - b2*I2[t],
    J'[t] == p1*I2[t] - b3*J[t],
    A'[t] == p2*J[t] - b4*A[t] ,
    S[0] == 50,
    I1[0] == 50,
    I2[0] == 50,
    A[0] == 10,
    J[0] == 10},
   {S, I1, I2, J, A},
   {t, 0, tmax}];
{f1, f2, f3, f4, f5} = SIIJA;

st = Style[#, 15, Black] &;

Plot[{f1[t], f2[t], f3[t], f4[t], f5[t]}, {t, 0, tmax}, 
 PlotStyle -> {Blue, Purple, Yellow, Orange, Red}, Frame -> True, 
 FrameLabel -> st /@ {"Time", "Density"}, 
 PlotLegends -> 
  Placed[LineLegend[{Blue, Purple, Yellow, Orange, Red}, {"S(t)", 
     "I1(t)", "I2(t)", "J(t)", "A(t)" }, 
    LegendFunction -> Framed], {0.85, 0.65}], ImageSize -> 500]

but this doesn't seem correct..

EDIT:

Our "full" system:

    tmax = 1000;
\[Beta]1 = 0.00001;
\[Beta]2 = 0.0006;
\[Beta]3 = 0.00001;
p = 0.9;
q = 0.8;
r = 0.3;
\[Xi]1 = 0.8;
\[Xi]2 = 0.9;
\[Epsilon] = 0.0002;
p1 = 0.01;
p2 = 0.03;
\[Alpha] = 0.01;
\[Mu] = 0.55;
\[Nu] = 0.01;
SIIJA = NDSolveValue[{
    S'[t] == \[Mu] - \[Beta]1*S[t]*I2[t] - \[Beta]2*S[t]*
       J[t] - \[Beta]3*S[t]*A[t] - \[Nu]*S[t],
    I1'[t] == 
     p*\[Beta]1*S[t]*I2[t] + q*\[Beta]2*S[t]*J[t] + 
      r*\[Beta]3*S[t]*A[t] + \[Xi]1 *J[t] - (\[Epsilon] + \[Nu])*I1[t],
    I2'[t] == (1 - p) \[Beta]1*S[t]*I2[t] + (1 - q)*\[Beta]2*S[t]*
       J[t] + (1 - r)*\[Beta]3*S[t]*A[t] + \[Epsilon]*I1[t] + \[Xi]2*
       J[t] - (p1 + \[Nu])*I2[t],
    J'[t] == p1*I2[t] - (\[Xi]1 + \[Xi]2 + p2 + \[Nu])*J[t],
    A'[t] == p2*J[t] - (\[Alpha] + \[Nu])*A[t] ,
    S[0] == 50,
    I1[0] == 50,
    I2[0] == 50,
    J[0] == 10,
    A[0] == 10},
   {S, I1, I2, J, A},
   {t, 0, tmax}];
{f1, f2, f3, f4, f5} = SIIJA;

st = Style[#, 15, Black] &;

Plot[{f1[t], f2[t], f3[t], f4[t], f5[t]}, {t, 0, tmax}, 
 PlotStyle -> {Blue, Purple, Yellow, Orange, Red}, Frame -> True, 
 FrameLabel -> st /@ {"Time", "Density"}, 
 PlotLegends -> 
  Placed[LineLegend[{Blue, Purple, Yellow, Orange, Red}, {"S(t)", 
     "I1(t)", "I2(t)", "J(t)", "A(t)" }, 
    LegendFunction -> Framed], {0.85, 0.75}], ImageSize -> 750]

EDIT 2:

Their plots(in the paper attached) are different to the plots I got when I modelled the same system, is there a reason why?

https://reader.elsevier.com/reader/sd/pii/S0307904X12002442?token=7F92FD800A67CDC6EBADE29BD17894CAE44E83139073043B925BCFE37615DC4CD3483D225B3EF9B55F9CE61DCD117F59&originRegion=eu-west-1&originCreation=20211130120104

tmax = 2000;
\[Beta]1 = 0.0001;
\[Beta]2 = 0.006;
p = 0.9;
q = 0.8;
\[Xi]1 = 0.8;
\[Xi]2 = 0.9;
\[Epsilon] = 0.002;
p1 = 0.01;
p2 = 0.03;
\[Alpha] = 0.01;
\[Mu] = 0.01;
\[Nu] = 0.55;
SIIJA = NDSolveValue[{
    S'[t] == \[Nu] - \[Beta]1*S[t]*I2[t] - \[Beta]2*S[t]*J[t] - \[Mu]*
       S[t],
    I1'[t] == 
     p*\[Beta]1*S[t]*I2[t] + 
      q*\[Beta]2*S[t]*J[t] + \[Xi]1 *J[t] - (\[Epsilon] + \[Mu])*I1[t],
    I2'[t] == (1 - p) \[Beta]1*S[t]*I2[t] + (1 - q)*\[Beta]2*S[t]*
       J[t] + \[Epsilon]*I1[t] + \[Xi]2*J[t] - (p1 + \[Mu])*I2[t],
    J'[t] == p1*I2[t] - (\[Xi]1 + \[Xi]2 + p2 + \[Mu])*J[t],
    A'[t] == p2*J[t] - (\[Alpha] + \[Mu])*A[t] ,
    S[0] == 50,
    I1[0] == 50,
    I2[0] == 50,
    J[0] == 20,
    A[0] == 20},
   {S, I1, I2, J, A},
   {t, 0, tmax}];
{f1, f2, f3, f4, f5} = SIIJA;

st = Style[#, 15, Black] &;

Plot[{f1[t], f2[t], f3[t], f4[t], f5[t]}, {t, 0, tmax}, 
 PlotStyle -> {Blue, Purple, Yellow, Orange, Red}, Frame -> True, 
 FrameLabel -> st /@ {"Time", "Density"}, 
 PlotLegends -> 
  Placed[LineLegend[{Blue, Purple, Yellow, Orange, Red}, {"S(t)", 
     "I1(t)", "I2(t)", "J(t)", "A(t)" }, 
    LegendFunction -> Framed], {0.85, 0.45}], ImageSize -> 500]

Here is the reproduction number:

r0 = ((\[Beta]1 b3 (\[Epsilon] p + 
          b1 (1 - p)) + \[Beta]2 p1 (\[Epsilon] q + 
          b1 (1 - q))) (\[Nu]/\[Mu]))/(b1 b2 b3 - (\[Epsilon] \[Xi]1 +
         b1 \[Xi]2) p1);

EDIT 3

tmax = 2000;
\[Beta]1 = 0.0001;
\[Beta]2 = 0.006;
p = 0.3;
q = 0.4;
\[Xi]1 = 0.001;
\[Xi]2 = 0.003;
\[Epsilon] = 0.0002;
p1 = 0.01;
p2 = 0.03;
\[Alpha] = 0.01;
\[Mu] = 0.01;
\[Nu] = 0.55;
SIIJA = NDSolveValue[{
    S'[t] == \[Nu] - \[Beta]1*S[t]*I2[t] - \[Beta]2*S[t]*J[t] - \[Mu]*
       S[t],
    I1'[t] == 
     p*\[Beta]1*S[t]*I2[t] + 
      q*\[Beta]2*S[t]*J[t] + \[Xi]1 *J[t] - (\[Epsilon] + \[Mu])*I1[t],
    I2'[t] == (1 - p) \[Beta]1*S[t]*I2[t] + (1 - q)*\[Beta]2*S[t]*
       J[t] + \[Epsilon]*I1[t] + \[Xi]2*J[t] - (p1 + \[Mu])*I2[t],
    J'[t] == p1*I2[t] - (\[Xi]1 + \[Xi]2 + p2 + \[Mu])*J[t],
    A'[t] == p2*J[t] - (\[Alpha] + \[Mu])*A[t] ,
    S[0] == 50,
    I1[0] == 50,
    I2[0] == 50,
    J[0] == 20,
    A[0] == 20},
   {S, I1, I2, J, A},
   {t, 0, tmax}];
{f1, f2, f3, f4, f5} = SIIJA;

st = Style[#, 15, Black] &;

Plot[{f1[t], f2[t], f3[t], f4[t], f5[t]}, {t, 0, tmax}, 
 PlotStyle -> {Blue, Purple, Yellow, Orange, Red}, Frame -> True, 
 FrameLabel -> st /@ {"Time", "Density"}, 
 PlotLegends -> 
  Placed[LineLegend[{Blue, Purple, Yellow, Orange, Red}, {"S(t)", 
     "I1(t)", "I2(t)", "J(t)", "A(t)" }, 
    LegendFunction -> Framed], {0.85, 0.75}], ImageSize -> 500]

We know a priori that $N(2000)$ should be approximately $55$. however from the model above if we do f1[2000] + f2[2000] + f3[2000] + f4[2000] + f5[2000] it gives 51.3631. This is what I was trying to explain to Ulrich.

$\endgroup$
31
  • $\begingroup$ Correct the typo, double I2, in SIIJA : {S, I2, I2, J, A} $\endgroup$
    – Akku14
    Nov 29 '21 at 14:55
  • $\begingroup$ @Akku14 Thanks, however $A$ isn't still showing up? Also, I don't know how to include the condition $N=S+I_1+I_2+J+A$ $\endgroup$
    – Math
    Nov 29 '21 at 14:57
  • $\begingroup$ There is a typo in your Plot-command: Change f[5] to f5[t] $\endgroup$ Nov 29 '21 at 14:59
  • 1
    $\begingroup$ @Math S'[0] == -604.5 so it's not surprising that S[t] quickly approaches zero. $\endgroup$
    – Chris K
    Nov 29 '21 at 15:46
  • 1
    $\begingroup$ @UlrichNeumann How do I verify that $S+I_1+I_2+J+A$ equals to $N$ at each point in time? $\endgroup$
    – Math
    Nov 29 '21 at 17:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.