## Appendix B Global Stability Result Proof  
The global stability proof of the endemic point \(E_{f}\) (Proposition 3) is given as follows. Consider the following Lyapunov function:  
\[\mathcal{L}(x,y,s,v,z)=s+\frac{d_{2}+k_{2}}{k_{2}}y+\frac{d_{3}(d_{2}+k_{2})}{ ak_{2}}v+\frac{p}{c}\frac{d_{2}+k_{2}}{k_{2}}z.\]  
The time derivative is given by:  
\[\dot{\mathcal{L}}(x,y,s,v,z) =\dot{s}+\frac{d_{2}+k_{2}}{k_{2}}\dot{y}+\frac{d_{3}(d_{2}+k_{2}) }{ak_{2}}\dot{v}+\frac{p}{c}\frac{d_{2}+k_{2}}{k_{2}}z.\] \[\dot{\mathcal{L}}(x,y,s,v,z) =\frac{k_{1}xv}{x+v}-(d_{2}+k_{2})s+\frac{d_{2}+k_{2}}{k_{2}}(k_ {2}s-d_{3}y)-\frac{d_{2}+k_{2}}{k_{2}}pyz\] \[+\frac{d_{3}(d_{2}+k_{2})}{ak_{2}}(ay-d_{4}v)+\frac{p}{c}\frac{d_ {2}+k_{2}}{k_{2}}(cyz-bz)\] \[\dot{\mathcal{L}}(x,y,s,v,z) =\frac{k_{1}xv}{x+v}-\frac{d_{3}d_{4}(d_{2}+k_{2})}{ak_{2}}v- \frac{bp}{c}\frac{d_{2}+k_{2}}{k_{2}}z\] \[\dot{\mathcal{L}}(x,y,s,v,z) \leq k_{1}v-\frac{d_{3}d_{4}(d_{2}+k_{2})}{ak_{2}}v\] \[\leq\frac{d_{3}d_{4}(d_{2}+k_{2})}{ak_{2}}\left(R_{0}-1\right)v.\]  
If \(R_{0}<1\), then \(\dot{\mathcal{L}}\leq 0\). Moreover, \(\dot{\mathcal{L}}\leq 0\) when \(v=0\). The largest compact invariant is  
\[E=\{(x,y,s,v,z)\in\mathbb{R}_{*}^{+}|v=0\},\]  
according to [PERSON]'s invariance principle, \(\lim_{+\infty}v(t)=0\), the limit system of equations is  
\[\left\{ \begin{array}{l}\dot{x}=\lambda-d_{1}x,\\ \dot{y}=-d_{2}s-k_{2}s,\\ \dot{s}=k_{2}s-d_{3}y-pyz,\\ \dot{z}=byz-cz.\end{array} \right.\]  
We define another Lyapunov function and for simplicity, using the same notation,  
\[\mathcal{L}(x,s,y,z)=\frac{1}{x_{0}}\left(x-x_{0}-x_{0}\ln\frac{x}{x_{0}} \right)+s+\frac{d_{2}+k_{2}}{k_{2}}y+\frac{p}{c}\frac{d_{2}+k_{2}}{k_{2}}z.\]  
Since \(x_{0}=\frac{\lambda}{d_{1}}\), then  
\[\dot{\mathcal{L}}(x,s,y,z)=d_{1}\left(2-\frac{x}{x_{0}}-\frac{x_{0}}{x} \right)-\frac{d_{3}(d_{2}+k_{2})}{k_{2}}s-\frac{pb}{c}z,\]  
since the arithmetic mean is greater than or equal to the geometric mean, it follows that  
\[2-\frac{x}{x_{0}}-\frac{x_{0}}{x}\leq 0;\]  
therefore, \(\dot{\mathcal{L}}\leq 0\) and the equality holds if \(x=x_{0}\) and \(s=y=z=0\), which completes the proof. The global stability proof of the endemic point \(E_{1}\) (Proposition 4) is given as follows. We employ the following Lyapunov function:  
\[\mathcal{L}(x,y,s,v,z) =x-x_{1}-\int_{x_{1}}^{x}\frac{(d_{2}+k_{2})s_{1}}{\frac{k_{1}w_{ 1}}{u+v_{1}}}du+s-s_{1}-s_{1}\ln\frac{s}{s_{1}}\] \[+\frac{d_{2}+k_{2}}{k_{2}}(y-y_{1}-y_{1}\ln\frac{y}{y_{1}})+\frac {d_{3}(d_{2}+k_{2})}{ak_{2}}(v-v_{1}-v_{1}\ln\frac{v}{v_{1}})\] \[+\frac{p}{c}\frac{d_{2}+k_{2}}{k_{2}}z.\]  
We have  
\[\mathcal{L}(x,y,s,v,z) =\dot{x}-(d_{2}+k_{2})s_{1}\frac{x+v_{1}}{k_{1}xv_{1}}\dot{x}+\dot {s}-\frac{s_{1}}{s}\dot{s}+\frac{d_{2}+k_{2}}{k_{2}}(\dot{y}-\frac{y_{1}}{y} \dot{y})\] \[+\frac{d_{3}(d_{2}+k_{2})}{ak_{2}}(\dot{v}-\frac{v_{1}}{v}\dot{v })+\frac{p}{c}\frac{d_{2}+k_{2}}{k_{2}}\dot{z}.\]  
On the other hand, we have  
\[\left\{ \begin{array}{l}\lambda=d_{1}x_{1}+(d_{2}+k_{2})s_{1},\\ \frac{k_{1}x_{1}v_{1}}{x_{1}+v_{1}}=(d_{2}+k_{2})s_{1},\\ \frac{s_{1}}{v_{1}}=\frac{d_{3}d_{4}}{ak_{2}},\,\,\,\frac{y_{1}}{v_{1}}=\frac {d_{4}}{a},\,\,\,\frac{s_{1}}{y_{1}}=\frac{d_{3}}{k_{2}}.\end{array}