An epidemic model with transport-related infection incorporating awareness and screening

Abstract

In this paper, an SWEIQR epidemic model with transport-related infection is proposed. The model considers inter-patch travel with entry-departure screening. The reproduction number, $R_{ed}^\phi$, is computed and analyzed with respect to awareness and screening parameters. The analytic computations show that the disease-free equilibrium in the absence of travel is globally asymptotically stable when $R_\omega \le 1$ and unstable otherwise. The trans-critical bifurcation occurs at $R_\omega =1$ and the locally stable endemic equilibrium point appears if $R_\omega >1$ near to $R_\omega =1$. The numerical simulations are performed to verify the analytical computation and explore the dynamic behavior with respect to different model parameters. The result shows that disseminating awareness through the population reduces the spread of disease. Furthermore, the full model results show that the departure screening may reduce the spread of disease in each patch.

1 Introduction

The movement of the population is a common phenomenon in global human society. Travel for trade, education, or tourism led to more frequent contact with the people. Movement is the major cause for the transmission of several diseases throughout the world. The past global trends show that many sexually transmitted or infectious diseases like Spanish flu, SARS, HIV/AIDS, cholera, Nipah and Ebola are easily transported from one region to the other. The population dispersal is the sole cause of Covid-19 pandemic. The reductions in the population flow, interpersonal contact, and travel bans would have been helpful to reduce transmission if implemented timely [9, 36]. Several measures such as spreading awareness, quarantining/isolation of infected/exposed individuals, screening/testing, and vaccination may control the endemic/pandemic. Various epidemic models have been proposed by considering such measures to predict and prevent the spread of infectious diseases.

Different researchers have studied the effect of disseminating awareness through media on infectious disease dynamics [19, 20, 25]. They considered the media as one dynamic variable, and their findings show that dissemination of awareness programs through media has a possibility to slow down the disease prevalence, but it does not affect the epidemic threshold of the disease due to media alerts. Yang et al. (see, [35]) proposed two different models incorporating awareness. The first model considered awareness programs as one class, and in the second model, they extend it with two distinct susceptible classes (aware and unaware). Their study concludes that the disease transmission rate decreases as the number of awareness programs increase. Furthermore, the second model exhibits the backward bifurcation. Sahu et al. (see, [24]) also indicate that media awareness helps to mitigate the disease burden over time by lowering the level of infection, but it doesn’t affect the effective reproduction number of the model. Agaba et al. (see, [2]) studied an ${\varvec{{SIR}}}$ mathematical model considering the impact of awareness on the dynamics of infectious disease. Their model incorporates dissemination of individual awareness arising from direct contacts between unaware and aware individuals and public awareness released from information campaigns. Misra et al. (see, [21]) have also studied an epidemic model with awareness programs by media. In this model, using optimal control theory, they have identified an optimal implementation rate of awareness campaigns so that disease can be controlled with minimal possible expenditure on awareness campaigns. Kumar et al. proposed an SIRZ epidemic model with information density (Z) in [12]. The study concluded that the effect of information-related behavioral responses played a crucial role in the absence of other controls. The models reviewed above were not considered individual movement between regions.

Commercial globalization and population movements are the main causes of the international spread of microorganisms. The mass movement of large numbers of people creates new opportunities for the spread and establishment of common or novel infectious diseases [28]. The globalization of infectious disease epidemiology will require the corresponding development of integrated programs to anticipate and manage these diseases in response to an increasingly mobile patient population [10]. Many researchers studied the mathematical modeling with population dispersal or transport-related infectious disease between two regions to understand the dynamic behavior of the disease. Mishra et al. (see, [18]) studied the dynamics of secondary dengue infection in two patches, and their study found that immigration in a patch increases the basic reproduction in the respective patch and vice-versa. Wang et al. (see, [31,32,33]) studied an epidemic model with population dispersal between two patches; their study indicates that the disease may spread when the population migrates in between these patches. A rabies transmission model with dog population dispersal is investigated by Liu et al. (see, [14]). In this study, they considered a two patch SEIRS epidemic model to analyze the impact of travel on the spatial spread of dog rabies between two patches. Their finding shows that when border control is properly implemented, it could stop the spread of rabies infection between patches.

The transport system among regions is one of the main factors affecting the outbreak of diseases. Sattenspiel L. and Dietz K. studied a model of spreading infectious disease among areas incorporating mobility process. The SIR epidemic model was formulated and summarized the behavior of a range of mobility from complete isolation to permanent migration between regions. They discuss how mobility and disease transmission interact with one another [26]. In this study, they also focus on the age-related epidemic model with West Indian island of Dominica measles data. However, the threshold of the model was not computed, the stability of equilibrium numbers and numerical simulations were also not done. Fei Xu et al. constructed an SIRS model to study the spread of infection in several locations through a transportation system. They analyzed the spread of disease in identical and non-identical cities. The reproduction number was computed, the stability of equilibrium points and numerical simulations were discussed. Their investigations indicate that the transportation efficiency, improving sanitation, and ventilation of the public transportation system have a positive effect of decreasing the chance of an outbreak occurring [34]. Even though the researchers studied the impact of transportation systems connecting regions on the disease expansion, they did not consider awareness, quarantine, and screening that affect disease dynamics.

Moreover, Some epidemic models with population dispersal have been considered with transport-related infection. Takeuchi et al. (see, [29]) proposed an SIS mathematical model to demonstrate the dynamics of the disease between two cities due to transport-related infection in the population dispersal. Further, their analysis shows that transport-related infection intensifies the disease to spread. Liu et al. (see, [16]) proposed an SIQS epidemic model to study the effect of transport-related infection and screening. Their analysis indicates that screening has the possibility to eradicate the spread of disease led by transport-related infection even when the disease is endemic in both isolated cities. Wan et al. (see, [30]) extended an epidemic model proposed by [29] to an SEIS model to investigate the effect of transport-related infection on the spread and control of infectious disease. Their study shows that even if infected individuals are prohibited from traveling, the movement of exposed individuals can bring infections from regions to each other, and transport-related infection may raise the disease to spread. [4], and [15] proposed the SIR and SIRS epidemic models respectively with transport-related infection; both studies indicate that they established global endemic equilibrium with sufficient conditions and the transport-related infection may lead the disease endemic even if both isolated regions are disease-free. Denphedtnong et al. proposed an SEIRS model with transport-related infection, and concluding that the transport-related infection intensifies the disease [6]. It may lead to an endemic situation in each region. They applied the model to the real data of 2003 SARS outbreak. These models investigate the dynamic behavior of the disease when individual movement is allowed between patches and the impact of transport-related infection on the spread of disease.

Public awareness campaigns influence the transmission of infection in the system. The role of departure and entry screening together with transport-related infection are crucial for the dispersal of endemic diseases. These should be integrated with disease dynamics. To our knowledge, no mathematical studies have attempted to portray the dynamics of a transport-related endemic model that includes awareness, quarantine and screening together. Accordingly, an SWEIQR model is proposed in this paper. The model incorporates the impact of awareness and quarantine together with entry/departure screening on the spread of disease led by transport. To understand the dynamic behavior of the proposed model, this paper is organized as follows: In Sect. 2, the proposed model is formulated and presents the description of parameters. The analysis of the travel and non-travel models is carried out in Sect. 3. The analysis includes computation of reproduction numbers, sensitivity index of parameters, and existence/stability of equilibrium points. In Sect. 4, the numerical simulations are presented, and the graphical results are discussed. Finally, the conclusion is presented in Sect. 5.

2 Model formulation

The model considers the total population N(t) be distributed in n patches. Let the population size in ith patch be $N_i(t)$ such that $N(t)=\sum _{i=1}^nN_i(t)$. The population in each patch is divided into six distinct classes: Susceptibles ($S_i$), aware ($W_i$), exposed ($E_i$), infected ($I_i$), quarantined ($Q_i$) and recovered ($R_i$). The susceptible are those individuals who are not infected but will become infected when they come in contact with infectious individuals. The aware individuals have knowledge about the disease and take precautionary measures as well as follow strict disease-specific appropriate behavior. Exposed individuals are infectious but do not infect others. Infected are those individuals who are infectious and may infect others. Quarantine individuals are those individuals who are isolated from the population due to infectiousness. Recovered are individuals who are rescued from infectiousness either with treatment or by their immunity. Accordingly,

$$\begin{aligned} N_i(t)=S_i(t)+W_i(t)+E_i(t)+I_i(t)+Q_i(t)+R_i(t),\qquad (i=1,2, 3, . . ., n). \end{aligned}$$

The model considers the following basic assumptions:

• (i) It is assumed that the disciplined behavior (awareness) provides them perfect protection and will not get an infection. The transition from a susceptible class to a non-infected aware class is considered as $W_i$. Once an aware individual becomes reluctant to follow the specified norms, he loses protection against the disease and again become susceptible.

• (ii) Disease is transmitted with the incidence rate ($\beta _i S_iI_i/N_i, \quad i = 1,2,3,...,n$ within the patch i).

• (iii) The movement of individuals between the n-patches are restricted. The susceptible, aware, exposed, and recovered individuals are allowed to travel out from ith patch proportional to their respective sizes with constant rate $\alpha _i$. These traveling individuals from ith patch will move to jth patch at a rate $\alpha _{ji}$ such that $\alpha _i=\sum _{j=1}^{n}\alpha _{ji}$,    (for $i=1,2,3...,n$, $j\ne i$).

• (iv) The travelers are screened at the time of entry/departure from ith patch. Let $e_i$ be the probability of detecting infectivity at entry time in the ith patch and $d_i$ be the probability of detecting infectivity at departure.

• (v) A fraction $\theta _{ji}$ ( $\alpha _{ji} > \theta _{ji}$) of infected individuals travel out of the ith patch to jth patch. They may some how evade the screening or the test may be false negative. Let the fraction $d_i$ of infective travelers are stopped at the exit point. Therefore, $(1-d_i )\theta _{ji}$ fraction of infected individuals of ith patch travel to jth patch. Let us denote $(1-d_i)\theta _i$ as $\sum _{j=1}^{n}(1-d_i)\theta _{ji}$,   (for $i=1,2,3...,n$, $i\ne j$) be the infective individuals who traveled out of ith patch.

• (vi) For individuals in patch i travel to patch $j \quad (i\ne j)$, disease is transmitted with the incidence rate $(1-d_i)\phi \alpha _{ji}S_i\theta _{ji} I_i/N_i$    $(i,\,j=1,2,3,...,n$    and $i\ne j)$, where $\phi$ is the infection transmission rate due to transport and $1-d_i$ is the fraction of non-detected infected individuals at departure time and going out from patch i.

• (vii) At the time of travel there is no birth and death and the recovered individuals do not lose immunity.

• (viii) Quarantined individuals are not allowed to travel and they do not contact with anyone.

• (ix) The model parameters such as the recruitment rate ($\varLambda _i$), death rates ($\mu _i$, $\delta _i$ and $\kappa _i$) disease transmission rate ($\beta _i$), awareness rate ($\omega _i$), recovery rates ($\gamma _i$ and $\sigma _i$), the transition rate from exposed to infectious class ($\rho _i$), loss of awareness rate ($\xi _i$) and quarantine rate ($\epsilon _i$) in the ith patch, ($i=1, 2, 3, ... ,n$) are non negative constants.

Based on the assumptions and description of parameters, the following deterministic model of non linear system of ordinary differential equations is formulated.

$$\begin{aligned} \frac{dS_i}{dt}= & {} \varLambda _i -\frac{\beta _i S_iI_i}{N_i}-(\omega _i+\mu _i)S_i+\xi _iW_i-\alpha _iS_i\nonumber \\&+\sum _{j=1,i\ne j}^{n}\left[ \alpha _{ij} S_j -\frac{(1-d_j)\phi \alpha _{ij}\theta _{ij}S_jI_j}{N_j}\right] \end{aligned}$$

(2.1)

$$\begin{aligned} \frac{dW_i}{dt}= & {} \omega _i S_i -(\xi _i + \mu _i)W_i -\alpha _iW_i+\sum _{j=1,i\ne j}^{n}\alpha _{ij}W_j \end{aligned}$$

(2.2)

$$\begin{aligned} \frac{dE_i}{dt}= & {} \frac{\beta _i S_iI_i}{N_i}-(\rho _i + \mu _i)E_i-\alpha _iE_i\nonumber \\&+\sum _{j=1,i\ne j}^{n}\alpha _{ij} E_j+\phi \sum _{j=1,i\ne j}^{n}\frac{(1-d_j)\alpha _{ij}\theta _{ij}S_jI_j}{N_j} \end{aligned}$$

(2.3)

$$\begin{aligned} \frac{dI_i}{dt}= & {} \rho _iE_i -(\epsilon _i +\mu _i +\delta _i +\gamma _i +\theta _i)I_i+(1-e_i)\sum _{j=1,i\ne j}^{n}(1-d_j)\theta _{ij} I_j \end{aligned}$$

(2.4)

$$\begin{aligned} \frac{dQ_i}{dt}= & {} \epsilon _iI_i-(\sigma _i + \mu _i + \kappa _i)Q_i+d_i\theta _i I_i+e_i\sum _{j=1,i\ne j}^{n}(1-d_j)\theta _{ij} I_j \end{aligned}$$

(2.5)

$$\begin{aligned} \frac{dR_i}{dt}= & {} \gamma _iI_i+\sigma _iQ_i -\mu _iR_i-\alpha _iR_i+\sum _{j=1,i\ne j}^{n}\alpha _{ij} R_j \end{aligned}$$

(2.6)

The following initial conditions are associated with the system:

$$\begin{aligned}&S_i(0)> 0,\quad W_i(0)\ge 0, \quad E_i(0)\ge 0, \quad I_i(0)\ge 0, \end{aligned}$$

(2.7)

$$\begin{aligned}&Q_i(0)\ge 0,\quad R_i(0)\ge 0,\quad (i=1,2,3,...,n). \end{aligned}$$

(2.8)

Table 1 Description of parameters used in the model

The first four equations are independent of the last two components, so for computational ease, we consider the sub-system model equation (2.1)–(2.4) with initial conditions given in (2.7). The total population at time t in the ith patch of the sub-system is ${\mathcal {M}}_i(t)$, such that

$$\begin{aligned} {\mathcal {M}}_i(t)=S_i(t)+W_i(t)+E_i(t)+I_i(t)\le N_i(t), \quad i=1,2,3,...,n. \end{aligned}$$

2.1 Basic properties of the model

The dynamics of the sub-system will be studied in the biologically feasible closed set defined by:

$$\begin{aligned} \varOmega =\left\{ (S_i,W_i,E_i,I_i)\in R_+^{4n}:\sum _{i=1}^{n}(S_i+W_i+E_i+I_i)=\sum _{i=1}^{n}{\mathcal {M}}_i\le \frac{\sum _{i=1}^{n}\varLambda _i}{\mu }\right\} \nonumber \\ \end{aligned}$$

(2.9)

Lemma 1

The region $\varOmega$ given in (2.9) with initial conditions given in (2.7) is positively invariant with respect to the sub-system model (2.1)–(2.4).

Proof

Adding the differential equations (2.1) - (2.4) in the model provides the derivative of total population of each patch given by:

$$ \begin{aligned}&\frac{d{\mathcal {M}}_i(t)}{dt}=\varLambda _i-\mu _i {\mathcal {M}}_i(t)-(\epsilon _i+\gamma _i+\delta _i) I_i(t)-\alpha _i S_i+\sum _{j=1}^{n}\alpha _{ij}S_j -\alpha _i W_i \\&+\sum _{j=1}^{n}\alpha _{ij}W_j-\alpha _i E_i+\sum _{j=1}^{n}\alpha _{ij}E_j \\&-\,\,\theta _i I_i+(1-e_i)\sum _{j=1}^{n}(1-d_j)\theta _{ij}I_j,\quad (for \quad i=1,2,3,...,n\quad \& \quad i\ne j) \end{aligned}$$

The grand total population is ${\mathcal {M}}(t)=\sum _{i=1}^{n}{\mathcal {M}}_i(t),\quad i=1,2,3,...,n$ and its rate of change becomes:

$$\begin{aligned}&\frac{d{\mathcal {M}}(t)}{dt}=\sum _{i=1}^{n}\left[ \varLambda _i-\mu _i {\mathcal {M}}_i(t)-(\epsilon _i+\gamma _i+\delta _i)I_i(t)-e_i\sum _{j=1}^{n}(1-d_j)\theta _{ij}I_j\right] \\&-\sum _{i=1}^{n}\alpha _i S_i+\sum _{i=1}^{n}\sum _{j=1}^{n}\alpha _{ij}S_j-\sum _{i=1}^{n}\alpha _iW_i\\&+\sum _{i=1}^{n}\sum _{j=1}^{n}\alpha _{ij}W_j-\sum _{i=1}^{n}\alpha _iE_i+\sum _{i=1}^{n}\sum _{j=1}^{n}\alpha _{ij}E_j \\&-\sum _{i=1}^{n}\theta _iI_i+\sum _{i=1}^{n}\sum _{j=1}^{n}(1-d_j)\theta _{ij}I_j,\quad (for \quad i\ne j) \\&\Longrightarrow \frac{d{\mathcal {M}}(t)}{dt}=\sum _{i=1}^{n}\left[ \varLambda _i-\mu _i {\mathcal {M}}_i(t)-(\epsilon _i+\gamma _i+\delta _i)I_i(t)-e_i\sum _{j=1}^{n}(1-d_j)\theta _{ij}I_j\right] \\&-\sum _{i=1}^{n}\alpha _i S_i+\sum _{i=1}^{n}\sum _{j=1}^{n}\alpha _{ji}S_i-\sum _{i=1}^{n}\alpha _iW_i \\&+\sum _{i=1}^{n}\sum _{j=1}^{n}\alpha _{ji}W_i-\sum _{i=1}^{n}\alpha _iE_i+\sum _{i=1}^{n}\sum _{j=1}^{n}\alpha _{ji}E_i \\&-\sum _{i=1}^{n}\theta _iI_i+\sum _{i=1}^{n}\sum _{j=1}^{n}(1-d_i)\theta _{ji}I_i,\quad (for \quad i\ne j) \end{aligned}$$

From the basic assumptions,

$$\begin{aligned} \alpha _i =\sum _{j=1}^{n}\alpha _{ji} \text {and} (1-d_i)\theta _i =\sum _{j=1}^{n}(1-d_i)\theta _{ji} (\text {for} i\ne j) \end{aligned}$$

since $(1-d_i)\theta _i\le \theta _i$, we have:

$$\begin{aligned} \frac{d{\mathcal {M}}(t)}{dt}\le & {} \sum _{i=1}^{n}\left[ \varLambda _i-\mu _i {\mathcal {M}}_i(t)-(\epsilon _i+\gamma _i+\delta _i)I_i(t)-e_i\sum _{j=1}^{n}(1-d_j)\theta _{ij}I_j\right] \\\le & {} \sum _{i=1}^{n}\varLambda _i-\mu \sum _{i=1}^{n} {\mathcal {M}}_i(t)= \sum _{i=1}^{n}\varLambda _i-\mu {\mathcal {M}}(t) ,\quad where \quad \mu =\min _{1\le i\le n}\{\mu _i\} \end{aligned}$$

It follows that

$$\begin{aligned} {\mathcal {M}}(t)\le \frac{\sum _{i=1}^{n}\varLambda _i}{\mu }-\left( \frac{\sum _{i=1}^{n}\varLambda _i}{\mu }-{\mathcal {M}}(0)\right) e^{-\mu t} \end{aligned}$$

Where, ${\mathcal {M}}(0)$ is the initial grand total population. If ${\mathcal {M}}(0)\ge \sum _{i=1}^{n}\varLambda _i/\mu$, then either the solution enters $\varOmega$ in finite time, or ${\mathcal {M}}(t)$ approaches to $\sum _{i=1}^{n}\varLambda _i/\mu$ and ${\mathcal {M}}(t)\le \sum _{i=1}^{n}\varLambda _i/\mu$ as $t\rightarrow \infty$ if ${\mathcal {M}}(0)\le \sum _{i=1}^{n}\varLambda _i/\mu$, which implies the solution of the system is bounded. Thus, the region $\varOmega$ is positively invariant for the system given in (2.1)–(2.4).

Theorem 1

Solutions of the state variables in the model (2.1)–(2.4) with non-negative initial conditions remain non-negative for all time $t\ge 0$.

Proof

For the given initial conditions of the model, it can be shown that all solutions of the system (2.1)–(2.4) are positive for all time $t\ge 0$. By assuming a contradiction [11, 17],

there exists a time $t_{i1}$, for $i=1,2,3,...,n$ such that

$$\begin{aligned} S_i(t_{i1})=0,\quad S_i(t)\ge 0,\quad W_i(t)\ge 0,\quad E_i(t)\ge 0,\quad I_i(t)\ge 0, \quad 0\le t\le t_{i1}\nonumber \\ \end{aligned}$$

(2.10)

there exists a $t_{i2}$ such that

$$\begin{aligned} W_i(t_{i2})=0,\quad S_i(t)\ge 0,\quad W_i(t)\ge 0,\quad E_i(t)\ge 0,\quad I_i(t)\ge 0, \quad 0\le t\le t_{i2}\nonumber \\ \end{aligned}$$

(2.11)

there exists a $t_{i3}$ such that

$$\begin{aligned} E_i(t_{i3})=0,\quad S_i(t)\ge 0,\quad W_i(t)\ge 0,\quad E_i(t)\ge 0,\quad I_i(t)\ge 0, \quad 0\le t\le t_{i3}\nonumber \\ \end{aligned}$$

(2.12)

there exists a $t_{i4}$ such that

$$\begin{aligned} I_i(t_{i4})=0,\quad S_i(t)\ge 0,\quad W_i(t)\ge 0,\quad E_i(t)\ge 0,\quad I_i(t)\ge 0, \quad 0\le t\le t_{i4}\nonumber \\ \end{aligned}$$

(2.13)

Considering $(1-d_j)\phi \theta _{ij} <1$, ($i,j=1,2,3,...,n\quad \& \quad i\ne j$), from equation (2.1) in the model system and the case in (2.10), we have:

$$ \begin{aligned}&\frac{dS_i(t_{i1})}{dt}=\varLambda _i -\frac{\beta _i S_i(t_{i1})I_i(t_{i1})}{N_i(t_{i1})}-(\omega _i +\mu _i)S_i(t_{i1})+\xi _i W_i(t_{i1})- \alpha _iS_i(t_{i1})\\&\qquad +\sum _{j=1,i\ne j}^{n}\left[ \alpha _{ij} S_j(t_{i1})-\frac{\alpha _{ij}(1-d_j)\phi \theta _{ij}S_j(t_{i1})I_j(t_{i1})}{N_j(t_{i1})}\right] \\&\quad =\varLambda _i+\xi _i W_{i}(t_{i1})+\sum _{j=1}^{n}\left[ \alpha _{ij} S_j(t_{i1})-\frac{\alpha _{ij}(1-d_j)\phi \theta _{ij}S_j(t_{i1})I_j(t_{i1})}{N_j(t_{i1})}\right] \ge 0,\\&\quad \quad i=1,2,3,...,n \& i\ne j \end{aligned}$$

It leads to have $S_{i}(t)<0$ for $t\ge 0$, which contradicts the assumption that $S_{i}(t)\ge 0$ ($i=1,2,3,...,n$) for $t\ge 0$.

From the second equation (2.2) in the model system and the case in (2.11), we have:

$$ \begin{aligned} \frac{dW_i(t_{i2})}{dt}= & {} \omega _iS_{i}(t_{i2}) -(\xi _i+\mu _i)W_{i}(t_{i2})-\alpha _iW_{i}(t_{i2})+\sum _{j=1,i\ne j}^{n}\alpha _{ij}W_j(t_{i2}) \\= & {} \omega _i S_{i}(t_{i2})+\sum _{j=1}^{n}\alpha _{ij}W_j(t_{i2})\ge 0,\quad i=1, 2, 3, ... , n \quad \& \quad i\ne j \end{aligned}$$

It means that $W_{i}(t)<0$ for $t\ge 0$, which is a contradiction that $W_{i}(t)\ge 0$ ($i=1,2,3,...,n$) for $t\ge 0$ from our assumption. Similarly, it can be shown that $E_{i}(t)\ge 0$, $I_{i}(t)\ge 0$ ($i=1,2,3,...,n$) for all $t\ge 0$.

Thus, for the given initial conditions in the domain $\varOmega$ the solutions $S_{i}(t)$, $W_{i}(t)$, $E_{i}(t)$ and $I_{i}(t)$ (for $i=1,2,3,...,n$) remain non-negative for all time $t\ge 0$.

3 Model analysis

In this section, the dynamical system without inter-patch migration is discussed first. The full model with inter-patch migration is discussed next with $n=2$. The analytical computation with different parameters of each patch is complicated. So, we assume that the rate of biological parameters (say ${\mathcal {P}}_i$) is the same in the ith patch for simplicity (i.e ${\mathcal {P}}_i$=${\mathcal {P}}$).

3.1 Model without migration

In the absence of travel (migration), $\alpha =0$, $\theta =0$ and considering the same parameters in each patch, the model equations (2.1)–(2.4) are written as:

$$\begin{aligned} \frac{dS_{i}}{dt}= & {} \varLambda -\frac{\beta S_{i}I_{i}}{N_i} -(\omega +\mu )S_{i}+\xi W_{i} \end{aligned}$$

(3.1)

$$\begin{aligned} \frac{dW_{i}}{dt}= & {} \omega S_{i} -(\xi + \mu )W_{i} \end{aligned}$$

(3.2)

$$\begin{aligned} \frac{dE_{i}}{dt}= & {} \frac{\beta S_{i}I_{i}}{N_i}-(\rho + \mu )E_{i} \end{aligned}$$

(3.3)

$$\begin{aligned} \frac{dI_{i}}{dt}= & {} \rho E_{i} -(\epsilon +\mu +\delta +\gamma )I_{i} \end{aligned}$$

(3.4)

With associated initial conditions:

$$\begin{aligned} S_i(0)> 0, \quad W_i(0)\ge 0, \quad E_i(0)\ge 0, \quad I_i(0)\ge 0, \end{aligned}$$

(3.5)

From Theorem 1, it is easy to establish that the positive invariant region for the system (3.1)–(3.4) is

$$\begin{aligned} \varOmega ^*=\left\{ (S_{i},W_{i},E_{i},I_{i})\in R_{+}^{4}:{\mathcal {M}}_i\le \frac{\varLambda }{\mu }\right\} \end{aligned}$$

(3.6)

3.1.1 Basic reproduction number

The system (3.1)–(3.4) has a disease free equilibrium point $E^0$ given by:

$$\begin{aligned} E^{0}=\left( \frac{\varLambda (\xi +\mu )}{\mu (\xi +\mu +\omega )},\,\frac{\varLambda \omega }{\mu (\xi +\mu +\omega )},\,0,\,0\right) \end{aligned}$$

(3.7)

The reproduction number $R_\omega$ of the system can be computed by the next generation matrix [8] as follows:

In the system (3.1)–(3.4) $X_{i}=\left( E_{i},I_{i}\right) ^T$ and $Y_{i}=\left( S_{i},W_{i}\right) ^T$ are classified as diseased and non-diseased compartments of ith patch respectively. The diseased compartments of the system can be written as:

$$\begin{aligned} \frac{dX_{i}}{dt}={\mathcal {F}}(X_{i})-{\mathcal {V}}(X_{i}) \end{aligned}$$

Where,

$$\begin{aligned} {\mathcal {F}}(X_{i})=\left( \begin{array}{c} \frac{\beta I_{i}S_{i}}{N_i}\\ \\ 0 \end{array}\right) \quad {\mathcal {V}}(X_{i})=\left( \begin{array}{c} (\rho +\mu )E_{i}\\ \\ -\rho E_{i}+ (\epsilon +\mu +\delta +\gamma )I_{i} \end{array}\right) \end{aligned}$$

At disease free equilibrium $E^0$ in (3.7) the Jacobian matrix of ${\mathcal {F}}(X_i)$ and ${\mathcal {V}}(X_i)$ are respectively computed as:

$$\begin{aligned} F=\left( \begin{array}{cc} 0&{}\frac{\beta (\xi +\mu )}{(\xi +\mu +\omega )}\\ \\ 0 &{} 0 \end{array} \right) V=\left( \begin{array}{cc} \rho +\mu &{} 0\\ \\ rho &{} \epsilon +\delta +\mu +\gamma \end{array}\right) \end{aligned}$$

Thus, the reproduction number related to awareness $R_\omega$ of the model in each patch is the spectral radius $\rho (FV^{-1})$, which is given by:

$$\begin{aligned} R_\omega =\frac{\beta \rho (\xi +\mu )}{(\rho +\mu )(\xi +\mu +\omega )(\epsilon +\delta +\mu +\gamma )} \end{aligned}$$

(3.8)

and in the absence of awareness ($\omega$=0), $R_\omega$ is going to be the basic reproduction number $R_0$, which is:

$$\begin{aligned} R_0=\frac{\beta \rho }{(\rho +\mu )(\epsilon +\delta +\mu +\gamma )} \end{aligned}$$

(3.9)

From (3.8) and (3.9) we observe that $R_\omega \le R_0$. It follows, disseminating awareness through the population may reduce the infection to disease free state of each patch.

3.1.2 Existence of endemic equilibrium

Let the system in the model (3.1)–(3.4) has an equilibrium point $\Gamma ^{*}=(S^{*}_{i},W^{*}_{i},E^{*}_{i},I^{*}_{i})$. Set the system equations equal to zero. The solution gives the following equilibrium values interms of the incidence rate $\lambda ^*_{i}$,

$$\begin{aligned} \begin{aligned}&S^{*}_{i}=\frac{\varLambda (\xi +\mu )}{(\lambda ^*_i+\mu )(\xi +\mu )+\mu \omega },E^{*}_{i}=\frac{\lambda ^*_iS^*_i}{(\rho +\mu )},\\&W^{*}_{i}=\frac{\varLambda \omega }{(\lambda ^*_i+\mu )(\xi +\mu )+\mu \omega },I^{*}_{i}=\frac{\lambda ^*_iS^*_i\rho }{(\gamma +\delta +\epsilon +\mu )(\mu +\rho )} \end{aligned} \end{aligned}$$

(3.10)

The incidence rate $\lambda ^*_i$ is given by:

$$\begin{aligned} \lambda ^*_i=\frac{\beta I^*_i}{N^*_i}, \quad N^*_i=S^*_i+W^*_i+E^*_i+I^*_i+Q^*_i+R^*_i \end{aligned}$$

(3.11)

From (3.11), the following quadratic equation is obtained:

$$\begin{aligned} a(\lambda ^{*}_{i})^{2}+b\lambda ^{*}_{i}=0 \quad \Rightarrow \quad \lambda ^{*}_{i}=0 \quad \text {or} \quad \lambda ^{*}_{i}=\frac{-b}{a} \end{aligned}$$

Where,

$$\begin{aligned}&a=(\xi +\mu )[\mu (\epsilon +\mu +\delta )(\kappa +\mu )+\mu (\epsilon +\kappa +\mu )\rho +(\delta \mu +(\epsilon +\mu )(\mu +\rho ))\sigma \\&\quad +\gamma (\mu +\rho )(\kappa +\mu +\sigma )] \end{aligned}$$

and

$$\begin{aligned}&b=\chi (1-R_\omega )\\&\Rightarrow \quad \lambda ^{*}_{i}=\frac{\chi }{a}(R_\omega -1) \end{aligned}$$

where,

$$\begin{aligned} \chi =\mu (\kappa +\mu +\sigma )(\xi +\mu +\omega )(\rho +\mu )(\epsilon +\gamma +\mu +\delta ) \end{aligned}$$

It follows that, for $R_\omega <1$ the system (3.1)–(3.4) has one equilibrium point which is disease free $E^0$ defined in (3.7). But it has a unique endemic equilibrium point $\Gamma ^*$ whenever $R_\omega >1$. Substituting $\lambda ^{*}_{i}$ and simplifying, yields the following lemma.

Lemma 2

The system (3.1)–(3.4) has a unique endemic equilibrium point $\Gamma ^{*}=(S^{*}_{i},W^{*}_{i},E^{*}_{i},I^{*}_{i})$ provided that $R_\omega >1$.

$$\begin{aligned} S^{*}_{i}= & {} \frac{a\varLambda (\xi +\mu )}{a\mu \omega +(\xi +\mu )(a\mu +\chi (R_\omega -1))}, \\ \\ W^{*}_{i}= & {} \frac{a\varLambda \omega }{a\mu \omega +(\xi +\mu )(a\mu +\chi (R_\omega -1))}, \\ \\ E^{*}_{i}= & {} \frac{\varLambda (\xi +\mu )\chi (R_\omega -1)}{(\mu +\rho )(a\mu \omega +(\xi +\mu )(a\mu +\chi (R_\omega -1)))},\\ \\ I^{*}_{i}= & {} \frac{\varLambda \rho (\xi +\mu )\chi (R_\omega -1)}{(\gamma +\delta +\epsilon +\mu )(\mu +\rho )(a\mu \omega +(\xi +\mu )(a\mu +\chi (R_\omega -1))} \end{aligned}$$

3.1.3 Stability of equilibrium points

In this section, the local stability of disease free and endemic equilibrium points are discussed.

Theorem 2

(Local stability of DFE) For the system (3.1)–(3.4) in the non travel model, the disease free equilibrium point $E^{0}$ is locally asymptotically stable if $R_\omega <1$. It is unstable if $R_\omega >1$.

The proof is presented in Appendix A.

At disease free point, the Jacobian matrix (see Appendix A.1) has a simple zero eigenvalue and other eigenvalues are negative for $R_\omega =1$. Following this, the disease free equilibrium point of the system (3.1)–(3.4) is non-hyperbolic. Linearization do not show the stability of non-hyperbolic equilibrium points, so it can be analyzed using central manifold theory [1, 3]; it is particularly presented in Appendix B. The computation in Appendix C indicates that the system exhibits a forward bifurcation at $\beta =\beta ^*$ for $R_\omega =1$. Accordingly the disease free equilibrium point changes its stability from locally stable for $R_\omega <1$ to unstable for $R_\omega >1$. Consequently, the following result is established.

Theorem 3

(Local stability of endemic equilibrium) The system (3.1)–(3.4) exhibits a forward bifurcation at $\beta =\beta ^*$ (for $R_\omega =1$) and a locally stable positive endemic equilibrium will appear whenever $R_\omega >1$ near to $R_\omega =1$.

Theorem 4

(Global stability of DFE) The disease-free equilibrium point $E^0$ of the non travel model system (3.1)–(3.4) is globally-asymptotically stable (GAS) if $R_\omega \le 1$.

Proof

To proof the global stability of disease free equilibrium point (3.7), we consider a positive definite function $V_i(t)$ in $\varOmega ^*$. For a positive value a, $V_i(t)$ in the ith patch is defined by:

$$\begin{aligned} V_i(t)=aE_i(t)+I_i(t). \end{aligned}$$

(3.12)

$dV_i(t)/dt\le 0$, when $R_0\le 1$, $R_0$ is given in (3.9). Since $R_\omega \le R_0$, By Lasalle’s invariance principle [13] the disease free equilibrium point is globally asymptotically stable if $R_\omega \le 1$ (see the details of the proof in Appendix D).

3.2 Full model with migration

In this section, full model given in (2.1) - (2.6) with inter-patch travel for $n=2$ is discussed to study the effect of transport-related infection and entry-departure screening. The infectious individuals do not move at the same rate as non-infectious individuals. The quarantine and recovered classes do not appear in the first four equations of the model system (2.1)–(2.6), then the full system is reduced to the following sub-system:

$$\begin{aligned} \frac{dS_i}{dt}= & {} \varLambda _i -\frac{\beta _i S_iI_i}{N_i}-(\omega _i+\mu _i)S_i+\xi _iW_i-\alpha _iS_i\nonumber \\&+\sum _{j=1,i\ne j}^{n}\left[ \alpha _{ij} S_j -\frac{(1-d_j)\phi \alpha _{ij}\theta _{ij}S_jI_j}{N_j}\right] \end{aligned}$$

(3.13)

$$\begin{aligned} \frac{dW_i}{dt}= & {} \omega _i S_i -(\xi _i + \mu _i)W_i -\alpha _iW_i+\sum _{j=1,i\ne j}^{n}\alpha _{ij}W_j \end{aligned}$$

(3.14)

$$\begin{aligned} \frac{dE_i}{dt}= & {} \frac{\beta _i S_iI_i}{N_i}-(\rho _i + \mu _i)E_i-\alpha _iE_i+\sum _{j=1,i\ne j}^{n}\alpha _{ij} E_j\nonumber \\&+\phi \sum _{j=1,i\ne j}^{n}\frac{(1-d_j)\alpha _{ij}\theta _{ij}S_jI_j}{N_j} \end{aligned}$$

(3.15)

$$\begin{aligned} \frac{dI_i}{dt}= & {} \rho _iE_i -(\epsilon _i +\mu _i +\delta _i +\gamma _i +\theta _i)I_i+(1-e_i)\sum _{j=1,i\ne j}^{n}(1-d_j)\theta _{ij} I_j \end{aligned}$$

(3.16)

With associated initial conditions:

$$\begin{aligned} S_i(0)> 0, W_i(0)\ge 0, E_i(0)\ge 0, I_i(0)\ge 0, (i=1,2). \end{aligned}$$

(3.17)

The model has a disease free equilibrium point ($E_\phi ^0$) in the closed region $\varOmega$ defined in (2.9) and since the parameters ${\mathcal {P}}_i$=${\mathcal {P}}$ for $i=1,2$, it is given by:

$$\begin{aligned} E_\phi ^0=\left( \frac{\varLambda (\xi +\mu )}{\mu (\xi +\mu +\omega )},0,0,0,\frac{\varLambda (\xi +\mu )}{\mu (\xi +\mu +\omega )},0,0,0\right) \end{aligned}$$

(3.18)

3.2.1 The reproduction number of transport-related infection with screening

The reproduction number of transport-related infection with screening is computed by next generation matrix approach. The matrix F and V of the reduced sub-system (3.13) - (3.16) at the disease free equilibrium point $E_\phi ^0$ in (3.18) are given by:

$$\begin{aligned} F=\left( \begin{array}{cccc} 0 &{} \frac{\beta (\xi +\mu )}{(\xi +\mu +\omega )} &{} 0 &{} \frac{\phi \alpha \theta (1-d)(\xi +\mu )}{(\xi +\mu +\omega )}\\ 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} \frac{\phi \alpha \theta (1-d)(\xi +\mu )}{(\xi +\mu +\omega )} &{} 0 &{} \frac{\beta (\xi +\mu )}{(\xi +\mu +\omega )} \\ 0 &{} 0 &{} 0 &{} 0 \end{array}\right) \end{aligned}$$

(3.19)

and

$$\begin{aligned} V=\left( \begin{array}{cccc} \alpha +\mu +\rho &{} 0 &{} -\alpha &{} 0 \\ rho &{} \epsilon +\gamma +\mu +\delta +\theta &{} 0 &{} -(1-e)(1-d)\theta \\ alpha &{} 0 &{} \alpha +\mu +\rho &{} 0 \\ 0 &{} -(1-e)(1-d)\theta &{} -\rho &{} \epsilon +\gamma +\mu +\delta +\theta \end{array}\right) \end{aligned}$$

(3.20)

Now the screening related reproduction number $R_{ed}^\phi$ is given by:

$$\begin{aligned} R_{ed}^\phi =\frac{(\mu +\xi )\rho (\beta +(1-d)\alpha \theta \phi )}{(\mu +\xi +\omega )(\mu +\rho )(\epsilon +\gamma +\mu +\delta +(d+e(1-d))\theta )} \end{aligned}$$

(3.21)

In the absence of entry and departure screenings (i.e, when $d=0$ and $e=0$) equation (3.21) becomes the transport-related infection reproduction number given by:

$$\begin{aligned} R^\phi _0=\frac{(\mu +\xi )\rho (\beta +\alpha \theta \phi )}{(\mu +\xi +\omega )(\mu +\rho )(\epsilon +\gamma +\mu +\delta )} \end{aligned}$$

(3.22)

Let’s denote $R^\phi$ and $R_T$ as:

$$\begin{aligned}&R^\phi =R_\omega \frac{\epsilon +\gamma +\mu +\delta }{\epsilon +\gamma +\mu +\delta +(d+e(1-d))\theta }, \quad R_\omega \text { as given in } (3.8)\\&R_T=\frac{(\mu +\xi )\rho (1-d)\alpha \theta \phi }{(\mu +\xi +\omega )(\mu +\rho )(\epsilon +\gamma +\mu +\delta +(d+e(1-d))\theta )} \end{aligned}$$

It follows that, $R_{ed}^\phi =R^\phi +R_T$.

Note that, $R_T=0$ if the transport-related infection is zero (i.e $\phi =0$). Further, it is observed that

$$\begin{aligned} \frac{\epsilon +\gamma +\mu +\delta }{\epsilon +\gamma +\mu +\delta +(d+e(1-d))\theta }<1. \end{aligned}$$

Accordingly, $R_{ed}^{\phi }< R_\omega$. It means that, with proper entry and exit screening, travel may be allowed in the absence of travel-related infection ($\phi =0$) and it may not have adverse effect on disease dynamics.

When $\phi \ne 0$, then $R_\omega \le R^\phi _0$ and $R_{ed}^\phi \le R^\phi _0$. It means that the proper screening reduces the risk of further spread of infection.

3.2.2 Stability analysis for a full infection model

The reduced sub-system model (3.13)–(3.16) has a unique positive endemic equilibrium point $\Gamma _{\phi }^*$=$(S_i^*,W_i^*,E_i^*,I_i^*)$ for $i=1,2$ provided that $R_{ed}^\phi >1$.

$$\begin{aligned} S_i^*= & {} \frac{\varLambda (\mu +\xi )}{\mu (\mu +\xi +\omega )+(\mu +\xi )\lambda ^*_{1i}},\quad E_i^*=\frac{ S_i^*\lambda ^*_{1i}}{\mu +\rho }, \\ \\ W_i^*= & {} \frac{\varLambda \omega }{\mu (\mu +\xi +\omega )+(\mu +\xi )\lambda ^*_{1i}}, \quad I_i^*\nonumber \\= & {} \frac{\rho S^*_i\lambda ^*_{1i}}{(\epsilon +\gamma +\delta +\mu +(d+e(1-d))\theta )(\mu +\rho )}. \end{aligned}$$

Note that,

$$\begin{aligned} \lambda ^*_{1i}=\frac{\mu (\kappa +\mu +\sigma )\chi _1(R^\phi _{ed}-1)}{(\mu +\xi )(\chi _2+\chi _3)} \end{aligned}$$

(3.23)

$$\begin{aligned} \chi _1= & {} (\mu +\rho )(\mu +\xi +\omega )(\epsilon +\gamma +\mu +\delta +(d+e(1-d))\theta )\\ \\ \chi _2= & {} \mu (\delta +\mu +\epsilon +(d+e(1-d)) \theta )(\kappa +\mu )+\mu (\epsilon +\kappa +\mu +(d+e(1-d))\theta )\rho \\ \\ \chi _3= & {} (\delta \mu +(\epsilon +\mu +(d+e(1-d))\theta )(\mu +\rho ))\sigma +\gamma (\mu +\rho )(\kappa +\mu +\sigma ) \end{aligned}$$

The equilibrium densities at the endemic point $\Gamma _\phi ^*$ are computed as

$$\begin{aligned} S_i^*= & {} \frac{\varLambda (\mu +\xi )(\chi _2+\chi _3)}{\mu (\mu +\xi +\omega )(\chi _2+\chi _3)+\mu (\kappa +\mu +\sigma )\chi _1(R_{ed}^\phi -1)}, \\ W_i^*= & {} \frac{\varLambda \omega (\chi _2+\chi _3)}{\mu (\mu +\xi +\omega )(\chi _2+\chi _3)+\mu (\kappa +\mu +\sigma )\chi _1(R_{ed}^\phi -1)}, \\ E_i^*= & {} \frac{ \mu (\kappa +\mu +\sigma )\chi _1(R_{ed}^\phi -1)S^*_i}{(\mu +\rho )(\mu +\xi )(\chi _2+\chi _3)}, \\ I_i^*= & {} \frac{\rho \mu (\kappa +\mu +\sigma )\chi _1(R_{ed}^\phi -1)S^*_i}{(\epsilon +\gamma +\delta +\mu +(d+e(1-d))\theta )(\mu +\rho )(\mu +\xi )(\chi _2+\chi _3)}. \end{aligned}$$

The local stability results of equilibrium points of the system (3.13)–(3.16) is established in Theorems 5 and 6 below.

Theorem 5

(Local stability of DFE) For the system (3.13)–(3.16) in the reduced full model, the disease free equilibrium point $E^0_\phi$ is locally asymptotically stable if $R_{ed}^\phi <1$, and it is unstable if $R_{ed}^\phi >1$.

The proof is given in Appendix E. To study the stability behavior of the endemic equilibrium, we have to investigate the type of bifurcation at $R_{ed}^\phi =1$. Using central manifold theory (see, Appendix B) and the algebraic computation in Appendix F, we have established the following result.

Theorem 6

(Local stability of endemic equilibrium) If $R_{ed}^\phi >1$, then the endemic equilibrium point $\Gamma ^*_\phi$ of the reduced full model system (3.13) - (3.16) is locally asymptotically stable near to $R_{ed}^\phi =1$.

Now the relationships among the three reproduction numbers $R_\omega$, $R^\phi _0$ and $R^\phi _{ed}$ are discussed. It is clear that in the absence of both screenings (i.e, $d=0$ and $e=0$), $R^\phi _{ed}=R^\phi _0$.

Further, the following observations from equations (3.8), (3.21) and (3.22) can also be easily made:

• (i) $\frac{\partial R^\phi _{ed}}{\partial e}\le 0$, for $0\le e\le 1$

• (ii) $e=\frac{(\epsilon +\gamma +\mu +\delta )(1-d)\alpha \phi -\beta d}{\beta (1-d)}$    $\Longleftrightarrow$    $R^\phi _{ed}=R_\omega$,

• (iii) $\frac{((\epsilon +\gamma +\mu +\delta )(1-d)\alpha \phi -\beta d)}{\beta (1-d)}<e\le 1$    $\Longrightarrow$    $R^\phi _{ed}<R_\omega$

Table 2 Values of parameters used in the numerical simulation

Fig. 1

The bifurcation diagram in $\beta -e$ space with $\omega =0.12$, $\epsilon =0.23$, $\phi =0.45$ and other parameter values as in Table 2, a when $d=0$ and b when $d=0.2$

The two parameter bifurcation diagram with respect to infectivity rate ($\beta$) and entry-screening (e) are drawn in Fig. 1 for $d=0$ and $d\ne 0$. In the absence of departure screening ($d=0$) the parameter space in Fig. 1a has been divided in to five regions $(A)-(E)$ by the lines $R_\omega =1$, $R^\phi _0=1$ and $R^\phi _{ed}=1$. Their epidemiological meanings in each region are stated bellow:

(A):

$R_\omega >1$, $R^\phi _0>1$ and $R^\phi _{ed}>1$: A stable endemic equilibrium exists in both the patches and the transport-related infection increases endemicity even the presence of entry screening. It follows that, use of entry screening is not sufficient to eliminate the disease.

(B):

$R_\omega >1$, $R^\phi _0>1$, and $R^\phi _{ed}<1$: Endemic equilibrium exists and the transport-related infection will raise endemicity. But both patches may become disease free with entry screening.

(C):

$R_\omega <1$, $R^\phi _{ed}<1$, but $R^\phi _0>1$: The disease free equilibrium exists in both the patches. In absence of proper entry screening the transport-related infection may lead to endemicity.

(D):

$R_\omega <1$, $R^\phi _0<1$ and $R^\phi _{ed}<1$: The disease free equilibrium exists, and the transport-related infection may not lead to endemic case even in the absence of entry screening.

(E):

$R_\omega <1$, $R^\phi _0>1$ and $R^\phi _{ed}>1$. The disease free equilibrium exists. The transport-related infection leads the disease endemic and entry screening is not sufficient to eliminate it.

When $d\ne 0$ (particularly for $d=0.2$), due to this additional departure screening the region (E) gets vanished and the remaining four regions are present as shown in Fig. 1b. Furthermore, as d increases in $0<d\le 1$, the slope of the line $R_{ed}^\phi =1$ decrease. It can be observed that the area coverage of regions (A) and (E) gets less, and region (E) will disappear after some time. Regions (B) and (C) get increase. It means that the proper use of departure screening has its own role in eradicating the disease spreading caused by travel-related infection.

In the above analysis, the role of thresholds $R_{ed}^\phi$ and $R_\omega$ on the disease dynamics are emphasized. In Fig. 2 the variation of these thresholds with respect to awareness parameter $\omega$ and travel rate $\alpha$ are studied.

Fig. 2

The behavior of reproduction numbers related to awareness ($R_\omega$) and related to screening ($R^\phi _{ed}$) with respect to the parameters awareness ($\omega$) and travel rate of individuals ($\alpha$). When the parameter values are $\beta =0.75$, $\phi =0.45$, $d=0.2$, $e=0.5$ and other parameter values are as given in Table 2

Note that $R_\omega$ does not depend on travel parameters ($\alpha$, e, d and $\phi$), however $R_{ed}^\phi$ depend on awareness parameter $\omega$. It follows that, the following observations are made:

Both the thresholds decrease with increasing awareness. In other words, public awareness reduces the spread of disease. As travel rate ($\alpha$) increases the threshold $R_{ed}^\phi$ increases. The movement of population has an adverse effect on $R_{ed}^\phi$. For the chosen parameters ($\beta$, $\phi$, d and e) and the range space of $\omega$ and $\alpha$:

$$\begin{aligned} R_\omega >R_{ed}^\phi \end{aligned}$$

(3.24)

This is due to restriction (iii) on e. If entry screening is not monitored and condition (iii) is violated, then the inequality is reversed in the presence of travel-related infection ($\phi \ne 0$). If $\phi =0$, the restriction is trivially satisfied, and the above inequality holds for all combinations of nonzero e and d (refer Sect. 3.2.1).

3.3 Sensitivity analysis

Sensitivity analysis is important for determining which parameters are more affecting in reducing or expand the level of disease. In epidemiological models, the size of the reproduction number illustrates the power of the disease to raise or eradicate from the population. In this section, we use a direct differentiation method applied in many literature [5, 7, 23, 27] used the name as the normalized sensitivity index or the elasticity index. For a variable ${\mathcal {U}}$ and a parameter ${\mathcal {P}}$, it is given as:

$$\begin{aligned} \Upsilon ^{\mathcal {U}}_{\mathcal {P}}=\frac{\partial {\mathcal {U}}}{\partial {\mathcal {P}}}\times \frac{{\mathcal {P}}}{{\mathcal {U}}} \end{aligned}$$

(3.25)

This normalized sensitivity index gives the proportional rate of change of ${\mathcal {U}}$ as ${\mathcal {P}}$ changes. Now, the variable ${\mathcal {U}}$ represent $R^\phi _{ed}$ which is a reproduction number given in (3.21). The reproduction number is inversely affected by the parameters $\omega$, $\epsilon$, $\gamma$, $\mu$, $\delta$, e and d described in Table 1 with values in Table 2 and their normalized sensitivity indices are negative as shown in Table 3. The parameters $\xi$, $\rho$, $\beta$, $\alpha$ and $\phi$ have positive indices and affected the reproduction number directly.

Table 3 Normalized sensitivity indices of a reproduction number ($R^\phi _{ed}$) with respect to several parameters $\beta =0.75$, $\phi =0.45$, $d=0.2$, $e=0.5$ and other parameter values as in Table 2

4 Numerical simulation

In this section the non travel model system (3.1)–(3.4) and the reduced full model (3.13) - (3.16) are solved numerically for the parameters as specified in Table 2. Numerical simulations are performed in MatLab using ODE45 package. The dynamic behaviors of exposed, infectious, and quarantined populations are explored.

In the absence of travel, the patches are isolated. Due to limited resources and communication barriers, awareness may not be addressed through the overall population, and quarantine may not be accessible. In absence of awareness and internal quarantine ($\omega =\epsilon =0$), the basic reproduction number $R_\omega$ is computed as 7.2079 for the data in Table 2 with $\beta = 0.75$ using equation (3.8). The disease is endemic in isolated patches with ($E^*=288.59$, $I^*=832.18$, $Q^*=0.00$). The time series plot of Exposed (E), infectious (I) and quarantined (Q) are drawn in Fig. 3. The four cases are considered depending on awareness ($\omega$) and quarantine measure ($\epsilon$). The Fig. 3a confirms the convergence to endemic state when $\epsilon =0$ and $\omega =0$. In absence of awareness, the load of infection is reduced when flow rate to quarantine class $\epsilon = 0.35$, see Fig. 3b. Similarly, in Fig. 3c, the infection load is reduced with awareness, $\omega =0.3$ in absence of quarantine measures, $\epsilon =0$. In both these cases the endemic state is achieved. However, when awareness is combined with suitable quarantine measures, the disease can be eliminated. This is confirmed in Fig. 3d.

Fig. 3

The dynamical behavior of a model with no travel given in (3.1)–(3.4), $\beta =0.75$, a when $\omega =0, \epsilon =0$, b when $\omega =0, \epsilon =0.35$, c when $\omega =0.3, \epsilon =0$, d when $\omega =0.3, \epsilon =0.35$ and other parameter values are as given in Table 2. $E=$exposed, $I=$infectious and $Q=$quarantined populations

The full model (3.13)–(3.16) simulations with the same level of migration between the two patches are presented in Figs. 4, 5 and 6. The two patches are symmetric with respect to model parameters. The numerical simulations assume that the disease begins in patch 1 and no infectious individual is observed initially in patch 2. The Fig. 4 shows, the time series for exposed, infectious and quarantined populations considering no awareness and without transport-related infection for the choice of parameters such that $R_\omega >1$. The exposed, infectious, and quarantined populations of the two patches go to the endemic point ($E^*=200.59$, $I^*=121.80$, $Q=195.38$). The transport-related infection ($\phi =0.85$) increases the peak value and their endemicity level is ($E^*=249.07$, $I^*=151.25$, $Q=242.60$) as evident from Fig. 5 without awareness in the two patches. Accordingly, the disease invades in both patches. This happens when entry and departure screenings are carried out. These results are consistent with the analytical findings in Sect. 3. However, when the population is aware in both the patches with rate $\omega =0.3$, the reproduction number is computed as $R_{ed}^\phi =0.8054$. It follows that the exposed, infectious and quarantined populations tend to a stable disease-free state, since $R_{ed}^\phi <1$ when the travel-related infection is considered ($\phi =0.85$) [see Fig. 6]. Due to differences in initial conditions, the level and evolution pattern of the infected population initially are different in the two patches. The disease is ultimately eradicated from both patches. This is in line with the analytical result established in Theorem 5.

Fig. 4

The dynamical behavior of exposed, infectious and quarantined populations in the model system (3.13)–(3.16) when $\beta =0.75$, $d=0.2$, $e=0.5$, $\epsilon =0.15$, $\omega =0$, $\phi =0$ and other parameters are given in Table 2 resulting $R_{ed}^\phi =1.5179>1$. $E=$exposed, $I=$infectious and $Q=$quarantined populations. The respective solid lines refer patch 1 and doted lines refer patch 2

Fig. 5

The dynamical behavior of exposed, infectious and quarantined populations in the model system (3.13)–(3.16) when $\beta =0.75$, $d=0.2$, $e=0.5$, $\epsilon =0.15$, $\omega =0$, $\phi =0.85$ and other parameters are given in Table 2 resulting $R_{ed}^\phi =2.0133>1$. $E=$exposed, $I=$infectious and $Q=$quarantined populations. The respective solid lines refer patch 1 and doted lines refer patch 2

Fig. 6

The dynamical behavior of exposed, infectious and quarantined populations in the model system (3.13)–(3.16) when $\beta =0.75$, $d=0.2$, $e=0.5$, $\epsilon =0.15$, $\omega =0.3$, $\phi =0.85$ and other parameters are given in Table 2 resulting $R_{ed}^\phi =0.8054<1$. $E=$exposed, $I=$infectious and $Q=$quarantined populations. The respective solid lines refer patch 1 and doted lines refer patch 2

In Figs. (7, 8, 9 and 10) the effect of internal quarantine and different migration levels between the two patches in the dynamics of the full model (3.13)–(3.16) are discussed.

Fig. 7

The dynamical behavior of exposed, infectious and quarantined populations with no internal quarantine rate ($\epsilon _1=\epsilon _2=0$) in the model system (3.13)– (3.16) when $\beta =0.75$, $\theta _{ij}=\alpha _{ij}=0$, for $i, j=1, 2$, $i\ne j$ and other parameters are given in Table 2 resulting $R_{ed}^\phi = 2.8835>1$. $E=$exposed, $I=$infectious, $Q=$quarantined populations. The respective solid lines refer patch 1 and doted lines refer patch 2

Fig. 8

The dynamical behavior of exposed, infectious and quarantined populations with no internal quarantine rate ($\epsilon _1=\epsilon _2=0$) and unidirectional migration ($\theta _{12}=0.15$, $\alpha _{12}=0.5$, $\theta _{21}=\alpha _{21}=0$) in the model system (3.13)–(3.16) when $\beta =0.75$, $d=0.2$, $e=0.5$, $\phi =0.85$ and other parameters are given in Table 2. $E=$exposed and $I=$infectious populations. The respective solid lines refer patch 1 and doted lines refer patch 2

Fig. 9

The dynamical behavior of exposed, infectious and quarantined populations with no internal quarantine rate ($\epsilon _1=\epsilon _2=0$) and two side migration ($\theta _{12}=0.15$, $\alpha _{12}=0.5$, $\theta _{21}=0.25$, $\alpha _{21}=0.75$) in the model system (3.13)–(3.16) when $\beta =0.75$, $d=0.2$, $e=0.5$, $\phi =0.85$. Other parameters are given in Table 2. $E=$exposed, $I=$infectious and $Q=$ quarantined populations. The respective solid lines refer patch 1 and doted lines refer patch 2

In Fig. 7 the result of non-migration ($\theta _{ij}=\alpha _{ij}=0$ for $i,j=1,2$, $i\ne j$) without internal quarantine ($\epsilon _1=\epsilon _2=0$) is shown. For zero migration, our choice of screening and travel-induced infection does not affect the disease negatively or positively, but due to the absence of quarantine, the reproduction number exceeds unity ($R_\omega =2.8835>1$). It follows that the absence of internal quarantine leads to the disease endemic in the two patches. It is consistent with the result in Fig. 3c.

When migration is allowed only from patch 1 to patch 2 ($\theta _{12}=0.15$, $\alpha _{12}=0.5$ and $\theta _{21}=\alpha _{21}=0$) with no internal quarantine, undetected infectious individuals are going out from patch 1. It follows from Fig. 8 that the patch 1 becomes disease-free due to the emigration of infectious individuals. However, the disease is still endemic in patch 2, and the transport-related infection increases the endemic state to ($E^*=441.30$, $I^*=1272$). The appearance of quarantined population in the two patches is due to the detection at the departure/entry screenings time. In the two side migration ($\theta _{12}=0.15$, $\alpha _{12}=0.5$, $\theta _{21}=0.25$ and $\alpha _{21}=0.75$) it is seen in Fig. 9 that the disease remain endemic in both the patches. Further, the load of infection decreases to the endemic state ($E^*=128.4$, $I^*=159.4$) in patch 2 due to departure/entry screening-related quarantine. In Fig. 10, the possible reduction of a load of infection is observed in both the patches as the entry screening becomes more strict with increasing e.

Fig. 10

The dynamical behavior of infectious population with no internal quarantine rate ($\epsilon _1=\epsilon _2=0$) and two side migration ($\theta _{12}=0.15$, $\alpha _{12}=0.5$, $\theta _{21}=0.25$, $\alpha _{21}=0.75$) for increasing e in the model system (3.13)–(3.16) when $\beta =0.75$, $d=0.2$, $\phi =0.85$. other parameters are given in Table 2

Figures 11, 12 and 13 presents the dynamical behavior of the full model system for limited internal quarantine ($\epsilon _1=\epsilon _2=0.05$) and in different level of migration.

Fig. 11

The dynamical behavior of exposed, infectious and quarantined populations with internal quarantine rate ($\epsilon _1=\epsilon _2=0.05$) in the model system (3.13)–(3.16) when $\beta =0.75$, $\theta _{ij}=\alpha _{ij}=0$, for $i, j=1, 2$, $i\ne j$ and other parameters are given in Table 2 resulting $R_\omega =1.9475>1$. $E=$exposed, $I=$infectious and $Q=$ quarantined populations. The respective solid lines refer patch 1 and doted lines refer patch 2

Fig. 12

The dynamical behavior of exposed, infectious and quarantined populations with internal quarantine rate ($\epsilon _1=\epsilon _2=0.05$) and unidirectional migration ($\theta _{12}=0.15$, $\alpha _{12}=0.5$, $\theta _{21}=\alpha _{21}=0$) in the model system (3.13)–(3.16) when $\beta =0.75$, $d=0.2$, $e=0.5$, $\phi =0.85$ and other parameters are given in Table 2. $E=$exposed, $I=$infectious and $Q=$ quarantined populations

Fig. 13

The dynamical behavior of exposed, infectious and quarantined populations with internal quarantine rate and two side migration ($\theta _{12}=0.15$, $\alpha _{12}=0.5$, $\theta _{21}=0.25$, $\alpha _{21}=0.75$) in the model system (3.13)–(3.16) when $\beta =0.75$, $d=0.2$, $e=0.5$, $\phi =0.85$. Other parameters are given in Table 2. (a) when $\epsilon _1=\epsilon _2=0.05$, (b) when $\epsilon _1=\epsilon _2=0.35$. $E=$exposed, $I=$infectious and $Q=$ quarantined populations. The respective solid lines refer patch 1 and doted lines refer patch 2

When there is no migration between patches, the disease is endemic with their endemic value ($E^*=190.8$, $I^*=371.7$, $Q^*=76.44$) in both the patches as shown in Fig. 11. However, if people are allowed to move from patch 1 to patch 2 only, then patch 1 becomes disease-free but patch 2 remains endemic, and its load is also increasing. This is illustrated in Fig. 12. The bi-directional inter-patch movement decreases the infection load in patch 2, and both the patches are endemic [see Fig. 13a]. The two patches may become disease-free when more effort is applied to internal quarantine ($\epsilon _1=\epsilon _2=0.35$) [see Fig. 13b].

In Figs.14, 15, 16 and 17 it is seen that, the dynamical behavior of exposed, infectious and quarantined populations for the patches with different model parameters ($\beta _1=0.5$, $\beta _2=0.75$, $\omega _1=0.01$, $\omega _2=0.03$, $\epsilon _1=0.15$, $\epsilon _2=0.35$) with varying migration rates.

Fig. 14

The dynamical behavior of exposed, infectious and quarantined populations with no migration ($\theta _{ij}=\alpha _{ij}=0$, for $i, j=1, 2$, $i\ne j$) and different parameter values ($\beta _1=0.5$, $\beta _2=0.75$, $\omega _1=0.01$, $\omega _2=0.03$, $\epsilon _1=0.15$, $\epsilon _2=0.35$) and other parameters are given in Table 2. $E=$exposed, $I=$infectious and $Q=$ quarantined populations. The respective solid lines refer patch 1 and doted lines refer patch 2

Fig. 15

The dynamical behavior of exposed, infectious and quarantined populations with unidirectional migration from patch 1 to patch 2 ($\theta _{12}=0.15$, $\alpha _{12}=0.5$, $\theta _{21}=\alpha _{21}=0$) and different parameter values ($\beta _1=0.5$, $\beta _2=0.75$, $\omega _1=0.01$, $\omega _2=0.03$, $\epsilon _1=0.15$, $\epsilon _2=0.35$, $d=0.2$, $e=0.5$, $\phi =0.85$) and other parameters are given in Table 2. $E=$exposed, $I=$infectious and $Q=$ quarantined populations. The respective solid lines refer patch 1 and doted lines refer patch 2

Fig. 16

The dynamical behavior of exposed, infectious and quarantined populations with unidirectional migration from patch 2 to patch 1 ($\theta _{12}=\alpha _{12}=0$, $\theta _{21}=0.25$, $\alpha _{21}=0.75$) and different parameter values ($\beta _1=0.5$, $\beta _2=0.75$, $\omega _1=0.01$, $\omega _2=0.03$, $\epsilon _1=0.15$, $\epsilon _2=0.35$, $d=0.2$, $e=0.5$, $\phi =0.85$) and other parameters are given in Table 2. $E=$exposed, $I=$infectious and $Q=$ quarantined populations. The respective solid lines refer patch 1 and doted lines refer patch 2

Fig. 17

The dynamical behavior of exposed, infectious and quarantined populations with two side migration ($\theta _{12}=0.15$, $\alpha _{12}=0.5$, $\theta _{21}=0.25$, $\alpha _{21}=0.75$) and different parameter values ($\beta _1=0.5$, $\beta _2=0.75$, $\omega _1=0.01$, $\omega _2=0.03$, $\epsilon _1=0.15$, $\epsilon _2=0.35$, $d=0.2$, $e=0.5$, $\phi =0.85$) and other parameters are given in Table 2. $E=$exposed, $I=$infectious and $Q=$ quarantined populations. The respective solid lines refer patch 1 and doted lines refer patch 2

In the absence of migration, both the patches are endemic for the choice of parameters as shown in Fig. 14. When unidirectional migration is allowed from patch 1 to patch 2 ($\theta _{12}=0.15$, $\alpha _{12}=0.5$, $\theta _{21}=\alpha _{21}=0$), undetected infectious individuals may migrate to patch 2. It follows that the disease disappears from the community in patch 1. However, with infection flow from patch 1 the disease remains endemic in patch 2 and its endemic value goes to ($E^*=366.4$, $I^*=242.1$, $Q^*=348.5$) (refer Fig. 15). When the migration is from patch 2 to patch 1 ($\theta _{21}=0.25$, $\alpha _{21}=0.75$ and $\theta _{12}=\alpha _{12}=0$), patch 1 is endemic and patch 2 is disease free as shown in Fig. 16. If inter-patch migration is allowed between patch 1 and patch 2 ($\theta _{12}=0.15$, $\alpha _{12}=0.5$ $\theta _{21}=0.25$ and $\alpha _{21}=0.75$), migration may cause the disease to become endemic in both the patches as seen in Fig. 17. The presence of awareness and internal quarantine including detection from entry/departure screenings are insufficient to remove the disease from the community.

5 Conclusion

In this paper, an $n-$patch SWEIQR epidemic model is developed by considering departure-entry screening to study the spread of disease during inter-patch movement. The role of awareness, quarantine, and entry/departure screening on the spread of disease is investigated.

First, the model with $n=2$ and without movement is analyzed. The model admits a disease-free and an endemic equilibrium state. The disease-free equilibrium is locally stable when the reproduction number ($R_\omega$) is less than unity and unstable if it exceeds unity ($R_\omega >1$). The trans-critical bifurcation occurs and a local stable endemic equilibrium exists near $R_\omega =1$ if the reproduction number $R_\omega >1$. Furthermore, the disease-free equilibrium point is globally asymptotically stable if $R_\omega <1$. It is also observed that awareness in the population may reduce the infection load or bring the endemic state to a disease-free state, $R_\omega \le R_0$.

In model (3.13)–(3.16), the model with transport-related infection, awareness and screening may be applicable to SARS, Nipah, Swine-flu, or Covid-19 pandemic. The quarantine, screening, and creating awareness are essential tools to reduce spreading. Researchers indicate that entry screening has a possibility to eradicate the spread of disease led by transport-related infection. The analytical and numerical analysis of the model shows that the disease-free equilibrium is locally stable when $R^\phi _{ed}<1$ and it is unstable if $R^\phi _{ed}>1$. At $R^\phi _{ed}=1$ the trans-critical bifurcation point occur; consequently the local asymptotically stable endemic point appear when $R^\phi _{ed}>1$. It is also shown that awareness and departure screening may have the possibility to reduce the disease led by transport-related infection. The simulation results are in agreement with the field observations.

Appendices

Appendices

Proof of Theorem 2

Let $X_i=(S_i, W_I, E_i, I_i)^T$. The Jacobian matrix of the system (3.1)–(3.4) is given by:

$$\begin{aligned} J(X_i)= \begin{pmatrix} -\dfrac{\beta I_i}{N_i}-(\omega +\mu ) &{} \xi &{} 0 &{} -\dfrac{\beta S_i}{N_i}\\ \omega &{} -(\xi +\mu ) &{} 0&{} 0\\ \dfrac{\beta I_i}{N_i} &{} 0 &{} -(\rho +\mu ) &{} \dfrac{\beta S_i}{N_i}\\ 0 &{} 0 &{}\rho &{} -(\epsilon +\gamma +\mu +\delta ) \end{pmatrix} \end{aligned}$$

(A.1)

At the disease free equilibrium point ($E^{0}$), the Jacobian matrix A.1 becomes:

$$\begin{aligned} J(E^{0})= \begin{pmatrix} -(\omega +\mu ) &{} \xi &{} 0 &{} -\dfrac{\beta (\xi +\mu )}{(\xi +\mu +\omega )} \\ \omega &{} -(\xi +\mu ) &{} 0 &{} 0\\ 0 &{} 0 &{} -(\rho +\mu ) &{} \dfrac{\beta (\xi +\mu )}{(\xi +\mu +\omega )}\\ 0 &{} 0 &{}\rho &{} -(\epsilon +\gamma +\mu +\delta ) \end{pmatrix} \end{aligned}$$

(A.2)

The Characteristic polynomial for $J(E^0)$ is computed as

$$\begin{aligned} (\lambda +\mu )(\lambda +(\xi +\mu +\omega ))(\lambda ^2+A\lambda +B)=0 \end{aligned}$$

(A.3)

Two of its eigenvalues are: $\lambda _1=-\mu$, $\lambda _2=-(\xi +\mu +\omega )$

The remaining two eigenvalues are the roots of the quadratic equation:

$$\begin{aligned} \lambda ^{2}+A\lambda +B=0 \end{aligned}$$

(A.4)

$$\begin{aligned} A= & {} (\epsilon +\delta +\gamma +\rho +2\mu )>0\\ B= & {} (\rho +\mu )(\epsilon +\delta +\mu +\gamma )(1-R_\omega )>0, \quad for \quad R_\omega <1 \end{aligned}$$

This shows that all eigenvalues of the Jacobian matrix A.2 are either negative real or have negative real parts for $R_\omega <1$. Thus, The disease-free equilibrium point is locally stable on the system (3.1)–(3.4) if $R_\omega <1$. It is unstable if $R_\omega >1$, since the characteristic equation has one positive eigenvalue.

Central manifold theory

[3] Consider the following general system of ordinary differential equations with a bifurcation parameter $\phi$,

$$\begin{aligned} \dfrac{dx}{dt}=f(x,\phi ), f:{\mathbb {R}}^n\times {\mathbb {R}}\rightarrow {\mathbb {R}}^n\quad and\quad f\in {\mathbb {C}}^2 ({\mathbb {R}}^n\times {\mathbb {R}}) \end{aligned}$$

(B.1)

Without loss of generality, it is assumed that 0 is an equilibrium for the system B.1 , that is $f(0,\phi )\equiv 0$, for all $\phi$.

Assume that:

A1 : :

$A = D_xf(0, 0) =\dfrac{\partial f_i}{\partial x_j}(0, 0)$ is the linearization matrix of System B.1 around the equilibrium 0 with $\phi =0$. Zero is a simple eigenvalue of A and all other eigenvalues of A have negative real parts;

A2 : :

Matrix A has a non-negative right eigenvector w and a left eigenvector v corresponding to the zero eigenvalue.

Let $f_k$ be the $k^{th}$ component of f and

$$\begin{aligned}&a=\sum _{k,i,j=1}^nv_kw_iw_j\dfrac{\partial ^2f_k}{\partial x_i\partial x_j}(0, 0) \end{aligned}$$

(B.2)

$$\begin{aligned}&b=\sum _{k,i=1}^nv_kw_i\dfrac{\partial ^2f_k}{\partial x_i\partial \phi }(0, 0) \end{aligned}$$

(B.3)

The local dynamics of B.1 around 0 are totally determined by a and b given in B.2 and B.3 respectively.

• i. $a> 0, b > 0$. When $\phi <0$ with $|\phi |<<1$, 0 is locally asymptotically stable, and there exists a positive unstable equilibrium; when $0<\phi<<1$, 0 is unstable and there exists a negative and locally asymptotically stable equilibrium;

• ii. $a<0, b<0$. When $\phi <0$ with $|\phi |<< 1$, 0 is unstable; when $0<\phi<<1$, 0 is locally asymptotically stable, and there exists a positive unstable equilibrium;

• iii. $a>0, b<0$. When $\phi <0$ with $|\phi |<<1$, 0 is unstable, and there exists a locally asymptotically stable negative equilibrium; when $0<\phi<<1$, 0 is stable, and a positive unstable equilibrium appears;

• iv. $a<0, b>0$. When $\phi$ changes from negative to positive, 0 changes its stability from stable to unstable. Correspondingly a negative unstable equilibrium becomes positive and locally asymptotically stable.

Proof of Theorem 3

Set the state variables of the model (3.1) - (3.4) as $S_i=x_1$, $W_i=x_2$, $E_i=x_3$, $I_i=x_4$ and we can rewrite the system as follows:

$$\begin{aligned} \frac{dx_1}{dt}&=\varLambda -\dfrac{\beta x_1x_4}{N} -(\omega +\mu )x_1+\xi x_2:=f_1 \end{aligned}$$

(C.1)

$$\begin{aligned} \frac{dx_2}{dt}&=\omega x_1 -(\xi + \mu )x_2:=f_2 \end{aligned}$$

(C.2)

$$\begin{aligned} \frac{dx_3}{dt}&=\dfrac{\beta x_1x_4}{N}-(\rho + \mu )x_3:=f_3 \end{aligned}$$

(C.3)

$$\begin{aligned} \frac{dx_4}{dt}&=\rho x_3 -(\epsilon +\mu +\delta +\gamma )x_4:=f_4 \end{aligned}$$

(C.4)

For $R_\omega =1$ the corresponding $\beta =\beta ^*$ is assumed to be a bifurcation parameter.

It is note that,

$$\begin{aligned} \beta ^*=\dfrac{(\xi +\mu +\omega )(\rho +\mu )(\epsilon +\gamma +\mu +\delta )}{\rho (\xi +\mu )} \end{aligned}$$

The disease-free-equilibrium is $X^0=\left( x^0_1,x^0_2,x^0_3,x^0_4\right)$, where

$$\begin{aligned} x^0_1=\dfrac{\varLambda (\xi +\mu )}{\mu (\xi +\mu +\omega )}, x^0_2=\dfrac{\varLambda \omega }{\mu (\xi +\mu +\omega )}, x^0_3=x^0_4=0. \end{aligned}$$

The linearization matrix of system C.1 - C.4 around the disease-free-equilibrium when $\beta =\beta ^*$ is:

$$\begin{aligned} D_xf(X^0, \beta ^*)= \begin{pmatrix} -(\omega +\mu ) &{} \xi &{} 0 &{} -\dfrac{\beta ^*(\xi +\mu )}{(\xi +\mu +\omega )} \\ \omega &{} -(\xi +\mu ) &{} 0 &{} 0 \\ 0 &{} 0 &{} -(\rho +\mu ) &{} \dfrac{\beta ^*(\xi +\mu )}{(\xi +\mu +\omega )} \\ 0 &{} 0 &{}\rho &{} -(\epsilon +\gamma +\mu +\delta ) \end{pmatrix} \end{aligned}$$

(C.7)

It is clear that, $D_xf(X^0, \beta ^*)$ has a simple zero eigenvalue and a right eigenvector corresponding to the zero eigenvalue is $w=(w_1, w_2, w_3, w_4)^T$, where

$$\begin{aligned}&w_1= -\dfrac{(\mu +\xi )(\mu +\rho )(\gamma +\delta +\epsilon +\mu )}{\mu (\mu +\xi +\omega )},\quad w_2=-\dfrac{\omega (\mu +\rho )(\gamma +\delta +\epsilon +\mu )}{\mu (\mu +\xi +\omega )},\\&\quad w_3=(\gamma +\delta +\epsilon +\mu ),\quad w_4= \rho \end{aligned}$$

and the left eigenvector associated with the zero eigenvalue satisfying $w.v=1$ is $v=(v_1, v_2, v_3, v_4)$, where,

$$\begin{aligned} v_1= & {} 0,\quad v_2=0,\quad v_3=\dfrac{1}{(\gamma +\delta +\epsilon +2\mu +\rho )},\quad and\quad \nonumber \\ v_4= & {} \dfrac{(\mu +\rho )}{\rho (\gamma +\delta +\epsilon +2\mu +\rho )}. \end{aligned}$$

Based on the theoretical result given in B , the expressions for a and b at $(X^0, \beta ^*)$ are computed as

$$\begin{aligned} a&=\sum _{k,i,j=1}^4v_kw_iw_j\dfrac{\partial ^2f_k}{\partial x_i\partial x_j}(X^0, \beta ^*) \end{aligned}$$

(C.8)

$$\begin{aligned} b&=\sum _{k,i=1}^4v_kw_i\dfrac{\partial ^2f_k}{\partial x_i\partial \beta }(X^0, \beta ^*) \end{aligned}$$

(C.9)

Since $v_1=v_2=0$ and $v_3=v_4\ne 0$, equations C.8 and C.9 simplified as

$$\begin{aligned} a&=v_3\sum _{i,j=1}^4w_iw_j\dfrac{\partial ^2f_3}{\partial x_i\partial x_j}(X^0, \beta ^*)+v_4\sum _{i,j=1}^4w_iw_j\dfrac{\partial ^2f_4}{\partial x_i\partial x_j}(X^0, \beta ^*) \end{aligned}$$

(C.10)

$$\begin{aligned} b&=v_3\sum _{i=1}^4w_i\dfrac{\partial ^2f_3}{\partial x_i\beta }(X^0, \beta ^*)+v_4\sum _{i=1}^4w_i\dfrac{\partial ^2f_4}{\partial x_i\partial \beta }(X^0, \beta ^*) \end{aligned}$$

(C.11)

Substituting the eigenvectors and the computed partial derivatives of the system C.1 - C.4 at $(X^0, \beta ^*)$ in the formula for a in C.10 and b in C.11 , after some algebraic computation, yields

$$\begin{aligned}&a=-\dfrac{2\omega (\mu +\rho )^2(\epsilon +\gamma +\mu +\delta )^2}{\varLambda (\xi +\mu +\omega )(\epsilon +\gamma +2\mu +\delta +\rho )}<0 \text { and}\\&b=\dfrac{\rho (\mu +\xi )}{(\xi +\mu +\omega )(\epsilon +\gamma +2\mu +\delta +\rho )}>0 \end{aligned}$$

The values of a and b indicates that, at $\beta =\beta ^*$ for $R_\omega =1$, the system exhibits a trans-critical bifurcation (i.e the disease free equilibrium point changes its stability from locally stable for $R_\omega <1$ to unstable for $R_\omega >1$) and the endemic equilibrium $\Gamma ^*$ is locally stable.

Proof of Theorem 4

Let $V_i(t)$ be a positive definite function in $\varOmega ^*$. For a positive value a, $V_i(t)$ in the ith patch is defined in equation 3.12 and its derivative is:

$$\begin{aligned} \dfrac{dV_i(t)}{dt}&=\dfrac{d(a E_i(t))}{dt}+\dfrac{dI_i(t)}{dt}\\ \\&=a\left[ \dfrac{\beta S_i(t)I_i(t)}{N_i(t)}-(\rho + \mu )E_i(t)\right] + \left[ \rho E_i(t) -(\epsilon +\mu +\delta +\gamma )I_i(t)\right] \\ \\&\le \left[ a\beta -(\epsilon +\mu +\delta +\gamma )\right] I_i(t)-\left[ a(\rho +\mu )-\rho \right] E_i(t), \quad \nonumber \\&\quad Since\quad S_i(t)\le N_i(t)\\ \end{aligned}$$

Determine the value of a such that,

$$\begin{aligned}&a(\rho +\mu )-\rho =0 \\&\Longrightarrow a=\dfrac{\rho }{\rho +\mu } \end{aligned}$$

It follows that,

$$\begin{aligned} \dfrac{dV_i(t)}{dt}&\le (\epsilon +\mu +\delta +\gamma )\left( \dfrac{\beta \rho }{(\rho +\mu )(\epsilon +\mu +\delta +\gamma )}-1\right) I_i(t)\\&=(\epsilon +\mu +\delta +\gamma )\left( R_0-1\right) I_i(t), \quad since\quad R_0=\dfrac{\beta \rho }{(\rho +\mu )(\epsilon +\mu +\delta +\gamma )} \end{aligned}$$

Thus, for $R_0\le 1$, we have $dV_i(t)/dt\le 0$. It is also true for $R_\omega \le 1$, since $R_\omega \le R_0$ [see (3.8) and (3.9)]. Moreover, $dV_i(t)/dt=0$ if $I_i(t)=0$. Therefore, $V_i(t)$ is a Lyapunov function on $\varOmega ^*$ defined in (3.6). Substituting $I_i(t)=0$ in the system (3.1)–(3.4) and computing gives

$$\begin{aligned}&S_i(t)\rightarrow \dfrac{\varLambda (\xi +\mu )}{\mu (\xi +\mu +\omega )},\quad and \\&W_i(t)\rightarrow \dfrac{\varLambda \omega }{\mu (\xi +\mu +\omega )}\quad as\quad t\rightarrow \infty . \end{aligned}$$

Hence, the largest invariant set in $\left\{ (S_i, W_i, E_i, I_i)\in \varOmega ^*:\dfrac{dV_i}{dt}=0\right\}$ is the singleton $\{E^0\}$ whenever $R_\omega \le 1$. It follows that by Lasalle’s invariance principle, the solution of the system (3.1) - (3.4) with initial conditions (2.7) in $\varOmega ^*$ approaches $\{E^0\}$ as $t\longrightarrow \infty$.

Proof of Theorem 5

The Jacobian matrix of the system in a reduced full infection model (3.13)–(3.16) is given by the following block matrix:

$$\begin{aligned} J=\begin{pmatrix} A_1 &{} B_2\\ B_1 &{} A_2 \end{pmatrix} \end{aligned}$$

(E.1)

Where,

$$\begin{aligned}&A_1=\begin{pmatrix} -\dfrac{\beta I_1}{N_1}-(\alpha +\mu +\omega ) &{} \xi &{} 0 &{} -\dfrac{\beta S_1}{N_1}\\ \omega &{} -(\xi +\mu +\alpha ) &{} 0 &{} 0 \\ \dfrac{\beta I_1}{N_1} &{} 0 &{} -(\rho +\mu +\alpha ) &{} \dfrac{\beta S_1}{N_1}\\ 0 &{} 0 &{} \rho &{} -(\epsilon +\gamma +\mu +\delta +\theta ) \end{pmatrix} \nonumber \\\end{aligned}$$

(E.2)

$$\begin{aligned}&A_2=\begin{pmatrix} -\dfrac{\beta I_2}{N_2}-(\alpha +\mu +\omega ) &{} \xi &{} 0 &{} -\dfrac{\beta S_2}{N_2} \\ \omega &{} -(\xi +\mu +\alpha ) &{} 0 &{} 0\\ \dfrac{\beta I_2}{N_2} &{} 0 &{} -(\rho +\mu +\alpha ) &{} \dfrac{\beta S_2}{N_2}\\ 0 &{} 0 &{} \rho &{} -(\epsilon +\gamma +\mu +\delta +\theta ) \end{pmatrix} \nonumber \\ \end{aligned}$$

(E.3)

$$\begin{aligned}&B_1=\begin{pmatrix} \alpha -\dfrac{(1-d)\alpha \theta \phi I_1}{N_1} &{} 0 &{} 0 &{} -\dfrac{(1-d)\alpha \theta \phi S_1}{N_1}\\ 0 &{} \alpha &{} 0 &{} 0 \\ \dfrac{(1-d)\alpha \theta \phi I_1}{N_1} &{} 0 &{} \alpha &{} \dfrac{(1-d)\alpha \theta \phi S_1}{N_1}\\ 0 &{} 0 &{} 0 &{} (1-e)(1-d)\theta \end{pmatrix} \end{aligned}$$

(E.4)

and

$$\begin{aligned} B_2=\begin{pmatrix} \alpha -\dfrac{(1-d)\alpha \theta \phi I_2}{N_2} &{} 0 &{} 0 &{} -\dfrac{(1-d)\alpha \theta \phi S_2}{N_2}\\ 0 &{} \alpha &{} 0 &{} 0\\ \dfrac{(1-d)\alpha \theta \phi I_2}{N_2} &{} 0 &{} \alpha &{}\dfrac{(1-d)\alpha \theta \phi S_2}{N_2}\\ 0 &{} 0 &{} 0 &{} (1-e)(1-d)\theta \end{pmatrix} \end{aligned}$$

(E.5)

The equilibrium points of the two patches are the same (set as $S_1^*=S_2^*=S^*$, $W_1^*=W_2^*=W^*$, $E_1^*=E_2^*=E^*$, $I_1^*=I_2^*=I^*$). It follows that, $A_1=A_2$ (say $A^*$) and $B_1=B_2$ (say $B^*$), and the Jacobian matrix E.1 , at the equilibrium point becomes;

$$\begin{aligned} J=\begin{pmatrix} A^* &{} B^*\\ B^* &{} A^* \end{pmatrix} \end{aligned}$$

(E.6)

Where

$$\begin{aligned} A^*=\begin{pmatrix} -\dfrac{\beta I^*}{N^*}-(\alpha +\mu +\omega ) &{} \xi &{} 0 &{} -\dfrac{\beta S^*}{N^*}\\ \omega &{} -(\xi +\mu +\alpha ) &{} 0 &{} 0\\ \dfrac{\beta I^*}{N^*} &{} 0 &{} -(\rho +\mu +\alpha ) &{} \dfrac{\beta S^*}{N^*}\\ 0 &{} 0 &{} \rho &{} -(\epsilon +\gamma +\mu +\delta +\theta ) \end{pmatrix}\nonumber \\ \end{aligned}$$

(E.7)

and

$$\begin{aligned} B^*=\begin{pmatrix} \dfrac{\alpha -(1-d)\alpha \theta \phi I^*}{N^*} &{} 0 &{} 0 &{} -\dfrac{(1-d)\alpha \theta \phi S^*}{N^*} \\ 0 &{} \alpha &{} 0 &{} 0\\ \dfrac{(1-d)\alpha \theta \phi I^*}{N^*} &{} 0 &{} \alpha &{} \dfrac{(1-d)\alpha \theta \phi S^*}{N^*} \\ 0 &{} 0 &{} 0 &{} (1-e)(1-d)\theta \end{pmatrix} \end{aligned}$$

(E.8)

By [4, 29] the eigenvalues of the Jacobian matrix J given in E.6 at equilibrium point of the system (3.13)–(3.16) are identical with that of the eigenvalues of $A^*+B^*$ and $A^*-B^*$. Evaluating E.6 at disease free equilibrium gives:

$$\begin{aligned} J(E^0_\phi )=\begin{pmatrix} A^0 &{} B^0\\ B^0 &{} A^0 \end{pmatrix}. \end{aligned}$$

Where,

$$\begin{aligned} A^0=\begin{pmatrix} -(\alpha +\mu +\omega ) &{} \xi &{} 0 &{} -\dfrac{\beta (\xi +\mu )}{(\xi +\mu +\omega )}\\ \omega &{} -(\xi +\mu +\alpha ) &{} 0 &{} 0\\ 0 &{} 0 &{} -(\rho +\mu +\alpha ) &{} \dfrac{\beta (\xi +\mu )}{(\xi +\mu +\omega )}\\ 0 &{} 0 &{} \rho &{} -(\epsilon +\gamma +\mu +\delta +\theta ) \end{pmatrix} \end{aligned}$$

and

$$\begin{aligned} B^0=\begin{pmatrix} \alpha &{} 0 &{} 0 &{} -\dfrac{(1-d)\alpha \theta \phi (\xi +\mu )}{(\xi +\mu +\omega )} \\ 0 &{} \alpha &{} 0 &{} 0\\ 0 &{} 0 &{} \alpha &{} \dfrac{(1-d)\alpha \theta \phi (\xi +\mu )}{(\xi +\mu +\omega )} \\ 0 &{} 0 &{} 0 &{} (1-e)(1-d)\theta \end{pmatrix}\nonumber \\ \end{aligned}$$

Their sum, $A^0+B^0$ and difference, $A^0-B^0$ are denoted by:

$$\begin{aligned} A^0+B^0=\begin{pmatrix} -(\mu +\omega ) &{} \xi &{} 0 &{} -\dfrac{(\beta +(1-d)\alpha \theta \phi )(\xi +\mu )}{(\xi +\mu +\omega )}\\ \omega &{} -(\xi +\mu ) &{} 0 &{} 0\\ 0 &{} 0 &{} -(\rho +\mu ) &{} \dfrac{(\beta +(1-d)\alpha \theta \phi )(\xi +\mu )}{(\xi +\mu +\omega )} \\ 0 &{} 0 &{} \rho &{} -(\epsilon +\gamma +\mu +\delta +(d+e(1-d))\theta ) \end{pmatrix}\nonumber \\ \end{aligned}$$

(E.9)

and

$$\begin{aligned}&A^0-B^0= \nonumber \\&\begin{pmatrix} -(2\alpha +\mu +\omega ) &{} \xi &{} 0 &{} -\dfrac{(\beta -(1-d)\alpha \theta \phi )(\xi +\mu )}{(\xi +\mu +\omega )} \\ \omega &{} -(2\alpha +\xi +\mu ) &{} 0 &{} 0\\ 0 &{} 0 &{} -(2\alpha +\rho +\mu ) &{} \dfrac{(\beta -(1-d)\alpha \theta \phi )(\xi +\mu )}{(\xi +\mu +\omega )} \\ 0 &{} 0 &{} \rho &{} -(\epsilon +\gamma +\mu +\delta +(2-(d+e(1-d)))\theta ) \end{pmatrix}\nonumber \\ \end{aligned}$$

(E.10)

The eigenvalues of $A^0+B^0$ and $A^0-B^0$ are found to be the solutions of the following characteristic polynomial equations, respectively:

$$\begin{aligned} (\lambda +\mu )(\lambda +(\mu +\xi +\omega ))(\lambda ^2+m\lambda +n)=0 \end{aligned}$$

and

$$\begin{aligned} (\lambda +(2\alpha +\mu ))(\lambda +(2\alpha +\mu +\xi +\omega ))(\lambda ^2+p\lambda +q)=0 \end{aligned}$$

where,

$$\begin{aligned} m&=(\epsilon +\rho +\gamma +2\mu +\delta +(d+e(1-d))\theta )\quad \\ n&=(\mu +\rho )(\epsilon +\gamma +\mu +\delta +(d+e(1-d))\theta )(1-R_{ed}^\phi )\quad \\ p&=(2(\alpha +\mu )+\rho +\epsilon +\gamma +\delta +(2-(d+e(1-d)))\theta )\quad \\ q&=(2\alpha +\mu +\rho )(\epsilon +\gamma +\mu +\delta +(2-(d+e(1-d)))\theta )\\&\left\{ 1-\dfrac{(\mu +\xi )(\beta -(1-d)\alpha \theta \phi )\rho }{(\mu +\xi +\omega )(2\alpha +\mu +\rho )(\epsilon +\gamma +\mu +\delta +(2-(d+e(1-d)))\theta )}\right\} \end{aligned}$$

since $0\le e, d\le 1 \Longrightarrow 0\le d+e(1-d) \le 1, and \; 2-(d+e(1-d))\ge (d+e(1-d))$, it is clear to see that m, $p>0$ for all non-negative parameters and $n>0$ for $R_{ed}^\phi <1$.

It can also show that $q>0$, since

$$\begin{aligned}&\dfrac{(\mu +\xi )(\beta -(1-d)\alpha \theta \phi )\rho }{(\mu +\xi +\omega )(2\alpha +\mu +\rho )(\epsilon +\gamma +\mu +\delta +(2-(d+e(1-d)))\theta )}<R_{ed}^\phi ,\\&\quad R_{ed}^\phi \text { is defined in 3.21} \end{aligned}$$

$$\begin{aligned}&\Longrightarrow \left\{ 1-\dfrac{(\mu +\xi )(\beta -(1-d)\alpha \theta \phi )\rho }{(\mu +\xi +\omega )(2\alpha +\mu +\rho )(\epsilon +\gamma +\mu +\delta +(2-(d+e(1-d)))\theta )}\right\} \\&\qquad>(1-R_{ed}^\phi )\\ \\&\Longrightarrow q>(2\alpha +\mu +\rho )(\epsilon +\gamma +\mu +\delta +(2-(d+e(1-d)))\theta )(1-R_{ed}^\phi ),\quad \nonumber \\&\qquad for\quad R_{ed}^\phi <1 \end{aligned}$$

This shows that all eigenvalues of the matrices $A+B$ and $A-B$ are negative real or have negative real parts. It follows that the disease-free equilibrium point is locally stable for $R_{ed}^\phi <1$ and it is unstable when $R_{ed}^\phi >1$ with respect to the full model system (3.13)–(3.16).

Proof of Theorem 6

The model in this computation is assumed to be the full model with travel between two patches, and the same biological parameters are used in each patch. Set the state variables of the model (3.13)–(3.16) as $S_1=x_1$, $W_1=x_2$, $E_1=x_3$, $I_1=x_4$, $S_2=x_5$, $W_2=x_6$, $E_2=x_7$, $I_2=x_8$ and we can rewrite the system as follows:

$$\begin{aligned} \frac{dx_1}{dt}&=\varLambda -\dfrac{\beta x_1x_4}{N_1}-(\omega +\mu )x_1+\xi x_2 -\alpha x_1 +\alpha x_5 -\dfrac{(1-d)\phi \alpha \theta x_5x_8}{N_2}:=f_1 \end{aligned}$$

(F.1)

$$\begin{aligned} \frac{dx_2}{dt}&=\omega x_1 -(\xi + \mu )x_2 -\alpha x_2 +\alpha x_6:=f_2 \end{aligned}$$

(F.2)

$$\begin{aligned} \frac{dx_3}{dt}&=\dfrac{\beta x_1x_4}{N_1}-(\rho +\mu ) x_3-\alpha x_3+\alpha x_7 +\dfrac{(1-d) \phi \alpha \theta x_5x_8}{N_2}:=f_3 \end{aligned}$$

(F.3)

$$\begin{aligned} \frac{dx_4}{dt}&=\rho x_3 -(\epsilon +\mu +\delta +\gamma +\theta )x_4 + (1-e)(1-d)\theta x_8:=f_4 \end{aligned}$$

(F.4)

$$\begin{aligned} \frac{dx_5}{dt}&=\varLambda -\dfrac{\beta x_5x_8}{N_2}-(\omega +\mu )x_5+\xi x_6 -\alpha x_5 +\alpha x_1 -\dfrac{(1-d)\phi \alpha \theta x_1x_4}{N_1}:=f_5 \end{aligned}$$

(F.5)

$$\begin{aligned} \frac{dx_6}{dt}&=\omega x_5 -(\xi + \mu )x_6 -\alpha x_6 +\alpha x_2:=f_6 \end{aligned}$$

(F.6)

$$\begin{aligned} \frac{dx_7}{dt}&=\dfrac{\beta x_5x_8}{N_2}- (\rho +\mu ) x_7-\alpha x_7+\alpha x_3 + \dfrac{(1-d)\phi \alpha \theta x_1x_4}{N_1}:=f_7 \end{aligned}$$

(F.7)

$$\begin{aligned} \frac{dx_8}{dt}&=\rho x_7 -(\epsilon +\mu +\delta +\gamma +\theta )x_8 + (1-e)(1-d)\theta x_4:=f_8 \end{aligned}$$

(F.8)

Let’s consider $\beta$ is a bifurcation parameter, then for $R_{ed}^\phi =1$, we have

$$\begin{aligned} \beta ^*=\dfrac{(\mu +\xi +\omega )(\mu +\rho )(\epsilon +\gamma +\mu +\delta +(d+e(1-d))\theta )}{(\mu +\xi )\rho }-(1-d)\alpha \theta \phi \end{aligned}$$

and the disease free equilibrium point is $X^0=\left( x^0_1, x^0_2, x^0_3, x^0_4, x^0_5, x^0_6, x^0_7, x^0_8\right)$ where,

$$\begin{aligned} x^0_1=x^0_5=\dfrac{\varLambda (\xi +\mu )}{\mu (\xi +\mu +\omega )},\quad x^0_2=x^0_6=\dfrac{\varLambda \omega }{\mu (\xi +\mu +\omega )},\quad x^0_3=x^0_4=x^0_7=x^0_8=0. \end{aligned}$$

Then the Jacobian matrix of the system F.1–F.8 around the disease free equilibrium ($X^0$) and $\beta =\beta ^*$ is given by:

$$\begin{aligned}&D_xf(X^0,\beta ^*)\nonumber \\&\quad = \begin{pmatrix} -Q_1 &{} \xi &{} 0 &{} -\beta ^*Q_5&{}\alpha &{} 0 &{} 0 &{} -(1-d)\alpha \theta \phi Q_5\\ \omega &{} -Q_2 &{} 0 &{} 0&{}0 &{} \alpha &{} 0 &{} 0 \\ 0 &{} 0 &{} -Q_3 &{} \beta ^*Q_5&{}0 &{} 0 &{} \alpha &{} (1-d)\alpha \theta \phi Q_5\\ 0 &{} 0 &{} \rho &{} -Q_4&{}0 &{} 0 &{} 0 &{} (1-e)(1-d)\theta \\ \alpha &{} 0 &{} 0 &{} -(1-d)\alpha \theta \phi Q_5 &{}-Q_1 &{} \xi &{} 0 &{} -\beta ^*Q_5 \\ 0 &{} \alpha &{} 0 &{} 0&{}\omega &{} -Q_2 &{} 0 &{} 0\\ 0 &{} 0 &{} \alpha &{} (1-d)\alpha \theta \phi Q_5&{} 0 &{} 0 &{} -Q_3 &{} \beta ^*Q_5 \\ 0 &{} 0 &{} 0 &{} (1-e)(1-d)\theta &{}0 &{} 0 &{} \rho &{} -Q_4 \end{pmatrix}\nonumber \\ \end{aligned}$$

(F.9)

Where,

$$\begin{aligned}&Q_1=\alpha +\mu +\omega , Q_2=\alpha +\xi +\mu , Q_3=\alpha +\mu +\rho \\&Q_4=\epsilon +\gamma +\mu +\delta +\theta , Q_5=\dfrac{\xi +\mu }{\xi +\mu +\omega } \end{aligned}$$

The right (w) and left (v) eigenvectors of the matrix in F.9 with respect to the zero eigenvalue satisfying $w.v=1$ are given by $w=(w_1, w_2, w_3, w_4, w_5, w_6, w_7, w_8)$ and $v=(v_1, v_2, v_3, v_4, v_5, v_6, v_7, v_8)$ respectively. Where,

$$\begin{aligned} w_1=w_5=\dfrac{-(\mu +\xi )(\mu +\rho )(\epsilon +\gamma +\mu +\delta +(d+e(1-d))\theta )}{\mu (\mu +\xi +\omega )},\\ w_2=w_6=\dfrac{-\omega (\mu +\rho )(\epsilon +\gamma +\mu +\delta +(d+e(1-d))\theta )}{\mu (\mu +\xi +\omega )},\\ w_3=w_7=(\epsilon +\gamma +\mu +\delta +(d+e(1-d))\theta ),\\ w_4=w_8=\rho \end{aligned}$$

and

$$\begin{aligned} v_1=v_2=v_5=v_6=0,\\ v_3=v_7=\dfrac{1}{2(\rho +\epsilon +\gamma +2\mu +\delta +(d+e(1-d))\theta )},\\v_4=v_8=\dfrac{(\mu +\rho )}{2\rho (\rho +\epsilon +\gamma +2\mu +\delta +(d+e(1-d))\theta )} \end{aligned}$$

Substituting those eigenvectors given above in a and b, that are

$$\begin{aligned} a&=\sum _{k,i,j=1}^8v_kw_iw_j\dfrac{\partial ^2f_k}{\partial x_i\partial x_j}(X^0, \beta ^*) \end{aligned}$$

(F.10)

$$\begin{aligned} b&=\sum _{k,i=1}^8v_kw_i\dfrac{\partial ^2f_k}{\partial x_i\partial \beta }(X^0, \beta ^*) \end{aligned}$$

(F.11)

where $f_k$ and $x_i$, $i,k=1,2, . . ., 8$ are corresponding equations and state variables of the system F.1–F.8 respectively, after some direct computation and simplification of equations F.10 and F.11 we get:

$$\begin{aligned}&a=-\dfrac{2\mu \omega (\mu +\rho )^2(\epsilon +\gamma +\mu +\delta +(d+e(1-d))\theta )^2}{\varLambda (\rho +\epsilon +\gamma +2\mu +\delta +(d+e(1-d))\theta )}\text {and}\\&b=\dfrac{\rho (\mu +\xi )}{(\xi +\mu +\omega )(\rho +\epsilon +\gamma +2\mu +\delta +(d+e(1-d))\theta )} \end{aligned}$$

Now, it shows that $a<0$ and $b>0$ for $\beta =\beta ^*$ at $R_{ed}^\phi =1$. It follows, by central manifold theory in Appendix B , the forward bifurcation appears and the endemic equilibrium point is locally asymptotically stable near to $R_{ed}^\phi =1$.