Impact of human mobility on the transmission dynamics of infectious diseases

Abstract

Spatial heterogeneity is an important aspect to be studied in infectious disease models. It takes two forms: one is local, namely diffusion in space, and other is related to travel. With the advancement of transportation system, it is possible for diseases to move from one place to an entirely separate place very quickly. In a developing country like India, the mass movement of large numbers of individuals creates the possibility of spread of common infectious diseases. This has led to the study of infectious disease model to describe the infection during transport. An SIRS-type epidemic model is formulated to illustrate the dynamics of such infectious disease propagation between two cities due to population dispersal. The most important threshold parameter, namely the basic reproduction number, is derived, and the possibility of existence of backward bifurcation is examined, as the existence of backward bifurcation is very unsettling for disease control and it is vital to know from modeling analysis when it can occur. It is shown that dispersal of populations would make the disease control difficult in comparison with nondispersal case. Optimal vaccination and treatment controls are determined. Further to find the best cost-effective strategy, cost-effectiveness analysis is also performed. Though it is not a case study, simulation work suggests that the proposed model can also be used in studying the SARS epidemic in Hong Kong, 2003.

Appendices

Appendix 1

For local stability analysis, we use the Routh–Hurwitz criteria and we consider the Jacobian matrix of (3) at $E_{2}^{0}(S^{2},I^{2},R^{2})$

$$\begin{aligned} J=\left( \begin{array}{lll} a_{11} &{} \quad a_{12} &{} \quad a_{13} \\ a_{21} &{} \quad a_{22} &{} \quad 0 \\ 0 &{} \quad a_{32} &{} \quad a_{33} \\ \end{array} \right) . \end{aligned}$$

where $a_{11}=-\beta I^{2}-d, a_{12}=-\beta S^{2}, a_{13}=p, a_{21}=\beta I^{2}, a_{22}=\beta S^{2}-(d+\delta +m)-\frac{r u_{2}}{(1+b u_{2} I^{2})^{2}}, a_{23}=0, a_{31}=0, a_{32}=m+\frac{r u_{2}}{(1+b u_{2} I^{2})^{2}}, a_{33}=-(p+d)$.

Now let the characteristic equation of the Jacobian matrix J be $\lambda ^3+c_1\lambda ^2+c_2\lambda +c_3=0$, where

$$\begin{aligned}c_1&=-(a_{11}+a_{22}+a_{33}),\\ c_2&=a_{11}a_{22}+a_{22}a_{33}+a_{33}a_{11}-a_{12}a_{21}\,{{\rm and}}\\ c_3&=a_{12}a_{21}a_{33}-a_{13}a_{21}a_{32}-a_{11}a_{22}a_{33}.\end{aligned}$$

Then, as stated in the Routh–Hurwitz criteria, the eigenvalues of J have negative real parts if $c_{1}$, $c_{3}$ and $c_1c_2-c_3$ all are positive.

Now we calculate the coefficients as

$$\begin{aligned} c_{1}& {}= p+\delta +m+3d+\beta I^2+\frac{ru_2}{(1+bu_2I^2)^2}-\beta S^2,\\ c_2& {}= (\beta I^2+d)(p+d)+\frac{ru_2 \beta I^{2}}{(1+bu_2I^2)^2}\\&\quad +\, (d+\delta +m)(\beta I^2+2d+p)\\&\quad +\,(2d+p) \left( \frac{ru_2}{(1+bu_2I^2)^2} -\beta S^2 \right) \hbox { and}\\ c_3& {}= (p+d)\left( \beta ^{2} S^2 I^2 +p\beta I^2+(\beta I^2+d)(d+\delta +m)\right) \\&\quad +\,(\beta I^2+d) (p+d) \left( \frac{ru_2}{(1+bu_2I^2)^2}-\beta S^2 \right) . \end{aligned}$$

Also,

$$\begin{aligned} c_1c_2-c_3& {}= \left( \frac{ru_2}{(1+bu_2I^2)^2}-\beta S^2 \right) ^{2} (\beta I^{2}+p+2d)\\&\quad +\,\left( \frac{ru_2}{(1+bu_2I^2)^2}-\beta S^2 \right) \left( \beta ^{2} S^{2} I^{2}+(\beta I^{2}+\right. p\\&\quad \left. +\,2d)(p+2\delta +4d+2m+\beta I^{2})\right) \\&\quad +\,(p+\delta +m+3d+\beta I^{2})(d+\delta +m)\\&\quad (\beta I^{2}+p+2d)+\beta ^{2} I^{2} S^{2}\\&\quad (\delta +m+2d+\beta I^{2})+(p+d)((\beta I^{2}+d)\\&\quad (2d+\beta I^{2})+pd) \end{aligned}$$

Now it is easy to note that for $\frac{ru_2}{(1+bu_2I^2)^2}>\beta S^2$, $c_{1},c_{2},c_{3}$ all are positive and also $c_1c_2>c_3$. Thus, for $\frac{ru_2}{(1+bu_2I^2)^2}>\beta S^2$, Routh–Hurwitz criteria are satisfied. Hence, the endemic steady state $E_{2}^{0}$ is locally asymptotically stable for $R^{0}_{0}>1$ and $\frac{ru_2}{(1+bu_2I^2)^2}>\beta S^2$.

Appendix 2

The Jacobian matrix of system (5) at $E_2^1$ is given by

$$\begin{aligned} J(E_2^1)=\left( \begin{array}{ll} A_4 &{} \quad B \\ B &{} \quad A_4 \\ \end{array} \right) \end{aligned}$$

where

$$\begin{aligned} A_4=\left( \begin{array}{lll} -(\beta I_1^{1*}+d+\alpha ) &{} \quad -\,\beta S_1^{1*} &{} \quad p \\ \beta I_1^{1*} &{} \quad \beta S_1^{1*} -(d+\delta +m)-\frac{ru_2}{(1+bu_2I_1^{1*})^2} &{} \quad 0 \\ 0 &{} \quad m+\frac{ru_2}{(1+bu_2I_1^{1*})^2} &{} \quad -(d+p) \\ \end{array} \right) \end{aligned}$$

and B is the same as earlier. We observe that $A_4+B$ is similar as $J(E_2^0)$ and the detailed proof is similar as given in Appendix 1. Hence, $A_4+B$ is stable if $R^{1}_{0}>1$ and $\frac{ru_2}{(1+bu_2I_1^{1*})^2}>\beta S_1^{1*}$.

Now

$$\begin{aligned}A_4-B=\left( \begin{array}{lll} -(\beta I_1^{1*}+d+2\alpha ) &{} \quad -\beta S_1^{1*} &{} \quad p \\ \beta I_1^{1*} &{} \quad \beta S_1^{1*} -(d+\delta +m)-\frac{ru_2}{(1+bu_2I_1^{1*})^2} &{} \quad 0 \\ 0 &{} \quad m+\frac{ru_2}{(1+bu_2I_1^{1*})^2} &{} \quad -(d+p) \\ \end{array} \right) .\end{aligned}$$

Thus, it is enough to verify that the matrix $A_4-B$ fulfills the Routh–Hurwitz criteria. We already checked that $R_0^0=R_0^1$. Now let the characteristic equation of the matrix $A_4-B$ be $\lambda ^3+d_1\lambda ^2+d_2\lambda +d_3=0$, then

$$\begin{aligned} d_1& {}= 3d+2 \alpha +p+\delta +m+\beta I_1^{1*}+\left( \frac{ru_2}{(1+bu_2I_1^{1*})^2}-\beta S_1^{1*}\right) ,\\ d_2& {}= (\beta I_1^{1*}+d+2\alpha )(p+d)+(\beta I_1^{1*}+2d+2\alpha +p ) (d+\delta +m)\\&\quad +\,\frac{ru_2 \beta I_1^{1*}}{(1+bu_2I_1^{1*})^2}+(p+2d+2\alpha )\\&\quad \,\left( \frac{ru_2}{(1+bu_2I_1^{1*})^2}-\beta S_1^{1*}\right) \hbox { and}\\ d_3& {}= (p+d)\left( \beta ^{2} S_1^{1*} I_1^{1*} +p\beta I_1^{1*} +(\beta I_1^{1*}+d+2\alpha )(d+\delta +m)\right) \\&\quad +\,(\beta I_1^{1*}+d+2\alpha ) (p+d) \left( \frac{ru_2}{(1+bu_2I_1^{1*})^2}-\beta S_1^{1*} \right) . \end{aligned}$$

Also after some simplifications, we obtain

$$\begin{aligned} d_1d_2-d_3& {}= \left( \frac{ru_2}{(1+bu_2I_1^{1*})^2}-\beta S_1^{1*} \right) ^{2} (\beta I_1^{1*}+p+2d+2\alpha )\\&\quad +\,\left( \frac{ru_2}{(1+bu_2I_1^{1*})^2}-\beta S_1^{1*} \right) \left( \beta ^{2} S_1^{1*}I_1^{1*}\right. \\&\quad +\,(\beta I_1^{1*}+p+2d+2\alpha )\\&\quad \left. (p+2\delta +4d+2m+2\alpha +\beta I_1^{1*})\right) \\&\quad +\,(p+\delta +m+3d+\beta I_1^{1*}+2\alpha )(d+\delta +m)\\&\quad \,(\beta I_1^{1*}+p+2d+2\alpha )\\&\quad +\,\beta ^{2} S_1^{1*}I_1^{1*}(\delta +m+2d+\beta I_1^{1*}+2\alpha )\\&\quad +\,(p+d)\left( (\beta I_1^{1*}+d+2\alpha )(2d+\beta I_1^{1*}+2\alpha )\right. \\&\left. \quad +\,p(d+2\alpha )\right) \end{aligned}$$

Now it is easy to observe that for $\frac{ru_2}{(1+bu_2I_1^{1*})^2}>\beta S_1^{1*}$, $d_{i}>0, i=1,2,3$ and $d_1d_2>d_3$. Then, all the conditions of Routh–Hurwitz criteria are satisfied for $\frac{ru_2}{(1+bu_2I_1^{1*})^2}>\beta S_1^{1*}$. Hence, $A_4-B$ is stable for $R_0^1>1$, $\frac{ru_2}{(1+bu_2I_1^{1*})^2}>\beta S_1^{1*}$.

Therefore, combining the above two cases, we conclude that the endemic steady state $E_2^1$ is locally asymptotically stable if $R_{0}^{1}>1$ and $\frac{ru_2}{(1+bu_2I_1^{1*})^2}>\beta S_1^{1*}$. Hence, the theorem is proved.

Appendix 3

The Jacobian matrix of system (1) at $E_2^2$ is given by

$$\begin{aligned} J(E_2^2)=\left( \begin{array}{ll} P_3 &{} \quad Q_3 \\ Q_3 &{} \quad P_3 \\ \end{array} \right) , \end{aligned}$$

where

$$\begin{aligned} P_3=\left( \begin{array}{lll} -(d+\beta I_1^{2*}+\alpha ) &{} \quad -(\beta S_1^{2*}) &{} \quad p \\ \beta I_1^{2*} &{} \quad \beta S_1^{2*}-(d+\delta +m+\alpha +\frac{ru_2}{(1+bu_2I_1^{2*})^2}) &{} \quad 0 \\ 0 &{} \quad m+\frac{ru_2}{(1+bu_2I_1^{2*})^2} &{} \quad -(p+d+\alpha ) \\ \end{array} \right) \end{aligned}$$

and

$$\begin{aligned} Q_3=\left( \begin{array}{lll} \alpha -\gamma \alpha I_1^{2*} &{} \quad -\gamma \alpha S_1^{2*} &{} \quad 0 \\ \gamma \alpha I_1^{2*} &{} \quad \alpha +\gamma \alpha S_1^{2*} &{} \quad 0 \\ 0 &{} \quad 0 &{} \quad \alpha \\ \end{array} \right) . \end{aligned}$$

Now to check that the matrix $P_3+Q_3$ satisfies the Routh–Hurwitz criteria, we consider the characteristic equation of the matrix as $\lambda ^3+m_1\lambda ^2+m_2\lambda +m_3=0$, where

$$\begin{aligned} m_{1}& {}= p+\delta +m+3d+\left( \beta +\gamma \alpha \right) I_1^{2*}\\&\quad +\,\frac{ru_2}{(1+bu_2I_1^{2*})^2}-\left( \beta +\gamma \alpha \right) S_1^{2*},\\ m_2& {}= ((\beta +\gamma \alpha )I_1^{2*}+d)(p+d)+\frac{ru_2 (\beta +\gamma \alpha )I_1^{2*}}{(1+bu_2I_1^{2*})^2}\\&\quad +\,(d+\delta +m)((\beta +\gamma \alpha ) I_1^{2*}+2d+p)\\&\quad +\,(2d+p) \left( \frac{ru_2}{(1+bu_2I_1^{2*})^2}- (\beta +\gamma \alpha )S_1^{2*} \right) \hbox { and }\\ m_3& {}= (p+d)\left( (\beta +\gamma \alpha )^{2} S_1^{2*}I_1^{2*} +p(\beta +\gamma \alpha ) I_1^{2*} \right. \\&\quad \left. +\,(d+\delta +m) (d+(\beta +\gamma \alpha ) I_1^{2*})\right) \\&\quad +\,(p+d) (d+(\beta +\gamma \alpha ) I_1^{2*}) \left( \frac{ru_2}{(1+bu_2I_1^{2*})^2}- (\beta +\gamma \alpha )S_1^{2*} \right) . \end{aligned}$$

Also after some simplifications, we obtain

$$\begin{aligned}&m_1m_2-m_3\\&\quad =\left( \frac{ru_2}{(1+bu_2I_1^{2*})^2}-\beta S_1^{2*} \right) ^{2} ((\beta +\gamma \alpha ) I_1^{2*}+p+2d)\\&\qquad +\,\left( \frac{ru_2}{(1+bu_2I_1^{2*})^2}-\beta S_1^{2*} \right) \left( (\beta +\gamma \alpha )^{2} S_1^{2*}I_1^{2*}\right. \\&\qquad \left. +\,((\beta +\gamma \alpha ) I_1^{2*}+p+2d)(p+2\delta +4d+2m+(\beta +\gamma \alpha ) I_1^{2*})\right) \\&\qquad +\,(p+\delta +m+3d+(\beta +\gamma \alpha ) I_1^{2*})(d+\delta +m)\\&\qquad ((\beta +\gamma \alpha ) I_1^{2*}+p+2d)\\&\qquad +\,(\beta +\gamma \alpha )^{2} S_1^{2*} I_1^{2*}(\delta +m+2d+(\beta +\gamma \alpha ) I_1^{2*})\\&\qquad +\,(p+d)(((\beta +\gamma \alpha ) I_1^{2*}+d)(2d+(\beta +\gamma \alpha ) I_1^{2*})+pd). \end{aligned}$$

Now it is easy to note that if $\frac{ru_2}{(1+bu_2I_1^{2*})^2}>\left( \beta +\gamma \alpha \right) S_1^{2*}$, then $m_{i}>0,i=1,2,3$ and $m_{1}m_{2}>m_{3}$, i.e., the Routh–Hurwitz criteria are satisfied.

Now to check the eigenvalue of the matrix $P_3-Q_3$, we consider that the characteristic equation of the matrix $P_3-Q_3$ is $\lambda ^3+n_1\lambda ^2+n_2\lambda +n_3=0$ where

$$\begin{aligned} n_1& {}= p+\delta +m+3d+6\alpha +(\beta -\gamma \alpha )I_1^{2*}\\&\quad +\,\frac{ru_2}{\left( 1+bu_2I_2^{2*}\right) ^{2}}-(\beta -\gamma \alpha )S_1^{2*},\\ n_2& {}= (p+d+2\alpha ) \left( d+2\alpha +(\beta -\gamma \alpha ) I_1^{2*}\right) \\&\quad +\, (d+\delta +m+2\alpha ) \\& \quad \left( p+2d+4\alpha +(\beta -\gamma \alpha ) I_1^{2*}\right) \\&\quad +\,\frac{ru_2(\beta -\gamma \alpha ) I_1^{2*}}{\left( 1+bu_2I_1^{2*}\right) ^{2} }\\&\quad +\,\left( p+2d+4\alpha \right) \\& \quad \left( \frac{ru_2}{\left( 1+bu_2I_1^{2*}\right) ^{2} }\right. \\&\quad \left. -(\beta -\gamma \alpha ) S_1^{2*} \right) \hbox { and}\\ n_3& {}= (p+d+2\alpha ) \left( (\beta -\gamma \alpha )^{2}I_1^{2*}S_1^{2*}\right. \\&\quad +\,pI_1^{2*} (\beta -\gamma \alpha )+ (d+\delta +m+2\alpha )(d+2\alpha \\&\quad \left. +(\beta -\gamma \alpha )I_1^{2*})\right) \\&\quad +\, (p+d+2\alpha )(d+2\alpha +(\beta -\gamma \alpha )I_1^{2*}) \\& \quad \left( \frac{ru_2}{(1+bu_2I_1^{2*})^2}\right. \\&\quad \left. -(\beta -\gamma \alpha )S_1^{2*} \right) . \end{aligned}$$

Moreover, after some simplifications, we obtain

$$\begin{aligned}&n_1n_2-n_3\\&\quad =\left( \frac{ru_2}{(1+bu_2I_1^{2*})^2}-(\beta -\gamma \alpha ) S_1^{2*} \right) ^{2} \\&\qquad ((\beta -\gamma \alpha ) I_1^{2*}+p+2d+4\alpha )\\&\qquad +\,\left( \frac{ru_2}{(1+bu_2I_1^{2*})^2}-(\beta -\gamma \alpha ) S_1^{2*} \right) \left( (\beta -\gamma \alpha )^{2} S_1^{2*}I_1^{2*}\right. \\&\qquad +\,((\beta -\gamma \alpha ) I_1^{2*}+p+2d+4\alpha )\\&\left. \qquad (p+2\delta +4d+2m+8\alpha +(\beta -\gamma \alpha ) I_1^{2*})\right) \\&\qquad +\,(p+\delta +m+3d+6\alpha +(\beta -\gamma \alpha ) I_1^{2*})(d+\delta +m+2\alpha )\\&\qquad ((\beta -\gamma \alpha ) I_1^{2*}+p+2d+4\alpha )+(\beta -\gamma \alpha )^{2} S_1^{2*} I_1^{2*}\\&\qquad (\delta +m+2d+4\alpha +(\beta -\gamma \alpha ) I_1^{2*})\\&\qquad +\,(p+d+2\alpha )(((\beta -\gamma \alpha ) I_1^{2*}+2d+4\alpha )(d+2\alpha \\&\qquad +\,(\beta -\gamma \alpha ) I_1^{2*})+p(d+2\alpha )). \end{aligned}$$

Now if $\beta -\gamma \alpha >0$, and $\frac{ru_2}{(1+bu_2I_1^{2*})^2}>(\beta -\gamma \alpha )S_1^{2*}$, then $n_i>0, i=1,2,3$ and $n_1n_2>n_3$. Then, all the conditions of Routh–Hurwitz criteria are satisfied.

So, combining both the above two results, it may be concluded that the eigenvalues of $J(E_2^2)$ are all negative or have negative real part if $R_0^2>1$, $\beta >\gamma \alpha$ and $\frac{ru_2}{(1+bu_2I_1^{2*})^2}>(\beta +\gamma \alpha )S_1^{2*}$. Hence, the theorem is proved.