m_{s+2,s}I_{s}^\diamond & \dots & \dots & m_{n,s}I_{s}^\diamond\\
 m_{s,s+1}I_{s+1}^\diamond   & -\tilde{b}_{s+1,s+1} &  m_{s+2,s+1}I_{s+1}^\diamond & \dots &\dots & m_{n,s+1}I_{s+1}^\diamond\\
  m_{s,s+2}I_{s+2}^\diamond &   m_{s+1,s+2}I_{s+2}^\diamond & -\tilde{b}_{33} & \dots& \dots&    m_{n,s+2}I_{s+2}^\diamond\\
 \vdots& \vdots & \vdots& \ddots & \ddots& \vdots \\
m_{s,n}I_{n}^\diamond & m_{s+1,n}I_{n}^\diamond & m_{s+2,n}I_{n}^\diamond  & \hdots&\hdots &-\tilde{b}_{nn}
\end{pmatrix},$$ where $\tilde{b}_{kk}=\sum_{j=s}^nm_{kj}I_j^\diamond$ for $s\leq k\leq n$. Since $M_{22}$ is irreducible and $I^\diamond_i>0$ for all $s\leq i\leq n$, the matrix $\tilde B$ is irreducible. Moreover, $\tilde B$ is the Laplacian matrix of the graph interconnecting the stages $I_i$ for $s\leq i\leq n$. Hence, as previously stated, Kirchhoff’s matrix tree theorem affirms that the solution of $\tilde B\textbf{c}_2=0$ is such that $c_i=-C_{ii}\gg0$, where $C_{ii}$ is the cofactor of $i^{th}$ diagonal element of $\tilde B$. Hence $\textbf{c}_2\gg0$.

Let $\mathbf{I}_h^\diamond=(\mathbf{I}_1^\diamond,\mathbf{I}_2^\diamond)^T$ and consider the Lyapunov function candidate $ \mathcal V= \mathcal V_h+ \mathcal V_v$, where $$\mathcal V_h=c_1 \frac{\nu_h}{\alpha_h} E_h+\bigl\langle\textbf{c}_1   \, \mbox{\Large $\mid$} \, \mathbf{I}_1  \bigr\rangle
  +\sum_{i=s}^nc_{i}\int_{I_{i}^\diamond}^{I_{i}}\left(1-\frac{I_{i}^\diamond}{x}\right)dx,\quad\textrm{and}\quad  \mathcal V_v=c_vE_v+c_v\frac{\alpha_v}{\nu_v}I_v,$$ where $c_v=c_1\frac{\nu_h}{\alpha_h}\frac{a\beta_{vh}\Lambda_h}{\mu_hN_h}\frac{\nu_v}{(\mu_v+\delta_v)\alpha_v}$, $\textbf{c}_1$ is positive vector to be determined later. The derivative of $\mathcal V$ along the trajectories of System (\[ModelAD2Compact\]) is: $$\begin{aligned}
\label{LyapBoundary}
  \dot {\mathcal V_h} &=&c_1 \frac{\nu_h}{\alpha_h}\dot{E}_h+\bigl\langle\textbf{c}_1   \, \mbox{\Large $\mid$} \, \dot{\mathbf{I}}_1^\diamond  \bigr\rangle
  +\sum_{i=s}^nc_i\left(1-\frac{I_i^\diamond}{I_i}\right)\dot I_i+c_v\dot E_v+c_v\frac{\alpha_v}{\nu_v}\dot I_v\nonumber\\
  &=&c_1 a\,\beta_{vh}\frac{\nu_h}{\alpha_h}\dfrac{S_{h}I_{v}}{N_{h}}+\bigl\langle\textbf{c}_1   \, \mbox{\Large $\mid$} \, (-\textrm{diag}(\tilde\alpha)+M_{11})\mathbf{I}_1  \bigr\rangle+\sum_{i=s}^nc_i\left(1-\frac{I_i^\diamond}{I_i}\right)\dot I_i,    \end{aligned}$$ where $\tilde\alpha=(\alpha_1,\dots,\alpha_{s-1})$. Moreover, as 