The effect of immigrant communities coming from higher incidence tuberculosis regions to a host country

TL;DR

This paper introduces a comprehensive TB model analyzing how immigrant communities from high-incidence regions affect TB dynamics in a host country, revealing that semi-closed communities can aid TB control and identifying optimal population distributions.

Contribution

The paper develops a new 25-variable TB mathematical model incorporating seasonal population flux and community structure, providing insights into TB spread and control strategies in immigrant populations.

Findings

Semi-closed communities can be beneficial for TB control.

Transmission coefficient variations significantly impact infection numbers.

An optimal community ratio minimizes the basic reproduction number R0.

Abstract

We introduce a new tuberculosis (TB) mathematical model, with $25$ state-space variables where $15$ are evolution disease states (EDSs), which generalises previous models and takes into account the (seasonal) flux of populations between a high incidence TB country (A) and a host country (B) with low TB incidence, where (B) is divided into a community (G) with high percentage of people from (A) plus the rest of the population (C). Contrary to some beliefs, related to the fact that agglomerations of individuals increase proportionally to the disease spread, analysis of the model shows that the existence of semi-closed communities are beneficial for the TB control from a global viewpoint. The model and techniques proposed are applied to a case-study with concrete parameters, which model the situation of Angola (A) and Portugal (B), in order to show its relevance and meaningfulness. Simulations show that variations of the transmission coefficient on the origin country has a big influence on the number of infected (and infectious) individuals on the community and the host country. Moreover, there is an optimal ratio for the distribution of individuals in (C) versus (G), which minimizes the reproduction number $R_0$. Such value does not give the minimal total number of infected individuals in all (B), since such is attained when the community (G) is completely isolated (theoretical scenario). Sensitivity analysis and curve fitting on $R_0$ and on EDSs are pursuit in order to understand the TB effects in the global statistics, by measuring the variability of the relevant parameters. We also show that the TB transmission rate $\beta$ does not act linearly on $R_0$, as is common in compartment models where system feedback or group interactions do not occur. Further, we find the most important parameters for the increase of each EDS.

Tables

Table 1. Table 1 : Estimated parameters and initial conditions values for Portugal and Angola.

Symbol	Description	Portugal	Angola
$\beta$	Transmission coefficient	variable ($72.358yr{s}^{-1}$)	variable ($150yr{s}^{-1}$)
$\nu$	Proportion of pulmonary TB cases	$0.75$	$0.937$
$\mu$	Natural death rate	$1/80yr{s}^{-1}$	$1/51yr{s}^{-1}$
$\delta$	Rate at which individuals leave P compartment	$2yr{s}^{-1}$	$2yr{s}^{-1}$
$\varphi$	Fraction of infected population developing active TB	$0.05$	$0.05$
$\sigma$	Reinfection (exogenous) factor for latent	$0.5$	$0.5$
$\omega$	Rate of endogenous reactivation for latent infections	$0.0003yr{s}^{-1}$	$0.0003yr{s}^{-1}$
$\tau$	Rate at which infectious individuals enter treatment	$4.26yr{s}^{-1}$	$2.13yr{s}^{-1}$
$k$	Proportion of detected cases in a year	$0.87$	$0.79$
${\delta }_T$	Inverse of treatment length	$1.36yr{s}^{-1}$	$1.36yr{s}^{-1}$
${\varphi }_T$	Fraction of treatment default and failure	$0.04$	$0.219$
$\eta$	Recruitment rate for Portugal	$78672$	$1287900$
$d_T$	TB induced death rate for Portugal	$1/5yr{s}^{-1}$	$1/8yr{s}^{-1}$
$N\left(0\right)$	Initial total population	10560000	24300000
$𝒮\left(0\right)$	Initial susceptible population	$8947300$	$9618729$
$𝒫\left(0\right)$	Initial primary infected with TB population	$11000$	$24300$
$ℐ\left(0\right)$	Initial actively infected (and infectious) population	$500$	$16164$
$ℒ\left(0\right)$	Initial latent infected population	$1600000$	$14580000$
$𝒯\left(0\right)$	Initial under treatment population	$1200$	$60807$

Table 2. Table 2 : Curve fitting of R 0 T superscript subscript 𝑅 0 𝑇 R_{0}^{T} .

Parameter	Type	Curve Fitting	${\log }_{10}\left(SQR\right)$
${\beta }^A$	best	$R_0^T=0.009021048+0.0422905212{\beta }^A+4.39929508E-7{\left({\beta }^A\right)}^2$
		$-7.99262125E-10{\left({\beta }^A\right)}^3$	$-5.278627639$
	as in $R_0^A$	$R_0^T=0.042398419{\beta }^A$	$-2.2700918$
${\beta }^C$	best	$R_0^T=6.3577692+2.2968396E-005{\beta }^C+6.545004E-008{\left({\beta }^C\right)}^2$
		$+6.933227073E-12{\left({\beta }^C\right)}^3+8.83463551E-13{\left({\beta }^C\right)}^4$	$-7.055330529$
	as in $R_0^C$	$R_0^T=0.08671477{\beta }^C$	$0.8692731$
${\beta }^G$	best	$R_0^T=\frac{-3639.13063363+1285.78172120{\beta }^G-4.168861981{\left({\beta }^G\right)}^2}{-572.32092506+202.21340477{\beta }^G-0.65585384{\left({\beta }^G\right)}^2}$	$-7.0129799$
	as in $R_0^G$	$R_0^T=0.0551009{\beta }^G$	$1.097186$
$k^C$	best	$R_0^T=P_{89}\left(k^C\right)$	$-7.058232835$
	not best	$R_0^T=\frac{-27.075657+53.2692494k^C+215.9660609{\left(k^C\right)}^2}{-4.2640340+8.36667709k^C+33.97750881{\left(k^C\right)}^2}$	$-6.6542336$
${\varphi }_T^G$	best	$R_0^T=P_{21}\left({\varphi }_T^G\right)$	$-7.007144237$
	not best	$R_0^T=\frac{1103.5908788-1522.87667500{\varphi }_T^G}{173.537338-239.505221{\varphi }_T^G}$	$-7.002127$
$a^A$	best	$R_0^T=\frac{38.07598473+2747.0900415a^A+42528.207079{\left(a^A\right)}^2}{6.0248383+419.985378a^A+6317.049255{\left(a^A\right)}^2}$	$-4.27964$
$a^B$	best	$R_0^T=6.78217128-77.5566782a^B+2571.6270531{\left(a^B\right)}^2$
		$-68307.97202742{\left(a^B\right)}^3+1194841.268572{\left(a^B\right)}^4-12711602.9588141{\left(a^B\right)}^5$
		$+73922994.9730162{\left(a^B\right)}^6-179395541.509671{\left(a^B\right)}^7$	$-2.1925449$
$\zeta$	best	$R_0^T=P_{91}\left(\zeta \right)$	$-7.11446753$
	not best	$R_0^T=6.383321621-0.1147570445\zeta +0.1218747529{\zeta }^2$	$-3.0198063$

Table 3. Table 3 : Qualitative sensitivity analysis.

	${ℬ}_{+-}$	${ℳ}_{+-}$	${ℰ}_{+-}$	${ℰ}_{++}$	${\xi }_{𝒜}$ values, respectively
$S_A$	$a^B$	${\beta }^A$			$1.17E-01,2.64E-01$
$P_A$			$a^B$	${\beta }^A$	$2.34E-01,3.79E-01$
$I_A$			$a^B$	${\beta }^A$	$2.28E-01,4.23E-01$
$L_A$			${\beta }^A$		$1.26E-01$
$T_A$			$a^B$	${\beta }^A$	$2.18E-01,4.58E-01$
$S_C$			$a^A$		$1.33E-01$
$P_C$			$k^C,a^B$	${\beta }^A,{\beta }^C,a^A,\zeta$	$2.47E-01,1.82E-01,3.78E-01,2.84E-01,8.43E-01,1$
$I_C$			${\beta }^C,k^C,a^B$	${\beta }^A,a^A,\zeta$	$1.81E-01,3.72E-01,1.84E-01,3.76E-01,8.72E-01,9.87E-01$
$L_C$			$a^A$		$2.43E-01$
$T_C$			${\beta }^C,a^B$	${\beta }^A,a^A,\zeta$	$1.39E-01,1.73E-01,3.97E-01,9.28E-01,9.74E-01$
$S_G$		$a^A$		$\zeta$	$3.95E-01,2.88E-01$
$P_G$			${\beta }^A,{\varphi }_T^G,a^A$	${\beta }^G,k^C,\zeta$	$1.80E-01,1.54E-01,2.26E-01,7.64E-01,4.21E-01,3.90E-01$
$I_G$			${\beta }^A,{\beta }^G,{\varphi }_T^G,\zeta$	$k^C,a^A$	$1.97E-01,6.53E-01,1.97E-01,2.15E-01,5.15E-01,3.75E-01$
$L_G$		$a^A$	${\beta }^G$	$\zeta$	$2.38E-01,9.46E-02,2.59E-01$
$T_G$			${\beta }^A,{\varphi }_T^G$	$k^C,a^A,\zeta$	$2.02E-01,1.72E-01,2.82E-01,4.56E-01,1.94E-01$

Equations

$$\left\{\begin{array}[]{l}\dot{\mathcal{S}}=\eta N-\left(\lambda(t)+\mu\right)\mathcal{S},\\ \dot{\mathcal{P}}=\lambda(t)\mathcal{S}+\sigma\lambda(t)\mathcal{L}-\left(\delta+\mu\right)\mathcal{P},\\ \dot{\mathcal{I}}=\phi\delta\mathcal{P}+\omega\mathcal{L}+\phi_{T}\delta_{T}\mathcal{T}-\left(\tau k+\mu\right)\mathcal{I},\\ \dot{\mathcal{L}}=(1-\phi)\delta\mathcal{P}+(1-\phi_{T})\delta_{T}\mathcal{T}-\left(\sigma\lambda(t)+\omega+\mu\right)\mathcal{L},\\ \dot{\mathcal{T}}=\tau k\mathcal{I}-\left(\delta_{T}+\mu\right)\mathcal{T}.\end{array}\right.$$

$$\dot{\lambda }=\beta \nu \left(\dot{\mathcal{I}}{N}^{-1}-\mathcal{I}{N}^{-2}\dot{N}\right).$$

$$\dot{\mathcal{X}}=\left(\beta \nu \mathcal{I}\mathcal{A}+\mathcal{B}\right)\mathcal{X}+\mathcal{C},$$

$$\mathcal{A}=\left(\begin{array}[]{ccccc}-1&0&0&0&0\\ 1&0&0&\sigma&0\\ 0&0&0&0&0\\ 0&0&0&-\sigma&0\\ 0&0&0&0&0\end{array}\right),\quad\mathcal{B}=\left(\begin{array}[]{ccccc}-\mu&0&0&0&0\\ 0&-(\delta+\mu)&0&0&0\\ 0&\phi\delta&-(\tau k+\mu)&\omega&\phi_{T}\delta_{T}\\ 0&(1-\phi)\delta&0&-(\omega+\mu)&(1-\phi_{T})\delta_{T}\\ 0&0&\tau k&0&-(\delta_{T}+\mu)\end{array}\right),$$

$$\left\{\begin{array}[]{l}\dot{\mathcal{S}}=\eta N-\left(\lambda(t)+\gamma(t)+\mu\right)\mathcal{S}+\tilde{\gamma}(t)\left((1-\zeta)\tilde{\mathcal{S}}_{1}+\zeta\tilde{\mathcal{S}}_{2}\right),\\ \dot{\mathcal{P}}=\lambda(t)\mathcal{S}+\sigma\lambda(t)\mathcal{L}-\left(\delta+\gamma(t)+\mu\right)\mathcal{P}+\tilde{\gamma}(t)\left((1-\zeta)\tilde{\mathcal{P}}_{1}+\zeta\tilde{\mathcal{P}}_{2}\right),\\ \dot{\mathcal{I}}=\phi\delta\mathcal{P}+\omega\mathcal{L}+\phi_{T}\delta_{T}\mathcal{T}-\left(\tau k+\gamma(t)+\mu\right)\mathcal{I}+\tilde{\gamma}(t)\left((1-\zeta)\tilde{\mathcal{I}}_{1}+\zeta\tilde{\mathcal{I}}_{2}\right),\\ \dot{\mathcal{L}}=(1-\phi)\delta\mathcal{P}+(1-\phi_{T})\delta_{T}\mathcal{T}-\left(\sigma\lambda(t)+\omega+\gamma(t)+\mu\right)\mathcal{L}+\tilde{\gamma}(t)\left((1-\zeta)\tilde{\mathcal{L}}_{1}+\zeta\tilde{\mathcal{L}}_{2}\right),\\ \dot{\mathcal{T}}=\tau k\mathcal{I}-\left(\delta_{T}+\gamma(t)+\mu\right)\mathcal{T}+\tilde{\gamma}(t)\left((1-\zeta)\tilde{\mathcal{T}}_{1}+\zeta\tilde{\mathcal{T}}_{2}\right).\end{array}\right.$$

$$\left\{\begin{array}[]{l}\dot{\mathcal{S}}=\eta N-\left(\lambda+\gamma(t)+\mu\right)\mathcal{S}+\tilde{\gamma}(t)\left((1-\zeta)\tilde{\mathcal{S}}_{1}+\zeta\tilde{\mathcal{S}}_{2}\right),\\ \dot{\mathcal{P}}=\lambda\mathcal{S}+\sigma\lambda\mathcal{L}-\left(\delta+\gamma(t)+\mu\right)\mathcal{P}+\tilde{\gamma}(t)\left((1-\zeta)\tilde{\mathcal{P}}_{1}+\zeta\tilde{\mathcal{P}}_{2}\right),\\ \dot{\mathcal{I}}=\phi\delta\mathcal{P}+\omega\mathcal{L}+\phi_{T}\delta_{T}\mathcal{T}-\left(\tau k+\gamma(t)+\mu\right)\mathcal{I}+\tilde{\gamma}(t)\left((1-\zeta)\tilde{\mathcal{I}}_{1}+\zeta\tilde{\mathcal{I}}_{2}\right),\\ \dot{\mathcal{L}}=(1-\phi)\delta\mathcal{P}+(1-\phi_{T})\delta_{T}\mathcal{T}-\left(\sigma\lambda+\omega+\gamma(t)+\mu\right)\mathcal{L}+\tilde{\gamma}(t)\left((1-\zeta)\tilde{\mathcal{L}}_{1}+\zeta\tilde{\mathcal{L}}_{2}\right),\\ \dot{\mathcal{T}}=\tau k\mathcal{I}-\left(\delta_{T}+\gamma(t)+\mu\right)\mathcal{T}+\tilde{\gamma}(t)\left((1-\zeta)\tilde{\mathcal{T}}_{1}+\zeta\tilde{\mathcal{T}}_{2}\right),\\ \dot{N}=\left(\eta-\gamma(t)-\mu\right)N+\tilde{\gamma}(t)\left((1-\zeta)\tilde{N}_{1}+\zeta\tilde{N}_{2}\right).\end{array}\right.$$

$$\dot{S}$$

$$=\dot{\mathcal{S}}{N}^{-1}-S{N}^{-1}\left(\left(\eta -\gamma (t)-\mu \right)N+\tilde{\gamma }(t)\left((1-\zeta ){\tilde{N}}_1+\zeta {\tilde{N}}_2\right)\right)$$

$$=\dot{\mathcal{S}}{N}^{-1}-\left(\eta -\gamma (t)-\mu +\tilde{\gamma }(t)\left((1-\zeta ){\tilde{N}}_1+\zeta {\tilde{N}}_2\right)N^{-1}\right)S,$$

$$=\dot{\mathcal{S}}{N}^{-1}-\left(M(t)-\gamma (t)-\mu \right)S,$$

$$\left\{\begin{array}[]{l}\dot{S}=\eta-\left(\lambda+M(t)\right)S+\tilde{\gamma}(t)\left((1-\zeta)\tilde{S}_{1}+\zeta\tilde{S}_{2}\right),\\ \dot{P}=\lambda S+\sigma\lambda L-\left(\delta+M(t)\right)P+\tilde{\gamma}(t)\left((1-\zeta)\tilde{P}_{1}+\zeta\tilde{P}_{2}\right),\\ \dot{I}=\phi\delta P+\omega L+\phi_{T}\delta_{T}T-\left(\tau k+M(t)\right)I+\tilde{\gamma}(t)\left((1-\zeta)\tilde{I}_{1}+\zeta\tilde{I}_{2}\right),\\ \dot{L}=(1-\phi)\delta P+(1-\phi_{T})\delta_{T}T-\left(\sigma\lambda+\omega+M(t)\right)L+\tilde{\gamma}(t)\left((1-\zeta)\tilde{L}_{1}+\zeta\tilde{L}_{2}\right),\\ \dot{T}=\tau kI-\left(\delta_{T}+M(t)\right)T+\tilde{\gamma}(t)\left((1-\zeta)\tilde{T}_{1}+\zeta\tilde{T}_{2}\right),\\ \dot{\lambda}=\beta\nu\dot{I}=\beta\nu\left(\phi\delta P+\omega L+\phi_{T}\delta_{T}T-\left(\tau k+M(t)\right)I+\tilde{\gamma}(t)\left((1-\zeta)\tilde{I}_{1}+\zeta\tilde{I}_{2}\right)\right),\\ \dot{N}=\left(M(t)-\gamma(t)-\mu\right)N.\end{array}\right.$$

$${\gamma }_A(t)=a^A+b^{A}t+a^A{p}^A\cos \left({\theta }^{A}t\right),\text{mbox}and{\gamma }_B(t)=a^B+b^{B}t+a^B{p}^B\cos \left({\theta }^{B}t\right),$$

$$\left\{\begin{array}[]{l}\dot{\gamma}_{A}=z^{A},\\ \dot{z}_{A}=-(\theta^{A})^{2}(\gamma_{A}-a^{A}-b^{A}\,t),\end{array}\right.\quad\mbox{ and }\quad\left\{\begin{array}[]{l}\dot{\gamma}_{B}=z^{B},\\ \dot{z}_{B}=-(\theta^{B})^{2}(x-a^{B}-b^{B}\,t),\end{array}\right.$$

$${\dot{V}}_N=\mathcal{A}(t)V_N,$$

$$\mathcal{A}(t)=\left(\begin{array}[]{ccc}\eta^{A}-\mu^{A}-\gamma_{A}(t)&\gamma_{B}(t)(1-\zeta)&\gamma_{B}(t)\zeta\\ \gamma_{A}(t)(1-\zeta)&\eta^{C}-\mu^{C}-\gamma_{B}(t)&0\\ \gamma_{A}(t)\zeta&0&\eta^{C}-\mu^{C}-\gamma_{B}(t)\end{array}\right).$$

$$\left\{\begin{array}[]{l}\dot{S}_{A}=\eta^{A}-\left(\lambda_{A}+M_{A}\right)S_{A}+\gamma_{B}\left((1-\zeta)S_{C}+\zeta S_{G}\right),\\ \dot{P}_{A}=\lambda_{A}S_{A}+\sigma^{A}\lambda_{A}L_{A}-\left(\delta^{A}+M_{A}\right)P_{A}+\gamma_{B}\left((1-\zeta)P_{C}+\zeta P_{G}\right),\\ \dot{I}_{A}=\phi^{A}\delta^{A}P_{A}+\omega^{A}L_{A}+\phi^{A}_{T}\delta^{A}_{T}T_{A}-\left(\tau^{A}k^{A}+M_{A}\right)I_{A}+\gamma_{B}\left((1-\zeta)I_{C}+\zeta I_{G}\right),\\ \dot{L}_{A}=(1-\phi^{A})\delta^{A}P_{A}+(1-\phi^{A}_{T})\delta^{A}_{T}T_{A}-\left(\sigma^{A}\lambda_{A}+\omega^{A}+M_{A}\right)L_{A}+\gamma_{B}\left((1-\zeta)L_{C}+\zeta L_{G}\right),\\ \dot{T}_{A}=\tau^{A}k^{A}I_{A}-\left(\delta^{A}_{T}+M_{A}\right)T_{A}+\gamma_{B}\left((1-\zeta)T_{C}+\zeta T_{G}\right),\\ \dot{\lambda}_{A}=\beta^{A}\nu^{A}\left(\phi^{A}\delta^{A}P_{A}+\omega^{A}L_{A}+\phi^{A}_{T}\delta^{A}_{T}T_{A}-\left(\tau^{A}k^{A}+M_{A}\right)I^{A}+\gamma_{B}\left((1-\zeta)I_{C}+\zeta I_{G}\right)\right),\\ \dot{N}_{A}=\left(M_{A}-\gamma_{A}-\mu^{A}\right)N_{A},\\ \end{array}\right.$$

$$\left\{\begin{array}[]{l}\dot{S}_{G}=\eta^{C}-\left(\lambda_{G}+M_{G}\right)S_{G}+\gamma_{A}\zeta S_{A},\\ \dot{P}_{G}=\lambda_{G}S_{G}+\sigma^{C}\lambda_{G}L_{G}-\left(\delta^{C}+M_{G}\right)P_{G}+\gamma_{A}\zeta P_{A},\\ \dot{I}_{G}=\phi^{C}\delta^{C}P_{G}+\omega^{C}L_{G}+\phi^{G}_{T}\delta^{C}_{T}T_{G}-\left(\tau^{C}k^{C}+M_{G}\right)I_{G}+\gamma_{A}\zeta I_{A},\\ \dot{L}_{G}=(1-\phi^{C})\delta^{C}P_{G}+(1-\phi^{G}_{T})\delta^{C}_{T}T_{G}-\left(\sigma^{C}\lambda_{G}+\omega^{C}+M_{G}\right)L_{G}+\gamma_{A}\zeta L_{A},\\ \dot{T}_{G}=\tau^{C}k^{C}I_{G}-\left(\delta^{C}_{T}+M_{G}\right)T_{G}+\gamma_{A}\zeta T_{A},\\ \dot{\lambda}_{G}=\beta^{G}\nu^{G}\left(\phi^{C}\delta^{C}P_{G}+\omega^{C}L_{G}+\phi^{G}_{T}\delta^{C}_{T}T_{G}-\left(\tau^{C}k^{C}+M_{G}\right)I^{G}+\gamma_{A}\zeta I_{A}\right),\\ \dot{N}_{G}=\left(M_{G}-\gamma_{B}-\mu^{C}\right)N_{G},\\ \end{array}\right.$$

$$\left\{\begin{array}[]{l}\dot{S}_{C}=\eta^{C}-\left(\lambda_{C}+M_{C}\right)S_{C}+\gamma_{A}(1-\zeta)S_{A},\\ \dot{P}_{C}=\lambda_{C}S_{C}+\sigma^{C}\lambda_{C}L_{C}-\left(\delta^{C}+M_{C}\right)P_{C}+\gamma_{A}(1-\zeta)P_{A},\\ \dot{I}_{C}=\phi^{C}\delta^{C}P_{C}+\omega^{C}L_{C}+\phi^{C}_{T}\delta^{C}_{T}T_{C}-\left(\tau^{C}k^{C}+M_{C}\right)I_{C}+\gamma_{A}(1-\zeta)I_{A},\\ \dot{L}_{C}=(1-\phi^{C})\delta^{C}P_{C}+(1-\phi^{C}_{T})\delta^{C}_{T}T_{C}-\left(\sigma^{C}\lambda_{C}+\omega^{C}+M_{C}\right)L_{C}+\gamma_{A}(1-\zeta)L_{A},\\ \dot{T}_{C}=\tau^{C}k^{C}I_{C}-\left(\delta^{C}_{T}+M_{C}\right)T_{C}+\gamma_{A}(1-\zeta)T_{A},\\ \dot{\lambda}_{C}=\beta^{C}\nu^{C}\left(\phi^{C}\delta^{C}P_{C}+\omega^{C}L_{C}+\phi^{C}_{T}\delta^{C}_{T}T_{C}-\left(\tau^{C}k^{C}+M_{C}\right)I^{C}+\gamma_{A}(1-\zeta)I_{A}\right),\\ \dot{N}_{C}=\left(M_{C}-\gamma_{B}-\mu^{C}\right)N_{C},\\ \end{array}\right.$$

$$\left\{\begin{array}[]{l}\dot{\gamma}_{A}=z^{A},\\ \dot{z}_{A}=-(\theta^{A})^{2}(\gamma_{A}-a^{A}-b^{A}\,t),\\ \dot{\gamma}_{B}=z^{B},\\ \dot{z}_{B}=-(\theta^{B})^{2}(\gamma_{B}-a^{B}-b^{B}\,t),\end{array}\right.$$

$$M_A$$

$$M_C$$

$$M_G$$

$${\dot{N}}_A+{\dot{N}}_C+{\dot{N}}_G$$

$$R_0=\frac{\beta \nu \delta ({\delta }_T+\mu )(\varphi \mu +\omega )}{\mu (\delta +\mu )[(\mu +\omega )(\tau k+{\delta }_T+\mu )+{\delta }_T\tau k(1-{\varphi }_T)]},$$

$$R_{0}^{A}=$6.784924946$,\quad R_{0}^{C}=$1.116995163$,\quad R_{0}^{G}=$2.365451295$,$$

$$x=(P_A,P_C,P_G,I_A,I_C,I_G,L_A,L_C,L_G,T_A,T_C,T_G,S_A,S_C,S_G)\in {\mathbb{R}}_{+}^{15}.$$

$$\left\{\begin{array}[]{l}\eta^{A}-M_{A}S_{A}+a^{B}\left((1-\zeta)S_{C}+\zeta S_{G}\right)=0,\\ \eta^{C}-M_{C}S_{C}+a^{A}(1-\zeta)S_{A}=0,\\ \eta^{G}-M_{G}S_{G}+a^{A}\zeta S_{A}=0,\\ -\left(\delta^{A}+M_{A}\right)P_{A}+a^{B}\left((1-\zeta)P_{C}+\zeta P_{G}\right)=0,\\ -\left(\delta^{C}+M_{C}\right)P_{C}+a^{A}(1-\zeta)P_{A}=0,\\ -\left(\delta^{C}+M_{G}\right)P_{G}+a^{A}\zeta P_{A}=0,\\ \phi^{A}\delta^{A}P_{A}+\omega^{A}L_{A}+\phi^{A}_{T}\delta^{A}_{T}T_{A}=0,\\ \phi^{C}\delta^{C}P_{C}+\omega^{C}L_{C}+\phi^{C}_{T}\delta^{C}_{T}T_{C}=0,\\ \phi^{C}\delta^{C}P_{G}+\omega^{C}L_{G}+\phi^{G}_{T}\delta^{C}_{T}T_{G}=0,\\ (1-\phi^{A})\delta^{A}P_{A}+(1-\phi^{A}_{T})\delta^{A}_{T}T_{A}-\left(\omega^{A}+M_{A}\right)L_{A}+a^{B}\left((1-\zeta)L_{C}+\zeta L_{G}\right)=0,\\ (1-\phi^{C})\delta^{C}P_{C}+(1-\phi^{C}_{T})\delta^{C}_{T}T_{C}-\left(\omega^{C}+M_{C}\right)L_{C}+a^{A}(1-\zeta)L_{A}=0,\\ (1-\phi^{C})\delta^{C}P_{G}+(1-\phi^{G}_{T})\delta^{C}_{T}T_{G}-\left(\omega^{C}+M_{G}\right)L_{G}+a^{A}\zeta L_{A}=0,\\ -\left(\delta^{A}_{T}+M_{A}\right)T_{A}+a^{B}\left((1-\zeta)T_{C}+\zeta T_{G}\right)=0,\\ -\left(\delta^{C}_{T}+M_{C}\right)T_{C}+a^{A}(1-\zeta)T_{A}=0,\\ -\left(\delta^{C}_{T}+M_{G}\right)T_{G}+a^{A}\zeta T_{A}=0.\end{array}\right.$$

$$\left(\begin{array}[]{ccc}-\delta^{A}_{T}-M_{A}&a^{B}(1-\zeta)&a^{B}\zeta\\ a^{A}(1-\zeta)&-\delta^{C}_{T}-M_{C}&0\\ a^{A}\zeta&0&-\delta^{G}_{T}-M_{G}\end{array}\right)\left(\begin{array}[]{c}T_{A}\\ T_{C}\\ T_{G}\end{array}\right)=\left(\begin{array}[]{c}0\\ 0\\ 0\end{array}\right)\>\>\Rightarrow\>\>T_{A}=T_{C}=T_{G}=0.$$

$$M_A{S}_A$$

$$M_C{S}_C$$

$$M_G{S}_G$$

$$x_0=(0,0,0,0,0,0,0,0,0,0,0,0,1,1,1),$$

$$X_s=\{(0,0,0,0,0,0,0,0,0,0,0,0,S_A,S_C,S_G)\in {\mathbb{R}}^{15}:S_A,S_C,S_G\ge 0\}.$$

$$\mathcal{F}_{1-3}(x)=\left(\begin{array}[]{c}\beta^{A}\nu^{A}I_{A}\left(S_{A}+\sigma^{A}L_{A}\right)+a^{B}\left((1-\zeta)P_{C}+\zeta P_{G}\right)\\ \beta^{C}\nu^{C}I_{C}\left(S_{C}+\sigma^{C}L_{C}\right)+a^{A}(1-\zeta)P_{A}\\ \beta^{G}\nu^{G}I_{G}\left(S_{G}+\sigma^{C}L_{G}\right)+a^{A}\zeta P_{A}\end{array}\right),\quad\mathcal{F}_{j}(x)=0\mbox{ for }j\in\{4,\cdots,15\},$$

Full text

Keywords:

tuberculosis; mathematical model; flux of populations; sensitivity analysis; curve fitting; reproduction number.

Mathematics Subject Classification 2010:

92D30.

1 Introduction

Tuberculosis (TB) is an infectious disease caused by the Mycobacterium tuberculosis (Mtb). Following the World Health Organization (WHO), the (Mtb) is the second cause of death worldwide from a single infectious agent, after the human immunodeficiency virus [29]. TB is present in all regions of the world. Most of the estimated number of cases in 2013 occurred in Asia ($56\%$) and the African region ($29\%$); smaller proportions of cases occurred in the Eastern Mediterranean region ($8\%$), the European region (4%) and the region of the Americas ($3\%$) [30].

In TB spread, migration plays an important role, e.g., following the International Organization for Migration (IOM), TB is a social disease and migration, as a social determinant of health, increases TB-related morbidity and mortality among migrants and surrounding communities [10]. Migrants of specific legal and social status, such as workers, undocumented migrants, trafficked and detained persons, face particular TB vulnerabilities. Among migrant workers with a legal status, their access to TB diagnosis and care is subject to their ability to access health care services and health insurance coverage, provided either by the state or the employer. Illegal migrants face particular challenges such as fear of deportation that delay or limit their access to diagnostic and treatment services. Deportation while on treatment or poor compliance with treatment may lead to drug resistant infection and increased chances of spreading TB in countries of origin, transit and destination [10].

Mathematical models are an important tool in analyzing the spread and control of infectious diseases [7, 8]. There are many mathematical dynamic models for TB, see, e.g., [1, 3, 4, 6, 27] and references cited therein. There are also models dedicated to study TB transmission dynamics in immigrants and local population. Usually, these models divide the total population into two subgroups: immigrants and local subpopulation. Each subgroup is divided into several epidemiological compartments: susceptible, latent, infectious, recovered, or other, depending on the type of the model, see, e.g., [2, 11, 32, 33]. In general, compartment models written with ordinary differential equations tend to be nice approximations of the true scenario that have rather simple formulation, e.g., with five state-space variables and a (non)autonomous quadratic vector field, because of numerical and analytic limitations and the tradeoff between complexity and the relevant information that they can present. In particular, heterogeneous situations may be studied using such models. However, no interaction between individuals in the different groups are considered in such models. We are interested in understanding how the flux and distribution of individuals affects TB on a host country. As a case-study, we have considered the situation of Angola and Portugal, although the techniques may be applied to any similar situation.

Angola is the seventh-largest country in Southern Africa with a total population of approximately $24.3$ million [9]. WHO predicts that by 2017 the TB cases rate may rise significantly in Angola. A natural question is to try to understand how this may affect the rest of the world. According to Celestino Teixeira, the Coordinator of the Fight Against Tuberculosis Programme, in 2013 Angola reported a total of $60,807$ cases of TB in all forms, observing an increase of $11\%$ over the previous year [39]. Portugal is a country in Southwest Europe with a total population of approximately $10.5$ million [9]. In 2014, for the first time, the incidence of TB in Portugal was estimated to be lower than $20$ new cases per $100,000$ inhabitants, placing Portugal among the countries with low TB incidence. However, there are still some regions (Lisbon and Porto) with much higher TB incidences [17]. Portugal is a relevant geographically area of study for TB because its infection behaviour is not similar to the rest of Europe, in the sense that has higher incidence of tuberculosis. Aside from the independence period, Angola is characterized by a reduced emigration and is becoming gradually an attractive region, receiving migrants from different regions, including Portugal [19]. Following the Portuguese Emigration Observatory, in 2014 there were 126,356 Portuguese emigrants living in Angola [40]. According to the Organisation for Economic Co-operation and Development (OECD) [16], for the first time in five years, 2012 saw the number of long-term entry visas grow. Visas to Angolans doubled in 2012, mainly for study. According to the Portuguese Foreigners and Borders Service, in 2012 there were 20,177 Angolans citizens living in Portugal [21]. Although Angolans living in Portugal are dispersed throughout the country, there is a very high concentration in the district of Lisbon, followed by Setúbal and Porto [15].

In this paper, we propose and study a new mathematical model for TB that generalises the one proposed in [13]. We consider three different populations: people living in a high TB incidence country (A), people living in a low TB incidence country in a semi-closed community of the high incidence country natives (G), and the other persons living in the low incidence country (C). Each of these three groups of population are subdivided into the five epidemiological categories considered in the model from [13]. Our model considers the movement of persons from the high TB incidence country to the low TB incidence country and vice-versa. We assume that the individuals that arrive and depart from the low TB incidence country are split into the ones that enter/leave the semi-closed community of the high TB incidence country natives and the ones that enter/leave other regions of the low TB incidence country. Our model is quite different from [13] and other TB models in the literature, since it has internal transfer of individuals between the subgroups, high TB incidence country, semi-closed community of high TB incidence country natives and other persons living in the low TB incidence country. We consider a case study where the low TB incidence country is represented by Portugal and the high TB incidence country is represented by Angola.

The paper is organized as follows. In Section 2, we explain how we construct our model. The basic reproduction number is algebraically and numerically computed in Section 3 for the autonomous case. This section also includes a sensitivity analysis of the basic reproduction number with respect to TB transmission rates, transfer of individuals and ratio of individuals that stay in the community versus spread in the host country. Section 4 is devoted to numerical simulations, which help us to make a qualitative sensitivity analysis for each epidemiological category of the subgroups Angola, semi-closed community of Angola natives and other persons living in Portugal, when relevant TB parameters are perturbed. We end with Section 5 of conclusions and future work.

2 Mathematical model

We construct a model with three components, based on [13], where there exists seasonal flux of population between some of the components. The model from [13] divides the total population $N$ in five epidemiological compartments: susceptible individuals ($\mathcal{S}$) that never have been in contact with (Mtb), primary infected individuals ($\mathcal{P}$) that have been infected by (Mtb) but it is not certain if the disease will progress, actively infected and infectious individuals ($\mathcal{I}$) that are not yet in treatment, latent infected individuals ($\mathcal{L}$) and under treatment individuals ($\mathcal{T}$). Susceptible individuals become primary infected at a rate $\lambda=\beta\nu\mathcal{I}$ $yrs^{-1}$, where $\beta$ is the transmission coefficient and $\nu$ is the proportion of pulmonary TB cases. A proportion $\phi$ and $(1-\phi)$ of individuals in the class $\mathcal{P}$ is transferred to the class $\mathcal{I}$ and $\mathcal{L}$, respectively, at a rate $\delta\,yrs^{-1}$. Each year, a proportion $k$ of individuals in the class $\mathcal{I}$ is detected and start TB treatment at a rate $\tau\,yrs^{-1}$, entering the class $\mathcal{T}$. It is assumed that individuals in the class $\mathcal{T}$ are neither infectious nor susceptible to reinfection. A fraction $\phi_{T}$ of individuals in class $\mathcal{T}$ is transferred to class $\mathcal{I}$ due to either treatment failure or default, while the remaining $(1-\phi_{T})$ are successfully treated and enter in the class $\mathcal{L}$. The inverse of treatment length is denoted by $\delta_{T}$. In [13], birth and death rates are assumed equal, here we assume that they can be different and we denote the recruitment rate by $\eta\,yrs^{-1}$ and the death rate by $\mu\,yrs^{-1}$. The reinfection factor is denoted by $\sigma$ (see [13] for more details). Optimal control strategies for such model were studied in [20, 23, 24].

Let $\mathcal{S}\equiv\mathcal{S}(t)$, $\mathcal{P}\equiv\mathcal{P}(t)$, $\mathcal{I}\equiv\mathcal{I}(t)$, $\mathcal{L}\equiv\mathcal{L}(t)$, $\mathcal{T}\equiv\mathcal{T}(t)$, where $t$ represents time in years. The model described above is given by the following system of ordinary differential equations:

[TABLE]

We have $N=\mathcal{S}+\mathcal{P}+\mathcal{I}+\mathcal{L}+\mathcal{T}$ and $\lambda(t)=\beta\nu\mathcal{I}N^{-1}$. Then

[TABLE]

On the other hand, $\dot{N}=(\eta-\mu)N$, so if $\eta=\mu$ then the population is constant. The system can be written in a matrix form as

[TABLE]

where $\mathcal{X}=(\mathcal{S},\mathcal{P},\mathcal{I},\mathcal{L},\mathcal{T})$,

[TABLE]

and $\mathcal{C}=(\eta N,0,0,0,0)$. We can verify that the matrix $\lambda(t)\mathcal{A}+\mathcal{B}$ can be diagonalizable, so there is a semi-closed form solution for the problem (it is not closed a priori because $\lambda$ still depends on $I$ and $N$).

Suppose this system interacts with (a convex combination of) another two similar systems $\tilde{X}_{1}$ and $\tilde{X}_{2}$, in the following way: there exist functions $\gamma(t),\tilde{\gamma}(t)\in[0,1]$ and a value $\zeta\in[0,1]$ such that

[TABLE]

Adding $N=\mathcal{S}+\mathcal{P}+\mathcal{I}+\mathcal{L}+\mathcal{T}$ as a new state variable, we have

[TABLE]

Let $S=\mathcal{S}N^{-1}$, $P=\mathcal{P}N^{-1}$, $I=\mathcal{I}N^{-1}$, $L=\mathcal{L}N^{-1}$, $T=\mathcal{T}N^{-1}$. These variables now represent the percentage of the population in each state, i.e., $S+P+I+L+T=1$. Since

[TABLE]

with $M(t)\mathop{\stackrel{{\scriptstyle\mbox{\tiny def}}}{{=}}}\eta+\left((1-\zeta)\tilde{N}_{1}+\zeta\tilde{N}_{2}\right)\tilde{\gamma}(t)N^{-1}$, where the calculations for the other variables are similar, and adding $\lambda(t)=\beta\nu I$ as a new state variable, we have

[TABLE]

Using the above model, we consider different population groups: people living in a high incidence TB country (A) and people living in a low incidence TB country (B), where (B) is subdivided in a community (G) with high percentage of people from (A), and (C) is the rest of the population of (B). We consider that the values of $\beta$, $\nu$, $\phi_{T}$ of the group (G) are different from the values of the group (C). The flux of population follow the distribution functions $\gamma_{A}$, from (A) to (B), and $\gamma_{B}$, from (B) to (A). We assume that the persons that arrive and departure from (B) are split in the following proportions: $\zeta$ goes to (G) and $(1-\zeta)$ goes to (C), with $\zeta\in[0,1]$ a fixed percentage value in this model.

This model accounts for an average moving value of persons $a^{A}$, $a^{B}$ that increases/decreases in time by the slopes $b^{A}$, $b^{B}$ and has a seasonality variation modeled by $p^{A}$, $p^{B}$, $\theta^{A}$, $\theta^{B}$. The flux of population will be modeled by the following functions:

[TABLE]

for constants $a^{A},a^{B},b^{A},b^{B},p^{A},p^{B},\theta^{A},\theta^{B}\in\mathbb{R}$ chosen to ensure that $0\leq\gamma_{A}(t),\gamma_{B}(t)\leq 1$ for all $t$ of the simulation.

The flux of population $\gamma_{A}(t)$, $\gamma_{B}(t)$ can be incorporated as state-space variables. In our case, the functions $\gamma^{A}$, $\gamma^{P}$ are solutions of the system of ODEs

[TABLE]

which we add to the model (2.8)–(2.11), obtaining the complete model with $25$ state-space variables. Note that if $V_{N}=(N_{A},N_{C},N_{G})$, then

[TABLE]

where

[TABLE]

So the population evolution is only dependent on the moving distribution functions $\gamma^{A}$, $\gamma^{P}$, born rates $\eta$, and natural death rates $\mu$. Hence, we obtain the complete model composed by the four subsystems (2.8)–(2.11) composed by: (i) the variables of the high incidence TB country

[TABLE]

(ii) the variables associated with the community in the host country

[TABLE]

(iii) the variables related with the population of the host country excluding the community

[TABLE]

(iv) and the variables measuring the flux of population

[TABLE]

where for presentation convenience we define

[TABLE]

Note that

[TABLE]

Again, if $\eta^{A}=\mu^{A}$ and $\eta^{C}=\mu^{C}$, then the total population is constant. Moreover, if $b^{A}=b^{B}=p^{A}=p^{B}=0$, then system (2.8)–(2.11) is autonomous. For notation clarity, all parameters (i.e., constant values) have upper indices whereas state variables have lower indices.

3 Reproduction number and its sensitivity analysis for the autonomous case

The transmissibility of an infection can be asymptotically quantified by its reproduction number $R_{0}$ (for autonomous models), defined as the mean number of secondary infections seeded by a typical infective into a susceptible population. Since $R_{0}$ is a condition for the asymptotic stability of solutions around a free disease equilibrium point, this value determines a threshold: whenever $R_{0}>1$, a typical infective gives rise, on average, to more than one secondary infection, leading to an epidemic. In contrast, when $R_{0}<1$, infectious typically give rise, on average, to less than one secondary infection, and the prevalence of infection cannot increase.

A key point is that the model (2.8)–(2.11) is a priori nonautonomous, due to the flux of population $\gamma_{A}$ and $\gamma_{B}$. For such reason, from now on we assume that $\gamma_{A}(t)\equiv a^{A}$ and $\gamma_{B}(t)=a^{B}$, i.e., $b^{A}=b^{B}=p^{A}=p^{B}=0$ in (2.6), so that model (2.8)–(2.11) becomes autonomous and we can apply the standard method from [26]. A complete nonautonomous situation will be considered in a future work.

The reproducing number $R_{0}$ of system (2.1) can be analytically determined and, when $\eta=\mu$, is given by

[TABLE]

see, e.g., [13]. Hence, $R_{0}$ is proportional to $\beta$, $\nu$, $\phi$, $\phi_{T}$ ($0<\phi_{T}<1$) and inverse proportional to $\tau$ and $k$. In the no-transfer situation, i.e., $\gamma_{A}\equiv\gamma_{B}\equiv 0$, our model reduces to the disjoint coupling of the (sub)systems $(A)$, $(C)$ and $(G)$ similar to (2.1), so we can compute the reproduction numbers for the subsystems (using the fixed parameters from Table 1) in the no-transfer situation using (3.1), which gives

[TABLE]

where $R_{0}^{A}$, $R_{0}^{C}$ and $R_{0}^{G}$ denote the basic reproduction number for populations (A), (C) and (G), respectively, when they are complete independent from each others (no flux of population between the compartments). For the complete system (2.8)–(2.11) the basic reproduction number will be denoted by $R_{0}^{T}$. Note that the coupling of only $(C)$ and $(G)$ (again in the no-transfer situation and without the components associated to $(A)$) is known in the literature as a model for heterogeneous infection risk [5, 13].

The complete system (2.8)–(2.11), although a generalization of previous models, is quite different from systems like (2.1), by the fact that it has internal transfer of individuals between subsystems $(A)$ and $(C)$ and $(G)$, so it is not expected that $R_{0}^{T}$ follows the same expression (3.1). So its relevant to understand how $R_{0}^{T}$ is affected by variation of the parameters. In order to verify the validity and to obtain the value of $R_{0}^{T}$, depending on the parameters chosen, we follow the approach in [26].

Let $x$ represent the state-space variables (in a special order) that group the individuals in each disease state and group compartment, i.e.,

[TABLE]

Note that there exists an equilibrium point with $I_{A},I_{C},I_{G}=0$, if $\lambda_{A}=\lambda_{B}=\lambda_{C}=0$ and

[TABLE]

From the last three equations, we have

[TABLE]

In the same way we can see, from fourth to sixth equations, that $P_{A}=P_{C}=P_{G}=0$ and, from the other equations, that $L_{A}=L_{C}=L_{G}=0$. Since, $\eta^{A},\eta^{C}\neq 0$ and

[TABLE]

from the first three equations, we have $S_{A}=S_{C}=S_{G}=1$. Hence, the disease free equilibrium point (DFE) is unique and given by

[TABLE]

and it makes sense to define the set of all disease free states $X_{s}$ as

[TABLE]

In our model the individuals get the first contact with the infection in the states $P_{A},P_{C},P_{G}$. We have $m=12$ states where individuals have different degrees of infection and $3$ states free of disease. The vector field $X$ in (2.8)–(2.11) is now divided as $X=\mathcal{F}-(\mathcal{V}^{-}-\mathcal{V}^{+})$, where $\mathcal{F}$ is the rate of appearance of new infections, $\mathcal{V}^{+}$ is the rate of in-transfers of individuals by other means, and $\mathcal{V}^{-}$ is the rate of out-transfers of individuals by other means. We have

[TABLE]

[TABLE]

Note that $\mathcal{F}_{1-3}$ denotes the entries of $\mathcal{F}$ from $1$ to $3$. Then $\mathcal{F}$ and $\mathcal{V}=\mathcal{V}^{+}-\mathcal{V}^{-}$ satisfy the following assumptions:

• $(A_{1})$

if $x\geq 0$, then $\mathcal{F}(x)$, $\mathcal{V}^{+}(x)$, $\mathcal{V}^{-}(x)\geq 0$ (each function represents a direct transfer of individuals);

• $(A_{2})$

if $x_{i}=0$, then $\mathcal{V}_{i}^{-}(x)=0$ (if the compartment is empty, then there cannot be out-transfers of individuals);

• $(A_{3})$

$\mathcal{F}_{i}(x)=0$ for $i>12$;

• $(A_{4})$

if $x\in X_{s}$, then $\mathcal{F}_{i}(x)=0$ and $\mathcal{V}_{i}^{+}(x)=0$ for $1\leq i\leq 12$ (if the population is free of disease, then it will remain free of disease);

• $(A_{5})$

when $\mathcal{F}(x)=0$ we have that $DX(x_{0})$ is a Hurwitz matrix, i.e., all eigenvalues have negative real part (the equilibrium point $x_{0}$ is asymptotically stable).

Only assumption $(A_{5})$ creates some difficulty, since the other assumptions are evident. We numerically checked $(A_{5})$ (in all calculations made) using the Routh–Hurwitz criterion, which states that the matrix $A=DX(x_{0})$ is Hurwitz if and only if all the principal subdeterminants, of a special matrix constructed with the coefficients of the characteristic polynomial of $A$, are all strictly positive.

By Lemma 1 in [26], the derivatives $D\mathcal{F}(x_{0})$ and $D\mathcal{V}(x_{0})$ are partitioned as

[TABLE]

where $F$ and $V$ are $m\times m$-matrices. Hence, we have $F_{i,j}(x)=0$, if $i>m$ or $j>m$, and

[TABLE]

The critical threshold function $R_{0}^{T}$ is then given as the spectral radius of the matrix $A=FV^{-1}$. We have that $A$ has all entries zero except

[TABLE]

Considering the algebraic complexity of computing the spectral radius of $A$, in the next subsection we proceed numerically by understanding $R_{0}$ from the variation of the parameters.

3.1 Sensitivity analysis: numerical simulations

The values of the parameters $\beta$, $\nu$, $\mu$, $\delta$, $\phi$, $\sigma$, $\omega$, $\tau$, $k$, $\delta_{T}$ and $\phi_{T}$ estimated for Portugal, are based on the values proposed in [13], as well as the initial conditions $N(0)$, $\mathcal{S}(0)$, $\mathcal{P}(0)$, $\mathcal{L}(0)$, $\mathcal{I}(0)$, $\mathcal{T}(0)$. We assume that the Portuguese total population will decrease ($\eta<\mu N$), based on the projections for resident population in Portugal from Statistics Portugal [9] and the value for TB induced death that comes from [25].

We assume that the reference value for the transmission coefficient in Angola is $\beta=150$ based on [37]. According to the World Bank, the natural death rate in Angola is equal to $\mu=1/51\,yrs^{-1}$ [34]. The value for the TB induced death rate is based on [25]. The proportion of pulmonary TB cases in Angola is equal to $\nu=0.937$ and the fraction of treatment default and failure for individuals under treatment is equal to $\phi_{T}=0.219$ [36]. We assume that the reinfection factor $\sigma$ in Angola takes the value proposed in [13]. According to WHO, the proportion of detected cases in a year is equal to $k=0.79$ [29]. The rate at which infectious individuals enter treatment is estimated to be $\tau=2.13\,yrs^{-1}$. The values of the parameters $\delta$, $\phi$, $\omega$ and $\delta_{T}$ are taken from [13]. The recruitment rate value $\eta=1287900$ is based on the population projections from Population Reference Bureau [38]. The initial conditions $N(0)$, $\mathcal{S}(0)$, $\mathcal{P}(0)$, $\mathcal{L}(0)$, $\mathcal{I}(0)$, $\mathcal{T}(0)$ are based on data from [23, 35, 37]. All previous values are resumed in Table 1.

If we firstly keep all parameters fixed (see Table 1), we have

[TABLE]

Then we vary one of the parameters $\beta^{A}$, $\beta^{C}$, $\beta^{G}$, $k^{C}$, $\phi^{G}_{T}$, $a^{A}$, $a^{B}$, or $\zeta$ in the ranges

[TABLE]

where $\theta=0.2$. Each simulation gives a curve $x\mapsto R_{0}^{T}(x)$, where $x$ is one of the above parameters, for which we find a best fitting curve in one of the models

[TABLE]

for some constants $a_{0},\dots,a_{n},b_{0},b_{1},b_{2}\in\mathbb{R}^{N}$.

Table 2 shows several curve fittings for the map $x\mapsto R_{0}^{T}(x)$. By “best fitting” we mean a model, chosen between the above models (3.2), where the square root of the sum of squares of the residuals $SQR=\sqrt{\sum_{i}r_{i}^{2}}$ has a minimum value or is smaller than the number of significant digits in determining $R_{0}^{T}$, i.e., $10^{-8}$. The same procedure applied to $R_{0}^{A},R_{0}^{C},R_{0}^{G}$ gave results compatible with the analytic formula (3.1).

3.1.1 Variation of the TB transmission rates

(i.e., changing $\beta^{A}$, $\beta^{C}$ and $\beta^{G}$)

A variation of $20\%$ in the value of $\beta^{A}$ implies a variation of approximately $20\%$ to $R_{0}^{T}$. However, the same variation of $20\%$ in the values of $\beta^{C}$ and $\beta^{G}$ affects $R_{0}^{T}$ less than $1\%$. Contrary to (3.1), the parameters $\beta^{A},\beta^{C},\beta^{G}$ do not appear linearly in the calculation of $R_{0}^{T}$, although locally look similar to an affine function, see Fig. 3.

The variation of $\beta^{A}$ has also a significative impact on the community and the host country, namely, in the number of infected and infectious individuals after $5$ years, see Fig. 4. Defining $I_{X}(t,s)=\left.I_{C}(t)\right|_{\beta^{A}=s}$ with $X\in\{C,G\}$, we have

[TABLE]

An increase (decrease) of $20\%$ in $\beta^{A}$ implies a 5 years increase of approximately $20\%$ (decrease of $30\%$) in $I_{C}$ and $I_{G}$, respectively. This enforce the importance of additional effort to treat TB in countries with high TB incidence, not only because of their population health improvement, but also because of the implications on the health of individuals in other host countries.

3.1.2 Variation in the transfer of individuals (i.e., changing $a^{A}$ and $a^{B}$)

The transfer of individuals between (A) and (C)+(G) (i.e., (B)) is determined by the functions $\gamma_{A}(t)$ and $\gamma_{B}(t)$, which are here assumed to be equal to the parameters $a^{A}$ and $a^{B}$. From Fig. 5 it is clear, as expected, that an increment on the flux of individuals moving from areas of lower TB incidence to areas of higher TB incidence reduces $R_{0}^{T}$ and, on the contrary, an increment in the flux of individuals moving from areas of high TB incidence to areas of lower TB incidence increases $R_{0}^{T}$. Note that $R_{0}^{T}$ grows very fast for smaller values of $a^{A}$ and then tends to stabilize with the flux of persons coming from the high incidence TB area.

An interesting phenomena when varying $a^{A}$ appears in the variable $I_{G}$, i.e., the number of infected individuals in (G) (the community), see Fig. 6. It tells us that it is better for the community to have some moderate exchange of persons with the high incidence TB region. Such behavior and its reverse, after some time, seems to be related to the chosen value of $\zeta$ (discussed in the next subsection). It also imply that a careful study of the seasonality distribution of persons traveling between (A) and (B) may be more relevant for (G) than expected a priori. On the host country viewpoint, such phenomena is not noticed as one can see from the evolution of the total number of infected individuals in the host country, i.e., $I_{C}(t)\,N_{C}(t)+I_{G}\,N_{G}(t)$, see Fig. 6.

3.1.3 About the ratio of individuals that stay in the community

versus spread in the host country (i.e., changing $\zeta$)

In what follows we analyze the impact of the existence of a community of immigrants coming from a high incidence TB area on the host country, the country of origin and in the global situation. Recall that $\zeta$ is the percentage of persons traveling that come/go specifically to (G) versus the complementary (C). Hence, the situation $\zeta=0$ means that all persons traveling between Angola and Portugal all come/go to (C) and none to (G). On the contrary, $\zeta=1$ means that all persons traveling between Angola and Portugal all come/go to (G). From the analysis of Fig. 7 (right), it is clear that the existence of a community of immigrants coming from a high incidence TB area is convenient for the host country in order to better control TB spread. Regarding the point of view of Angola, a change in $\zeta$ is not significative as one can see, in Table 3, that $I_{A}$ is not affected by a change in $\zeta$.

On a global viewpoint, a change in $\zeta$ has a big impact on the reproduction number $R_{0}^{T}$, see Fig. 7 (left), for which the existence of communities turn to be also convenient. In fact, the function attains a minimum value that can be estimated from the approximated fitting by a parabolic function as

[TABLE]

see Table 2. Hence, we may say that the optimal value for $\zeta$ is approximately

[TABLE]

4 Numerical results and discussion

Regarding the sensitivity analysis, we numerically simulated the system (2.8)–(2.11) by considering all parameters fixed except one chosen parameter for which we consider three possible values according with

[TABLE]

where $\theta=0.2$ (i.e., a variation of $\pm 20\%$). The middle levels are the values considered when the parameters are fixed.

Considering that system (2.8)–(2.11) has $15$ relevant state-space variables and we are perturbing $8$ parameters (with $3$ levels), even with overlapping of the levels on the same graphic, such analysis implies the study of $360$ functions aggregated in $120$ graphics. We want to quantify and describe the qualitative behavior and difference between the evolutions, when comparing the different levels. Additionally, a direct visual interpretation of the plots may be biased since the plots are not in the same scale, which may give a quite erroneous filling of disparity between functions when, in fact, the difference may be in a small amount, e.g., see Fig. 8. To deal with such issues, in a precise and normalized way, we considered the following procedure.

Let $F_{Y,P,1}(t),F_{Y,P,2}(t),F_{Y,P,3}(t)$ be the evolution functions associated to one of the state-variables

[TABLE]

and to one of the three variation levels of a parameter $P\in\{\beta^{A},\beta^{C},\beta^{G},k^{C},\phi_{T}^{G},a^{A},a^{B},\zeta\}$. Let $\mathcal{T}>0$ denote the total time of simulation. Define

[TABLE]

for $t\in[0,\mathcal{T}]$ and $\mathcal{L}=\{1,2,3\}$. We divide the analysis of the graphics, like in Fig. 8, in three regions of time: beginning for $t\in\mathcal{B}=[0,\frac{1}{3}\mathcal{T}]$; middle when $t\in\mathcal{M}=[\frac{1}{3}\mathcal{T},\frac{2}{3}\mathcal{T}]$; and end when $t\in\mathcal{E}=[\frac{2}{3}\mathcal{T},\mathcal{T}]$. The time set for the complete graph is denoted by $\mathcal{A}=[0,\mathcal{T}]$. Hence, we define

[TABLE]

It is clear, from the linearity of the integral, that $\xi_{\mathcal{A}}=\xi_{\mathcal{B}}+\xi_{\mathcal{M}}+\xi_{\mathcal{E}}$. To understand what $\xi_{\mathcal{A}}$ measures, consider the hypothetical situation where $F_{Y,P,1}(t)\equiv m+\theta$, $F_{Y,P,2}(t)\equiv m$, and consider $F_{Y,P,3}(t)\equiv m-\theta$ for some $m\in\mathbb{R}$ and $\theta>0$. Then,

[TABLE]

So, although different, $\xi_{\mathcal{A}}$ is somehow similar to the variance over the average, which gives an indication of how much the functions are spread from the average value (between them in each instant of time). The definition of $\xi_{\mathcal{A}}$ is also invariant to scale factors, which is quite useful to eliminate erroneous interpretations of graphics, that may happen without such measuring tools.

For the qualitative description of the variability of the evolution functions, we introduced the following tagging notation based on concrete specifications:

(cases $\mathcal{A}_{--}$, $\mathcal{A}_{+-}$, $\mathcal{A}_{++}$) if $\max(\xi_{\mathcal{B}},\xi_{\mathcal{M}},\xi_{\mathcal{E}})<0.4$;

• 2.

(cases $\mathcal{B}_{--}$, $\mathcal{B}_{+-}$, $\mathcal{B}_{++}$) if $\mathcal{S}\neq\mathcal{A}$ and $\max(\xi_{\mathcal{B}},\xi_{\mathcal{M}},\xi_{\mathcal{E}})=\xi_{\mathcal{B}}$;

• 3.

(case $\mathcal{M}_{--}$, $\mathcal{M}_{+-}$, $\mathcal{M}_{++}$) if $\mathcal{S}\neq\mathcal{A}$ and $\max(\xi_{\mathcal{B}},\xi_{\mathcal{M}},\xi_{\mathcal{E}})=\xi_{\mathcal{M}}$;

• 4.

(cases $\mathcal{E}_{--}$, $\mathcal{E}_{+-}$, $\mathcal{E}_{++}$) if $\mathcal{S}\neq\mathcal{A}$ and $\max(\xi_{\mathcal{B}},\xi_{\mathcal{M}},\xi_{\mathcal{E}})=\xi_{\mathcal{E}}$;

• 5.

(cases $\mathcal{S}_{--}$ with $\mathcal{S}\in\left\{\mathcal{B},\mathcal{M},\mathcal{E},\mathcal{A}\right\}$) if $\xi_{\mathcal{A}}<0.01$;

• 6.

(cases $\mathcal{S}_{+-}$ with $\mathcal{S}\in\left\{\mathcal{B},\mathcal{M},\mathcal{E},\mathcal{A}\right\}$) if it is not $\mathcal{S}_{--}$ and $\xi_{\mathcal{A}}<0.25$;

• 7.

(cases $\mathcal{S}_{+-}$ with $\mathcal{S}\in\left\{\mathcal{B},\mathcal{M},\mathcal{E},\mathcal{A}\right\}$) if it is not $\mathcal{S}_{--}$ and $\mathcal{S}_{+-}$.

If $\xi_{\mathcal{A}}<0.1$, then we consider that the variation is not numerically significative, so it is not discussed. Table 3 resumes the sensitivity analysis, where the only tag behaviors that appear are $\mathcal{B}_{+-}$, $\mathcal{M}_{+-}$, $\mathcal{E}_{+-}$, and $\mathcal{E}_{++}$.

Table 3 is quite explanatory and shows relations between parameter perturbations and epidemiological compartments, in a mathematically precise and rather simple visual representation way. The variation of some parameters just gives the expected behavior, which shows that the proposed model is suitable for the situation under study. On the other hand, it also shows that some parameters that a priori we do not give much attention, as the distribution of persons between (G) and (C) (i.e., $\zeta$), play an important role in TB spread.

5 Conclusions

In this paper, we propose and analyze a new mathematical model for TB transmission that considers internal transfer of individuals. As a case-study, we consider a situation with three populations, namely, Angola (a country with high TB incidence), people living in a semi-closed community of Angola natives, and other persons living in Portugal (a country with low TB incidence). Each of the previous subsystems is divided into five epidemiological categories, which follow the TB transmission dynamics found in [13].

For the analysis and verification of the results presented in this paper, we developed a software tool, so-called sDL [42], that combines in the same framework the power of pre-processing systems (as m4 [12] and cpp [31]), a logical verification tool for classical and hybrid systems (as SMT [41] or KeYmaera [18]), a computer algebra system (as Maple [14]), and a numerical computing language (as Matlab [22]). The pre-processing systems allow the existence of a unique and general file, where constants and ODEs are defined in two hierarchical levels, in order to be used across all tools. The verification tool and the computer algebra system allowed to test the validity of some assumptions and verify the correctness of analytic/algebraic formulae. As expected, the numerical computing language allowed to do the numeric simulations and generate the corresponding graphics. Considering the potential of the software tool sDL, in a forthcoming publication, we intend to study real situations that are modeled by pure hybrid model systems, e.g., transmission coefficients that are discontinuous functions varying with climate and season conditions.

Simulations and sensitivity analysis show that variations of the transmission coefficient on the origin country has a big influence on the number of infected (and infectious) individuals on the community and the host country. This enforce the importance of an additional effort to treat TB and improve health conditions in countries with high TB incidence, since they remarkably affect (in long term) the health of individuals on other countries. As expected, an increment on the flux of individuals moving from areas of lower TB incidence to areas of higher TB incidence reduces the global reproduction number and an increment in the flux of individuals moving from areas of high TB incidence to areas of lower TB incidence increases the global reproduction number, but also introduce modifications in the evolution of each disease category that is not linearly proportional to flux rate. From the community point of view, it is better to have some moderate exchange of persons with the high incidence TB region. Seasonality distribution of persons traveling between Angola and Portugal has an important impact in the number of infected (and infectious) individuals in the community.

The main conclusion is that, contrary to some beliefs, the existence of a community of immigrants coming from a high incidence TB area seems to be convenient in a global point of view, as well as for the host country, in order to better control TB spread. On the other hand, it does not affect the TB incidence in the origin country of the immigrant community. By nonexistence of the community of immigrants we mean the situation where the individuals traveling are spread uniformly on the host country. As shown above, a key parameter in such analysis is the percentage of persons traveling from the high incidence TB area that will stay in the community. Such parameter has an optimal value for TB control, in the sense of minimizing the global reproduction number, that is near to $47\%$. The obtained results are valid under the hypothesis of a semi-closed community. Further studies are necessary for the situation without any flux restrictions.

Acknowledgments. Work partially supported by Portuguese funds through the Center for Research and Development in Mathematics and Applications (CIDMA) and the Portuguese Foundation for Science and Technology (FCT), within project UID/MAT/04106/2013. Rocha is also supported by the FCT project “DALI – Dynamic logics for cyber-physical systems: towards contract based design” with reference P2020-PTDC/EEI-CTP/4836/2014; Silva by the FCT post-doc fellowship SFRH/BPD/72061/2010; Silva and Torres by project TOCCATA, reference PTDC/EEI-AUT/2933/2014, funded by Project 3599 – Promover a Produção Científica e Desenvolvimento Tecnológico e a Constituição de Redes Temáticas (3599-PPCDT) and FEDER funds through COMPETE 2020, Programa Operacional Competitividade e Internacionalização (POCI), and by national funds through FCT. The authors are grateful to two referees for useful comments and suggestions.