I have a system of 6 ODE representing infected subpopulations of a model of epidemic
$$
\begin{cases}
\dfrac{de}{dt} =\dfrac{dy_2}{dt} = \phi y_1(y_4+y_5)-p_2 y_2\\
\dfrac{dq^*}{dt} =\dfrac{dy_3}{dt} = \mu y_2 -P_3 y_3\\
\dfrac{ di_s}{dt} =\dfrac{dy_4}{dt} = \psi y_2 +\alpha y_3 -p_4 y_4\\
\dfrac{di_a}{dt} =\dfrac{dy_5}{dt} = \rho y_2+\sigma y_3 -p_5 y_5 \\
\dfrac{dh}{dt} =\dfrac{dy_6}{dt} = \kappa y_4+l y_5+\nu y_7-p_6 y_6\\
\dfrac{dh_m}{dt} =\dfrac{dy_7}{dt} = m y_5+h y_4-p_7 y_7
\end{cases}
$$
Here $y_2,\dots ,y_7$ are the number of people in different compartments.
Other epidemic parameters along with $y_1$ on the RHS are assumed to be non-negative.
I'm trying to find the basic reproduction number $R_0 $ using a next-generation matrix approach.
The Jacobian matrix ($J$) for the infected sub-population is:
$$
J (y_2, y_3, \ldots, y_7)=
\begin{bmatrix}
-p_2 & 0 & \phi y_1 & \phi y_1 & 0 & 0 \\
\mu & -P_3 & 0 & 0 & 0 & 0 \\
\psi & \alpha & -p_4 & 0 & 0 & 0 \\
\rho & \sigma & 0 & -p_5 & 0 & 0 \\
0 & 0 & \kappa & l &-p_6 & \nu \\
0 & 0 & h & m & 0 & -p_7 \\
\end{bmatrix}
= T + Q ,
$$ where
$$
T=
\begin{bmatrix}
0 & 0 & \phi y_1 &\phi y_1 & 0 & 0\\
0 & 0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 0 & 0
\end{bmatrix}
$$ and
$$
Q=
\begin{bmatrix}
-p_2 & 0 & 0 & 0 & 0 & 0 \\
\mu & -P_3 & 0 & 0 & 0 & 0 \\
\psi & \alpha & -p_4 & 0 & 0 & 0 \\
\rho & \sigma & 0 & -p_5 & 0 & 0 \\
0 & 0 & \kappa & l & -p_6 & \nu \\
0 & 0 & h & m & 0 & -p_7 \\
\end{bmatrix}
$$
For a large population, $R_0$= the dominant eigenvalue of the $NGM=TQ^{-1}$ which can be determined conveniently using computer programming like Matlab.
The inverse matrix of $Q$ is obtained by Gauss Jordan Method:\
$$
\begin{split}
[Q|I] &=
\left(\begin{array}{cccccc|cccccc}
-p_2 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0\\
\mu & -p_3 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\
\psi & \alpha & -p_4 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0\\
\rho & \sigma & 0 & -p_5 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0\\
0 & 0 & \kappa & l & -p_6 & \nu & 0 & 0 & 0 & 0 & 1 & 0\\
0 & 0 & h & m & 0 & -p_7 & 0 & 0 & 0 & 0 & 0 & 1
\end{array}\right)\\
\\
&=
\begin{array}{l}
\left[\begin{array}{cccccc|cccccc}
1 & 0 & 0 & 0 & 0 & 0 & -\frac{1}{p_2 } & 0 & 0 & 0 & 0 & 0\\
0 & 1 & 0 & 0 & 0 & 0 & -\frac{\mu }{p_2 \,p_3 } & -\frac{1}{p_3 } & 0 & 0 & 0 & 0\\
0 & 0 & 1 & 0 & 0 & 0 & -\frac{\sigma_4 }{p_2 \,p_3 \,p_4 } & -\frac{\alpha }{p_3 \,p_4 } & -\frac{1}{p_4 } & 0 & 0 & 0\\
0 & 0 & 0 & 1 & 0 & 0 & -\frac{\sigma_3 }{p_2 \,p_3 \,p_5 } & -\frac{\sigma }{p_3 \,p_5 } & 0 & -\frac{1}{p_5 } & 0 & 0\\
0 & 0 & 0 & 0 & 1 & 0 & -\frac{\sigma_4 \,\sigma_2 }{p_2 \,p_3 \,p_4 \,p_6 \,p_7 }-\frac{\sigma_1 \,\sigma_3 }{p_2 \,p_3 \,p_5 \,p_6 \,p_7 } & -\frac{\alpha \,\sigma_2 }{p_3 \,p_4 \,p_6 \,p_7 }-\frac{\sigma \,\sigma_1 }{p_3 \,p_5 \,p_6 \,p_7 } & -\frac{\sigma_2 }{p_4 \,p_6 \,p_7 } & -\frac{\sigma_1 }{p_5 \,p_6 \,p_7 } & -\frac{1}{p_6 } & -\frac{\nu }{p_6 \,p_7 }\\
0 & 0 & 0 & 0 & 0 & 1 & -\frac{h\,\sigma_4 }{p_2 \,p_3 \,p_4 \,p_7 }-\frac{m\,\sigma_3 }{p_2 \,p_3 \,p_5 \,p_7 } & -\frac{\alpha \,h\,p_5 +m\,p_4 \,\sigma }{p_3 \,p_4 \,p_5 \,p_7 } & -\frac{h}{p_4 \,p_7 } & -\frac{m}{p_5 \,p_7 } & 0 & -\frac{1}{p_7 }
\end{array}\right]
\end{array}\\
\\
&=[I|Q^{-1}]
\end{split}
$$
where
- $\sigma_1 =m\nu +lp_7$
- $\sigma_2 =h\nu +\kappa p_7$,
- $\sigma_3 =\mu \sigma + p_3 \rho$,
- $\sigma_4 =\alpha \mu +p_3 \psi$
Therefore, $$ Q^{-1} = \left(\begin{array}{cccccc} -\frac{1}{p_2 } & 0 & 0 & 0 & 0 & 0\\ -\frac{\mu }{p_2 \,p_3 } & -\frac{1}{p_3 } & 0 & 0 & 0 & 0\\ -\frac{\sigma_4 }{p_2 \,p_3 \,p_4 } & -\frac{\alpha }{p_3 \,p_4 } & -\frac{1}{p_4 } & 0 & 0 & 0\\ -\frac{\sigma_3 }{p_2 \,p_3 \,p_5 } & -\frac{\sigma }{p_3 \,p_5 } & 0 & -\frac{1}{p_5 } & 0 & 0\\ -\frac{\sigma_4 \,\sigma_2 }{p_2 \,p_3 \,p_4 \,p_6 \,p_7 }-\frac{\sigma_1 \,\sigma_3 }{p_2 \,p_3 \,p_5 \,p_6 \,p_7 } & -\frac{\alpha \,\sigma_2 }{p_3 \,p_4 \,p_6 \,p_7 }-\frac{\sigma \,\sigma_1 }{p_3 \,p_5 \,p_6 \,p_7 } & -\frac{\sigma_2 }{p_4 \,p_6 \,p_7 } & -\frac{\sigma_1 }{p_5 \,p_6 \,p_7 } & -\frac{1}{p_6 } & -\frac{\nu }{p_6 \,p_7 }\\ -\frac{h\,\sigma_4 }{p_2 \,p_3 \,p_4 \,p_7 }-\frac{m\,\sigma_3 }{p_2 \,p_3 \,p_5 \,p_7 } & -\frac{\alpha \,h\,p_5 +m\,p_4 \,\sigma }{p_3 \,p_4 \,p_5 \,p_7 } & -\frac{h}{p_4 \,p_7 } & -\frac{m}{p_5 \,p_7 } & 0 & -\frac{1}{p_7 } \end{array}\right) $$ Now, $$ T Q^{-1} = \left(\begin{array}{cccccc} -\frac{\phi \,y_1 \,{\left(\alpha \,\mu \,p_5 +p_3 \,p_5 \,\psi +\mu \,p_4 \,\sigma +p_3 \,p_4 \,\rho \right)}}{p_2 \,p_3 \,p_4 \,p_5 } & -\frac{\phi \,y_1 \,{\left(\alpha \,p_5 +p_4 \,\sigma \right)}}{p_3 \,p_4 \,p_5 } & -\frac{\phi \,y_1 }{p_4 } & -\frac{\phi \,y_1 }{p_5 } & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 \end{array}\right) $$ The eigenvalues of $TQ^{-1}$ are repesctively $$ \left(\begin{array}{c} \lambda_1 \\ \lambda_2 \\ \lambda_3 \\ \lambda_4 \\ \lambda_5 \\ \lambda_6 \end{array}\right)=\left(\begin{array}{c} 0\\ 0\\ 0\\ 0\\ 0\\ -\frac{\phi \,y_1 \,{\left(\alpha \,\mu \,p_5 +p_3 \,p_5 \,\psi +\mu \,p_4 \,\sigma +p_3 \,p_4 \,\rho \right)}}{p_2 \,p_3 \,p_4 \,p_5 } \end{array}\right) $$ while the basic reproduction number is $$ R_0=-\frac{\phi \,y_1 \,{\left(\alpha \,\mu \,p_5 +p_3 \,p_5 \,\psi +\mu \,p_4 \,\sigma +p_3 \,p_4 \,\rho \right)}}{p_2 \,p_3 \,p_4 \,p_5 }. $$ I got the dominant eigenvalue of NGM as a negative term $$ R_0=-\frac{\phi \,y_1 \,{\left(\alpha \,\mu \,p_5 +p_3 \,p_5 \,\psi +\mu \,p_4 \,\sigma +p_3 \,p_4 \,\rho \right)}}{p_2 \,p_3 \,p_4 \,p_5 } $$ which is, I think not acceptable for an epidemic. How can I make this correct?