Reproduction number in SAIQR model

Question

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

Answer 1

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 Λ + μ (α + μ)))/(α Λ) *)