0
$\begingroup$

$$ \dot S (t)=\Lambda-\beta_1 S(t)A(t)-\beta_2S(t)I(t)-\mu S(t)$$ $$ \dot A (t)= \beta_1 S(t)A(t)+\beta_2S(t)I(t)-(\mu+\alpha) A(t)$$ $$ \dot I (t)=\alpha A(t)-(\mu+\theta)I(t)$$ $$ \dot Q (t)=\theta I(t)-(\mu+\mu_1+\gamma)Q(t)$$ $$\dot R (t)=\gamma Q(t)-\mu R(t)$$ On the above problem, where all parameters$$\Lambda, \beta_1,\beta_2, \mu, \mu_1, \alpha, \theta, \gamma$$ are known, I wanted to calculate the Basic reproduction number of the system of ordinary differential equations using next-generation matrix method(or even any other known method). However, I could not find any inbuilt function for this method. Can somebody point me in the right direction? Or can somebody help me to write MATLAB code for it?

thanks for any help.

$\endgroup$

1 Answer 1

2
$\begingroup$

The following derivation follows closely Vivas-Barber, A. et al. Dynamics of an "SAIQR" influenza model. 2014.
I have used Wolfram Mathematica for calculations as shown below. You should validate yourself that the calculations are correct.

First, we define the state vector $E=(S,A,I,Q,R)$ and the total population size $$N=S+A+I+Q+R.$$ Solving $\dot{E}=0$ for the trivial, disease-free equilibrium, we obtain: $$N = \frac{\Lambda -\mu_1 Q}{\mu}.$$ Next, we calculate the Jacobian of the system at $E_0 = \left((\Lambda -\mu_1 Q)/\mu,0,0,0,0\right)$:

$$\left( \begin{array}{ccccc} -\mu & -\frac{\beta_1 \Lambda }{\mu } & -\frac{\beta_2 \Lambda }{\mu } & 0 & 0 \\ 0 & -\alpha +\frac{\beta_1 \Lambda }{\mu }-\mu & \frac{\beta_2 \Lambda }{\mu } & 0 & 0 \\ 0 & \alpha & -\theta -\mu & 0 & 0 \\ 0 & 0 & \theta & -\gamma -\mu -\mu_1 & 0 \\ 0 & 0 & 0 & \gamma & -\mu \\ \end{array} \right)$$ and the corresponding eigenvalues: $$ \begin{array}{l} \lambda_1 = -\mu, \\ \lambda_2 = -\mu, \\ \lambda_3 = -\frac{1}{2 \mu }\left(\sqrt{2 \beta_1 \Lambda \mu (\theta -\alpha )+\mu \left(4 \alpha \beta_2 \Lambda +\mu (\alpha -\theta )^2\right)+\beta_1^2 \Lambda ^2}+\mu (\alpha +\theta +2 \mu )-\beta_1 \Lambda \right), \\ \lambda_4 = \frac{1}{2 \mu }\left(\sqrt{2 \beta_1 \Lambda \mu (\theta -\alpha )+\mu \left(4 \alpha \beta_2 \Lambda +\mu (\alpha -\theta )^2\right)+\beta_1^2 \Lambda ^2}-\mu (\alpha +\theta +2 \mu )+\beta_1 \Lambda \right), \\ \lambda_5 = -\gamma -\mu -\mu_1. \end{array}$$ Rephrasing the conditions $\lambda_3<0$ and $\lambda_4<0$ as $\mathfrak{R}_0 < 1$ leads to: $$\mathfrak{R}_0 = \frac{\alpha \beta_2 \Lambda }{(\theta +\mu ) (\mu (\alpha +\mu )-\beta_1 \Lambda )},$$ and a necessary condition on the parameters: $\beta_1 < (\alpha +\theta +2 \mu )\mu/\Lambda $.


Clear["Global`*"];

eqs = {
   S'[t] == Λ - β1 S[t] A[t] - β2 S[t] Ι[t] - μ S[t],
   A'[t] == β1 S[t] A[t] + β2 S[t] Ι[t] - (μ + α) A[t],
   Ι'[t] == α A[t] - (μ + θ) Ι[t],
   Q'[t] == θ Ι[t] - (μ + μ1 + γ) Q[t],
   R'[t] == γ Q[t] - μ R[t]
   };
vars = DeleteDuplicates@Cases[eqs, _Symbol[t], All];
params = DeleteDuplicates@Cases[eqs, _Symbol, All];

$Assumptions = Thread[params > 0];

Eliminate[{Total[Last /@ eqs] == 0, Ν == Total[vars]}, {R[t]}];
Ν0 = First@SolveValues[%, Ν];

J = D[Last /@ eqs, {vars}] /. S[t] -> Ν0 /. (vars -> 0 // Thread);
λ = Eigenvalues[J];

FullSimplify[Reduce[λ[[3]] < 0 && λ[[4]] < 0, Reals]]
(* β1 Λ < μ (α + θ + 2 μ) && β2 < ((θ + μ) (-β1 Λ + μ (α + μ)))/(α Λ) *)
$\endgroup$
1
  • $\begingroup$ Dear @Domen Thank you very much. This advice and your valuable reply helped me a lot and I am now studying more closely and checking what you have said. If there is a problem, I will definitely take your and other researchers' valuable time again. $\endgroup$ Commented Jul 20, 2023 at 11:00

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .