Time-Delay Dynamic Model and Cost-Effectiveness Analysis of Major Emergent Infectious Diseases with Transportation-Related Infections and Entry-Exit Screening

1. Introduction

In today’s society, infectious diseases are constantly threatening human health. Therefore, understanding the outbreak mechanism and transmission mode of infectious diseases will help shed light on future prevention and treatment strategies. When the first case of AIDS was covered in America in June 1981, more and more people were infected and even dying of AIDS throughout the 1980s, and there is still a large population of AIDS patients today [1]. In 2003, an outbreak of Severe Acute Respiratory Syndrome (SARS) occurred in Guangdong Province, China, characterized by its high infectivity and rapid disease progression, eventually leading to outbreaks in various regions of China and some cities worldwide [2]. In 2009, a widespread outbreak of H1N1 influenza occurred in the United States and rapidly spread to other countries and regions [3]. In March 2013, H7N9 was first discovered in eastern China [4]. A flare-up of H7N9 avian influenza quickly gained wide attention as new cases and high infection-related death rates were not controlled. At the end of 2019, COVID-19 rapidly swept the world, becoming one of the most serious diseases in human history [5]. An increasing number of sudden infectious diseases are emerging, some of which have become endemic, continuing to harm humans. In contrast, others have evolved into more insidious and harmful viruses through the continuous mutation of viruses.

Infectious disease modeling is a vital field that investigates the transmission routes and factors influencing the spread of emerging infectious diseases. With the rapid development of the global economy and the increasingly convenient means of transportation, population movements and people spreading the virus along the way are recognized as critical factors in disease transmission, which has attracted extensive attention from researchers. Perrin et al. [6] emphasized that the global dissemination of the HIV-1 pandemic fundamentally revolves around travel, implying the significance of travel in disease transmission. Air travel is also considered pivotal in spreading infectious diseases, as discussed by Mangili and Gendreau [7], who explored the potential for disease epidemics associated with air travel. In Beijing, China, Jia et al. [8] conducted a study revealing the impact of population mobility on tuberculosis’s prevalence, further highlighting the role of population mobility in disease transmission. In addition, the study of Chan et al. [9] revealed that person-to-person transmission of the novel coronavirus is possible within households, hospitals, and between cities. To better understand and predict transportation-related disease spread, Cui et al. [10] proposed an epidemiological model that describes this transmission mode. Similarly, Liu et al. [11] developed an SIQS infectious disease model incorporating transportation and entry screening, demonstrating the effectiveness of entry screening in mitigating transportation-related disease transmission. Many similar models of infectious diseases spread between two regions [12,13,14,15,16]. In addition to considering patch models for two regions, many authors also consider patch models for multiple regions. In 2003, Arino et al. [17] proposed a model for the spread of disease in a population of individuals traveling between n cities. The results show that ${\mathcal{R}}_0$ is a threshold, when ${\mathcal{R}}_0<1$, the disease disappears, and when ${\mathcal{R}}_0>1$, the disease reaches epidemic levels in all connected cities. Gao and Ruan [18] proposed a model with multiple patches to study how population movement affects the spread of malaria. Their results showed that travel between areas contributed to the spread of the disease in both regions. However, if travel rates continue to increase, the disease might disappear again in both areas. Sun et al. [19] established a patch model reflecting the population flow between Hubei and other regions. They estimate the epidemic situation in Hubei based on the daily reported epidemic data from Hubei and other regions, as well as the data on the population flow between Hubei and other regions. Zhang et al. [20] established a multi-patch model of HIV/AIDS with heterosexual transmission. The results indicate that if the disease vanishes in one region and spreads in another, it can spread or disappear in both regions depending on the individual mobility patterns.

Prosper et al. [21] proposed an optimal control model for malaria and find that the effectiveness of vaccination programs at the population level can be improved by actively seeking out and treating asymptomatic infections. Kang et al. [22] proposed a delayed avian influenza model incorporating slaughtering conditions. Their study indicates that optimal control, achieved by slaughtering susceptible (infected) avians and educating susceptible populations, can effectively minimize the number of infections in poultries and mankind. Moreover, the cost of implementing these control strategies can be minimized. Song et al. [23] studied the impact of delayed vaccination and isolation on COVID-19 transmission and conclude that the best isolation rate minimizes the total number of infections and expenditure on disease control. Singh et al. [24] have found that measures such as social distancing, lockdowns, and wearing masks can reduce the spread of the disease. These studies and models provide valuable insights and guidance for our in-depth understanding of the impact of population movements and traffic factors on disease transmission and how to control disease effectively.

For diseases such as tuberculosis, hepatitis A, hepatitis B, and influenza, the process of infection is unique because susceptible individuals do not immediately become ill after infection. Instead, they undergo an incubative period, where the pathogens replicate and spread within the host without showing noticeable symptoms. With this factor in mind, many researchers have proposed infectious disease models with delayed incubative periods in mathematical modeling and conducted extensive studies in this field [25,26,27,28,29,30,31,32]. Although many infectious disease models take into account the presence of incubative periods, we also need to be aware of the limitations of these models. The exact length of the incubative periods, variations in infectiousness, and interactions between individuals can all introduce uncertainty into the models. Therefore, further research is still needed to validate and improve the accuracy of these models by incorporating actual data and conducting field investigations. So it is also extremely interesting that we discuss the effect of time delay on measures of time-varying control.

Major emerging infectious diseases differ from normalized infectious diseases, as they have uncertainty, stronger stealthiness, more destructive power, and more severe harm to society and the economy. Therefore, studying the characteristics of sudden outbreaks and epidemics of infectious diseases can help reduce the harm of infectious diseases. From the above revelation, relatively limited studies use actual data to explore transportation-related infections between two regions and how to control this spread effectively. The rest of this paper is organized as follows. In the next section, we establish an SIQR time-delayed infectious disease model for emerging infectious diseases that considers transportation-related infection and entry-exit screening. We count the basic reproduction number and demonstrate that the endemic equilibrium is sole, the global asymptotic stability of the disease-free equilibrium, and the persistence of the model (1). In the section 3, we simulate the model using actual data from two regions in the USA and sensitivity analysis of ${\mathcal{R}}_0$. We investigate optimal control strategies in the section 4, and the cost-effectiveness analysis results are presented in the section 5. Discussions and conclusions are provided in the last two sections.

2. Time-Delay Model with Entry-Exit Screening and Transportation-Related Infections

2.1. Model Formulation

We posit that the two regions are identical, signifying that the population parameters in each region are comparable. This assumption holds merit when both regions exhibit similar population density levels, living conditions, economic factors, healthcare services, and disease transmission probabilities. We designate ${\theta }_e$ and ${\theta }_d$, where $0<{\theta }_e,{\theta }_d<1$, which represent the probabilities of successfully detecting an infected person for entry screening and exit screening, respectively. Moreover, we assume that susceptible individuals will not be mistakenly identified as infected. Once infected individuals are detected, they undergo isolation in designated hospitals for treatment. However, border screening differs from quarantine or isolation, as it explicitly addresses infections related to transportation. The former can thwart the spread of diseases from one region to another, while the latter effectively prevents the localized spread of diseases within a region. Nevertheless, in the absence of transportation-related infections or under extreme scenarios where ${\theta }_d=1$ or ${\theta }_e=1$, border screening equates to isolating or quarantining a certain proportion of infected individuals. We propose a hypothesis similar to that of many other authors [10,11,33].

After the patients recover, they have a temporary immunity to ensure they do not relapse for a certain period. We consider the disease’s incubation period, leading to the following time-delayed infectious disease model. In the model, $\tau$ stands for the disease’s incubation period, and $e^{-\mu \tau }{S}_i(t-\tau )I_i(t-\tau )$, $(i=1,2)$ denotes how many people have not died at time t after been infected at time $t-\tau$. $e^{-\mu \tau }$ is the probability that the infected human survives to time t (with the death rate $\mu$). We consider that $e^{-\mu \tau }$ approaches 1, $1-e^{-\mu \tau }$ approaches 0, so we ignore $1-e^{-\mu \tau }$ in our model. The model appears as follows:

$$\left\{\begin{array}{cc}\hfill \frac{d{S}_1}{dt}=& \mathsf{\Lambda }-\beta S_1{I}_1-\alpha S_1+\alpha S_2-\left(1-{\theta }_d\right)\gamma \alpha S_2{I}_2-\mu S_1,\hfill \\ \hfill \frac{d{I}_1}{dt}=& \beta S_1(t-\tau )I_1(t-\tau )e^{-\mu \tau }+\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\gamma \alpha S_2(t-\tau )I_2(t-\tau )e^{-\mu \tau }\hfill \\ & +\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\alpha I_2-(b+d+\alpha )I_1,\hfill \\ \hfill \frac{d{Q}_1}{dt}=& {\theta }_e\left(1-{\theta }_d\right)\gamma \alpha S_2{I}_2+{\theta }_d\alpha I_1+{\theta }_e\left(1-{\theta }_d\right)\alpha I_2-(c+f)Q_1,\hfill \\ \hfill \frac{d{R}_1}{dt}=& f{Q}_1+d{I}_1-\alpha R_1+\alpha R_2-\mu R_1,\hfill \\ \hfill \frac{d{S}_2}{dt}=& \mathsf{\Lambda }-\beta S_2{I}_2-\alpha S_2+\alpha S_1-\left(1-{\theta }_d\right)\gamma \alpha S_1{I}_1-\mu S_2,\hfill \\ \hfill \frac{d{I}_2}{dt}=& \beta S_2(t-\tau )I_2(t-\tau )e^{-\mu \tau }+\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\gamma \alpha S_1(t-\tau )I_1(t-\tau )e^{-\mu \tau }\hfill \\ & +\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\alpha I_1-(b+d+\alpha )I_2,\hfill \\ \hfill \frac{d{Q}_2}{dt}=& {\theta }_e\left(1-{\theta }_d\right)\gamma \alpha S_1{I}_1+{\theta }_d\alpha I_2+{\theta }_e\left(1-{\theta }_d\right)\alpha I_1-(c+f)Q_2,\hfill \\ \hfill \frac{d{R}_2}{dt}=& f{Q}_2+d{I}_2-\alpha R_2+\alpha R_1-\mu R_2.\hfill \end{array}\right.$$

(1)

Here $S_i$ and $I_i$$(i=1,2)$ represent the number of susceptible and infected individuals in region i (we use the results of the entry-exit screening as a criterion for whether the population has the disease), $Q_i$$(i=1,2)$ denotes such isolated infected individuals in each region, $R_i$$(i=1,2)$ represents the number of the removed individuals in region i. Table 1 lists the biological interpretations of all parameters. The per capita mortality rate for infected individuals is b, and c means the average death rate of isolated infected individuals. Because not only natural mortality but disease-caused mortality are included in it, we have $\mu \le b$, $\mu \le c$. In model (1), rate of disease (that is, how many new infection cases are added within a unit of time) $\beta S_i{I}_i$, $i=1,2$, within city i. The susceptible and infected individuals in each i leave for city j ($i\ne j,i,j=1,2$) at a per capita rate of $\alpha$. When the individuals in city i travel to the city j, disease is transmitted with the incidence rate $\gamma \alpha S_i{I}_i$, $i=1,2$, and the variable $\gamma$ represents the transmission rate associated with transportation.

From the biological perspective, the susceptible count during travel should not drop below zero, for all $S_i,I_i\ge 0$ ($i=1,2$).

$$\alpha S_i-\left(1-{\theta }_d\right)\gamma \alpha S_i{I}_i>0,(i=1,2).$$

The form of the initial conditions for model (1) is

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{S}_1\left(\delta \right)={\varphi }_1\left(\delta \right),I_1\left(\delta \right)={\varphi }_2\left(\delta \right),Q_1\left(\delta \right)={\varphi }_3\left(\delta \right),R_1\left(\delta \right)={\varphi }_4\left(\delta \right),\hfill \\ S_2\left(\delta \right)={\varphi }_5\left(\delta \right),I_2\left(\delta \right)={\varphi }_6\left(\delta \right),Q_2\left(\delta \right)={\varphi }_7\left(\delta \right),R_2\left(\delta \right)={\varphi }_8\left(\delta \right),\hfill \\ {\varphi }_i\left(\delta \right)\ge 0,\delta \in [-\tau ,0],{\varphi }_i\left(0\right)>0,{\varphi }_i\left(\delta \right)\in \mathbb{C}\left([-\tau ,0],{\mathbb{R}}_{+}^8\right),i=1,2,\cdots ,8.\hfill \end{array}\right.\end{array}$$

(2)

Figure 1. Major emerging infectious disease SIQR model based on infection associated with transportation-related infection and entry-exit screening.

Figure 1. Major emerging infectious disease SIQR model based on infection associated with transportation-related infection and entry-exit screening.

2.2. Stability Analysis

The disease-free equilibrium (DFE) in model (1) is denoted as

$$P^0=\left(S^0,0,0,0,S^0,0,0,0\right),S^0=\frac{\mathsf{\Lambda }}{\mu }.$$

The next generation matrix of the model (1) is $F{V}^{-1}$. According to the articles [34,35], the basic reproduction number of model (1) we get is as follows:

$$\begin{array}{c}\hfill {\mathcal{R}}_0=\rho \left(F{V}^{-1}\right)=\frac{(\beta +\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\gamma \alpha )e^{-\mu \tau }{S}^0}{b+d+{\theta }_d\alpha +{\theta }_e\left(1-{\theta }_d\right)\alpha }.\end{array}$$

The detailed calculation process of ${\mathcal{R}}_0$ is shown in .

Theorem 1. 

If ${\mathcal{R}}_0>1$, model (1) have a unique endemic equilibrium

$$P^{*}=\left(S_1,I_1,Q_1,R_1,S_2,I_2,Q_2,R_2\right).$$

Proof 

We use the method proposed by Wang et al. [12] to prove this theorem. We set $N_1=S_1+I_1$ and $N_2=S_2+I_2$. By equating the right side of the second and sixth equations in model (1) to zero, the following two equalities are obtained:

$$\begin{array}{c}\hfill \beta e^{-\mu \tau }{I}_1^2+\left[b+d+\alpha -\beta e^{-\mu \tau }{N}_1\right]I_1+\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\gamma \alpha e^{-\mu \tau }{I}_2^2\\ \hfill -\left[\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\gamma \alpha e^{-\mu \tau }{N}_2+\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\alpha \right]I_2=0,\end{array}$$

(3)

$$\begin{array}{c}\hfill \left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\gamma \alpha e^{-\mu \tau }{I}_1^2-\left[\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\gamma \alpha e^{-\mu \tau }{N}_2+\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\alpha \right]I_1\\ \hfill +\beta e^{-\mu \tau }{I}_2^2+\left[b+d+\alpha -\beta e^{-\mu \tau }{N}_1\right]I_2=0.\end{array}$$

(4)

Eq. (3) and (4) can be written as two ellipses,

$$\begin{array}{cc}\hfill & \frac{{\left(I_1-x_1\right)}^2}{a_1^2}+\frac{{\left(I_2-y_1\right)}^2}{b_1^2}=1,\hfill \\ \hfill & \frac{{\left(I_1-x_2\right)}^2}{a_2^2}+\frac{{\left(I_2-y_2\right)}^2}{b_2^2}=1.\hfill \end{array}$$

(5)

Where,

$$\begin{array}{cc}\hfill & x_1=\frac{1}{2}{N}_1\left(1-\frac{b+d+\alpha }{N_1\gamma e^{-\mu \tau }}\right),y_1=\frac{1}{2}{N}_2\left(1+\frac{1}{N_2\gamma e^{-\mu \tau }}\right),\hfill \\ \\ \hfill & x_2=\frac{1}{2}{N}_1\left(1+\frac{1}{N_1\gamma e^{-\mu \tau }}\right),y_2=\frac{1}{2}{N}_2\left(1-\frac{b+d+\alpha }{N_2\gamma e^{-\mu \tau }}\right),\hfill \end{array}$$

$$\begin{array}{ccc}\hfill & \begin{array}{cc}\hfill & a_1^2=\frac{\beta e^{-\mu \tau }{x}_1^2+\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\gamma \alpha e^{-\mu \tau }{y}_1^2}{\beta e^{-\mu \tau }},\hfill \\ \\ \hfill & b_1^2=\frac{\beta e^{-\mu \tau }{x}_1^2+\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\gamma \alpha e^{-\mu \tau }{y}_1^2}{\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\gamma \alpha e^{-\mu \tau }},\hfill \end{array}\hfill & \hfill \begin{array}{cc}\hfill & a_2^2=\frac{\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\gamma \alpha e^{-\mu \tau }{x}_2^2+\beta e^{-\mu \tau }{y}_2^2}{\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\gamma \alpha e^{-\mu \tau }},\hfill \\ \\ \hfill & b_2^2=\frac{\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\gamma \alpha e^{-\mu \tau }{x}_2^2+\beta e^{-\mu \tau }{y}_2^2}{\beta e^{-\mu \tau }}.\hfill \end{array}\end{array}$$

We know these two ellipses must pass through the origin, where $x_2>0$, $y_1>0$, while the positive or negative cases of $x_1$ and $y_2$ remain uncertain. Therefore, we consider $I_1$ as the horizontal axis and $I_2$ as the vertical axis to determine the intersection point of the two ellipses in the Cartesian coordinate system, as shown in Figure 2. From Figure 2, we observe that in the first quadrant, the two ellipses have one and only one intersection point. Thus, system (1) exists a unique positive equilibrium $P^{*}$ when ${\mathcal{R}}_0>1$. Of course, we cannot express the specific expression of $P^{*}$. There are also special positive equilibriums, such as the one proposed by Liu [33].

□

Theorem 2. 

If ${\mathcal{R}}_0<1$, the disease-free equilibrium $P^0$ of model (1) is locally asymptotically stable.

Proof 

The Jacobian matrix of model (1) at $P^0$ is:

$$J=\left(\begin{array}{cccc}A& O& B& O\\ C& M& D& N\\ B& O& A& O\\ D& N& C& M\end{array}\right),$$

Where, O is the zero matrix,

$A=\left(\begin{array}{cc}-\alpha -\mu & -\beta S^0\\ 0& \beta S^0{e}^{-\lambda \tau }{e}^{-\mu \tau }-\left(b+d+\alpha \right)\end{array}\right),$

$B=\left(\begin{array}{cc}\alpha & -\left(1-{\theta }_d\right)\gamma \alpha S^0\\ 0& \left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\gamma \alpha S_0{e}^{-\lambda \tau }{e}^{-\mu \tau }+\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\alpha \end{array}\right),$

$C=\left(\begin{array}{cc}0& {\theta }_d\alpha \\ 0& d\end{array}\right)$, $D=\left(\begin{array}{cc}0& {\theta }_e(1-{\theta }_d)\gamma \alpha S^0+{\theta }_e(1-{\theta }_d)\alpha \\ 0& 0\end{array}\right),$

$M=\left(\begin{array}{cc}-(c+f)& 0\\ f& -(\alpha +\mu )\end{array}\right),$$N=\left(\begin{array}{cc}0& 0\\ 0& -\alpha \end{array}\right).$

The equation of characteristics of model (1) at the disease-free equilibrium $P^0$ is

$\begin{array}{cc}\hfill J_1=& {(\lambda +\alpha +\mu )}^2{(\lambda +c+f)}^2(\lambda +\mu )(\lambda +2\alpha +\mu )[\lambda +d+b+\alpha -\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\alpha \hfill \\ & -\left(\beta +\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\gamma \alpha \right)S^0{e}^{-\mu \tau }{e}^{-\lambda \tau }\left]\right[\lambda +d+b+\alpha +\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\alpha \hfill \\ & -\left(\beta -\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\gamma \alpha \right)S^0{e}^{-\mu \tau }{e}^{-\lambda \tau }]=0.\hfill \end{array}$

The characteristic polynomial has negative roots ${\lambda }_1=-\mu$, ${\lambda }_2=-\left(2\alpha +\mu \right),{\lambda }_3=-(\alpha +\mu ),{\lambda }_4=-(c+f),$ and the remaining roots are obtained from the following:

$$\lambda +d+b+\alpha -\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\alpha =\left(\beta +\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\gamma \alpha \right)S^0{e}^{-\mu \tau }{e}^{-\lambda \tau },$$

(6)

$$\lambda +d+b+\alpha +\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\alpha =\left(\beta -\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\gamma \alpha \right)S^0{e}^{-\mu \tau }{e}^{-\lambda \tau }.$$

(7)

We assume that $\lambda$ has a nonnegative real part, modulo both sides of the above two equations (6) and (7), and the left side of the equations is

$$\begin{array}{c}\left|\lambda +d+b+\alpha -\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\alpha \right|>d+b+{\theta }_d\alpha +{\theta }_e\left(1-{\theta }_d\right)\alpha ,\\ \left|\lambda +d+b+\alpha +\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\alpha \right|>d+b+\alpha +\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\alpha >d+b+{\theta }_d\alpha +{\theta }_e\left(1-{\theta }_d\right)\alpha .\end{array}$$

The right side of the equations is

$$\begin{array}{c}\left|\left(\beta +\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\gamma \alpha \right)S^0{e}^{-\mu \tau }{e}^{-\lambda \tau }\right|<\left(\beta +\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\gamma \alpha \right)S^0{e}^{-\mu \tau },\\ \left|\left(\beta -\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\gamma \alpha \right)S^0{e}^{-\mu \tau }{e}^{-\lambda \tau }\right|<\left(\beta +\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\gamma \alpha \right)S^0{e}^{-\mu \tau }.\end{array}$$

When ${\mathcal{R}}_0<1$,

$$\left(\beta +\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\gamma \alpha \right)S^0{e}^{-\mu \tau }{e}^{-\lambda \tau }<b+d+{\theta }_d\alpha +{\theta }_e\left(1-{\theta }_d\right)\alpha .$$

When ${\mathcal{R}}_0<1$, $\lambda =0$ is not the root of model (1), so $\lambda$ does not have the nonnegative real part. This contradicts the assumption that $\lambda$ has negative real part, i.e., the disease-free equilibrium $P^0$ of model (1) is locally asymptotically stable. □

Theorem 3. 

If ${\mathcal{R}}_0<1$, the disease-free equilibrium $P^0$ of model (1) is globally asymptotically stable.

Proof 

Let us take the following Lyapunov functional into account

$$\begin{array}{cc}\hfill V\left(t\right)=& I_1+I_2+{\int }_{t-\tau }^t\beta S_1\left(\varphi \right)I_1\left(\varphi \right)e^{-\mu \tau }+\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\gamma \alpha S_2\left(\varphi \right)I_2\left(\varphi \right)e^{-\mu \tau }d\varphi \hfill \\ & +{\int }_{t-\tau }^t\beta S_2\left(\varphi \right)I_2\left(\varphi \right)e^{-\mu \tau }+\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\gamma \alpha S_1\left(\varphi \right)I_1\left(\varphi \right)e^{-\mu \tau }d\varphi .\hfill \end{array}$$

The time derivative of it along the solutions of model (1) is

$$\begin{array}{cc}\hfill \dot{V\left(t\right)}& =\left(\beta +\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\gamma \alpha \right)e^{-\mu \tau }\left(S_1{I}_1+S_2{I}_2\right)-\left(b+d+{\theta }_d\alpha +{\theta }_e\left(1-{\theta }_d\right)\alpha \right)\left(I_1+I_2\right)\hfill \\ \hfill & \le \left(\beta +\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\gamma \alpha \right)\frac{\mathsf{\Lambda }}{\mu }{e}^{-\mu \tau }\left(I_1+I_2\right)-\left(b+d+{\theta }_d\alpha +{\theta }_e\left(1-{\theta }_d\right)\alpha \right)\left(I_1+I_2\right)\hfill \\ \hfill & =\left(b+d+{\theta }_d\alpha +{\theta }_e\left(1-{\theta }_d\right)\alpha \right)\left(R_0-1\right)\left(I_1+I_2\right).\hfill \end{array}$$

When ${\mathcal{R}}_0<1$,

$$\left(b+d+{\theta }_d\alpha +{\theta }_e\left(1-{\theta }_d\right)\alpha \right)\left(R_0-1\right)\left(I_1+I_2\right)<0.$$

In delay differential systems, by LaSalle’s invariance principle, the disease-free equilibrium $P^0$ of model (1) is globally asymptotically stable. □

Figure 2. The four possible cases of two ellipses given by (5) when ${\mathcal{R}}_0>1$.

Figure 2. The four possible cases of two ellipses given by (5) when ${\mathcal{R}}_0>1$.

3. Data Fitting and Sensitivity Analysis of ${\mathcal{R}}_0$

3.1. Data Fitting

In order to better evaluate the effectiveness of transportation-related infection and exit-entry screenings in preventing and controlling emerging infectious diseases, we will choose actual case data to illustrate. We take the COVID-19 epidemic in the United States as an example. On January 20, 2020, the first confirmed case of COVID-19 was reported in the US, Washington [36], so we take Washington as the hypothetical region 1 and select the surrounding three states of Idaho, Nevada, and Oregon as the hypothetical region 2.

In 2020, the average life expectancy of the U.S. population was $77.0$ years [37], approximately 300 new births occur in Washington every day [38]. Therefore, we choose the range of $\mathsf{\Lambda }$ to be $(0,300)$, and the natural death rate $\mu$ is $1/365\times 77.0=3.56\times {10}^{-5}$. We choose data ranging from June 10, 2020, to September 10, 2020. In Washington, the number of daily deaths due to illness per $100,000$ people ranges from approximately 100 to 300 [39]. Therefore, the range for the average death rate of the infected b is $(0.001,0.003)$. The estimated hospitalization death rate of isolated infected individuals ranges from 5% and 20% [40], so we choose the range for the average death rate of isolated infected individuals c to be $(0.05,0.2)$. Lightly symptomatic patients typically recover in 7 to 14 days, while severely symptomatic patients may take 1 to 2 months to recover [41]. We choose the same range for the recovery rate per capita during isolation treatment f and the per capita rate of natural recovery d, which is $(0.016,0.14)$. The incubative period of COVID-19 is typically between 2 and 14 days, so the range for the incubative period $\tau$ is $(2,14)$ [42]. According to the actual prevention and control situation in the United States at that time [43,44], we make some assumptions on the parameters of the model (1). We assume that 10% of the people in Washington, Idaho, Nevada, and Oregon flow daily, and the upper limit of the migration rate is $0.1\times 1=0.1$. Considering that all people leave Washington to go to Idaho, Nevada, and Oregon, or leave Idaho, Nevada, and Oregon to travel to Washington at least once a year with a lower limit of mobility is $1/365=0.0027$, the range of $\alpha$ is $(0.0027,0.1)$. Assuming that the success rate of entry and exit screening is less than half, the ${\theta }_e$ and ${\theta }_d$ ranges are $(0.1,0.5)$. We obtain the parameter range of Table 2.

Since the data we use is cumulative confirmed cases from the Population Health Institute of the Robert Wood Johnson Foundation at the University of Wisconsin [45], before conducting data fitting, we need to add two cumulative confirmed cases compartments represented by $C_i$$(i=1,2)$ :

$$\begin{array}{cc}\hfill \frac{d{C}_1}{dt}=& \beta S_1(t-\tau )I_1(t-\tau )e^{-\mu \tau }+\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\gamma \alpha S_2(t-\tau )I_2(t-\tau )e^{-\mu \tau }\hfill \\ & +\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\alpha I_2,\hfill \\ \hfill \frac{d{C}_2}{dt}=& \beta S_2(t-\tau )I_2(t-\tau )e^{-\mu \tau }+\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\gamma \alpha S_1(t-\tau )I_1(t-\tau )e^{-\mu \tau }\hfill \\ & +\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\alpha I_1.\hfill \end{array}$$

The initial values $S_1\left(0\right)=7\times {10}^6,S_2\left(0\right)=9\times {10}^6,I_1\left(0\right)=300,I_2\left(0\right)=200,Q_1\left(0\right)=0,Q_2\left(0\right)=0,R_1\left(0\right)=0,R_2\left(0\right)=0$ are selected. We use the method proposed by Zhang [46] and Li et al. [47,48] to carry out parameter estimation and data fitting, and the estimated optimal parameter values can be seen in Table 2.

Figure 3 displays the accumulated confirmed cases in Washington, Idaho, Nevada, and Oregon under the actual data. The simulation results are consistent with the actual data, proving that model (1) is more suitable to simulate the outbreak and epidemic of emerging infectious diseases.

3.2. Sensitivity Analysis of ${\mathcal{R}}_0$

Sensitivity analysis is performed to discern the degree of correlation between the input parameters and the predicted ones, as well as to figure out the necessary degree of variation in the input parameters to achieve the expected value of the parameters. Generally speaking, the beginning of disease transmission is directly affected by the basic reproduction number ${\mathcal{R}}_0$. In this segment, in order to study the effectiveness of COVID-19 control strategies in the United States, we calculate the sensitivity index of ${\mathcal{R}}_0$, whose purpose is to survey the relative alterations in ${\mathcal{R}}_0$ as the model parameters changes. Then, we ponder the parameters from Table 2. Sensitivity analysis based on PRCC evaluated the effect of parameters on the response function, i.e., basic reproduction number ${\mathcal{R}}_0$. We take the sample size $N=10000$, in addition to $\lambda$, $\mu$, all parameters as the input variables, the value of ${\mathcal{R}}_0$ as output variables. It can be obtained from Figure 4 which shows that $\beta$ the most sensitive parameter then followed d, $\alpha$, ${\theta }_d$, ${\theta }_e$, $\gamma$, $\tau$.

We propose adequate measures for preventing and controlling COVID-19 based on sensitivity analysis. These measures include wearing masks to interrupt transmission pathways, implementing city lockdowns to control population movement, and timely medication for patients who are ill but have not yet been hospitalized. Implementing these measures can enhance the effectiveness of prevention and control efforts. In the next part, we are going to explore the best control ways.

4. Optimal Control Strategies

Owing to the previous ${\mathcal{R}}_0$ sensitivity analysis and combined with the reality of COVID-19 in America, we put our hands to seek an effective way to control COVID-19 infections. $u_1\left(t\right)$ aims to control the transmission rate of susceptible people, which is measured by preventive measures like wearing masks and reducing gatherings. $u_2\left(t\right)$ aims to control population migration, measured through policy interventions and city lockdowns. $u_3\left(t\right)$ aims to take medication and actively pursue treatment to enhance the natural recovery rate. With all the assumptions we mentioned above, the time-delayed model takes the control problem as follows:

$$\left\{\begin{array}{cc}\hfill \frac{d{S}_1}{dt}=& \mathsf{\Lambda }-\beta \left(1-u_1\right)S_1{I}_1-\alpha \left(1-u_2\right)S_1+\alpha \left(1-u_2\right)S_2-\left(1-{\theta }_d\right)\gamma \alpha \left(1-u_2\right)S_2{I}_2-\mu S_1,\hfill \\ \hfill \frac{d{I}_1}{dt}=& \beta \left(1-u_1\right)S_1(t-\tau )I_1(t-\tau )e^{-\mu \tau }+\left(1-{\theta }_d\right)\left(1-{\theta }_e\right)\gamma \alpha \left(1-u_2\right)S_2(t-\tau )I_2(t-\tau )e^{-\mu \tau }\hfill \\ \hfill & +\left(1-{\theta }_d\right)\left(1-{\theta }_e\right)\alpha \left(1-u_2\right)I_2-\left(b+d\left(1+u_3\right)+\alpha \left(1-u_2\right)\right)I_1,\hfill \\ \hfill \frac{d{Q}_1}{dt}=& {\theta }_e\left(1-{\theta }_d\right)\gamma \alpha \left(1-u_2\right)S_2{I}_2+{\theta }_d\alpha \left(1-u_2\right)I_1+{\theta }_e\left(1-{\theta }_d\right)\alpha \left(1-u_2\right)I_2-(c+f)Q_1,\hfill \\ \hfill \frac{d{R}_1}{dt}=& f{Q}_1+d\left(1+u_3\right)I_1-\alpha \left(1-u_2\right)R_1+\alpha \left(1-u_2\right)R_2-\mu R_1,\hfill \\ \hfill \frac{d{S}_2}{dt}=& \mathsf{\Lambda }-\beta \left(1-u_1\right)S_2{I}_2-\alpha \left(1-u_2\right)S_2+\alpha \left(1-u_2\right)S_1-\left(1-{\theta }_d\right)\gamma \alpha \left(1-u_2\right)S_1{I}_1-\mu S_2,\hfill \\ \hfill \frac{d{I}_2}{dt}=& \beta \left(1-u_1\right)S_2(t-\tau )I_2(t-\tau )e^{-\mu \tau }+\left(1-{\theta }_d\right)\left(1-{\theta }_e\right)\gamma \alpha \left(1-u_2\right)S_1(t-\tau )I_1(t-\tau )e^{-\mu \tau }\hfill \\ \hfill & +\left(1-{\theta }_d\right)\left(1-{\theta }_e\right)\alpha \left(1-u_2\right)I_1-\left(b+d\left(1+u_3\right)+\alpha \left(1-u_2\right)\right)I_2,\hfill \\ \hfill \frac{d{Q}_2}{dt}=& {\theta }_e\left(1-{\theta }_d\right)\gamma \alpha \left(1-u_2\right)S_1{I}_1+{\theta }_d\alpha \left(1-u_2\right)I_2+{\theta }_e\left(1-{\theta }_d\right)\alpha \left(1-u_2\right)I_1-(c+f)Q_2,\hfill \\ \hfill \frac{d{R}_2}{dt}=& f{Q}_2+d\left(1+u_3\right)I_2-\alpha \left(1-u_2\right)R_2+\alpha \left(1-u_2\right)R_1-\mu R_2.\hfill \end{array}\right.$$

The initial conditions

$$\begin{array}{c}\hfill \begin{array}{c}\hfill \left\{\begin{array}{c}{S}_1\left(\delta \right)={\varphi }_1\left(\delta \right),I_1\left(\delta \right)={\varphi }_2\left(\delta \right),Q_1\left(\delta \right)={\varphi }_3\left(\delta \right),R_1\left(\delta \right)={\varphi }_4\left(\delta \right),\hfill \\ S_2\left(\delta \right)={\varphi }_5\left(\delta \right),I_2\left(\delta \right)={\varphi }_6\left(\delta \right),Q_2\left(\delta \right)={\varphi }_7\left(\delta \right),R_2\left(\delta \right)={\varphi }_8\left(\delta \right),\hfill \\ {\varphi }_i\left(\delta \right)\ge 0,\delta \in [-\tau ,0],{\varphi }_i\left(0\right)>0,{\varphi }_i\left(\delta \right)\in \mathbb{C}\left([-\tau ,0],{\mathbb{R}}_{+}\right),i=1,2,\cdots ,8.\hfill \end{array}\right.\end{array}\end{array}$$

Regarding variables $u_i\left(t\right)(i=1,2,3)$ and the above control system, considering that the spread of emerging infectious diseases is related to transportation-related infections between two regions, we consider reducing the number of patients who are not successfully detected in entry and exit screenings between the two regions.

Our objective function $J(u_1\left(t\right),u_2\left(t\right),u_3\left(t\right))$ is defined as

$$\begin{array}{cc}\hfill J(u_1\left(t\right),u_2\left(t\right),u_3\left(t\right))={\int }_0^T& [A_1(1-{\theta }_d)(1-{\theta }_e)\alpha I_1+A_2(1-{\theta }_d)(1-{\theta }_e)\alpha I_2\hfill \\ & +\frac{1}{2}(B_1{u}_1^2\left(t\right)+B_2{u}_2^2\left(t\right)+B_3{u}_3^2\left(t\right))]dt,\hfill \end{array}$$

where the interval $[0,T]$ signifies the time duration of the control program, the quantities $A_i(i=1,2)$ represents the weight constants of patients in region 1 and region 2 who were not successfully detected. $B_i(i=1,2,3)$ represents the transmission rate of the susceptible population. In every region i, the susceptible and infected individuals move to region j at a per capita rate $\alpha$ ($i,j=1,2,i\ne j$) and the per capita rate of natural recovery from disease. Constants $A_i,(i=1,2)$, $B_i,(i=1,2,3)$ signify a measure of the relative cost of the interventions over the interval $[0,T]$. The expression $(1-{\theta }_d)(1-{\theta }_e)\alpha I_i$, $(i=1,2)$ refers to the patients who are not successfully detected between two regions. Therefore, this term is included in the objective function. In addition, we should consider the costs of wearing masks, city lockdowns, and purchasing medications. Our goal is to find the optimal control $(u_1,u_2,u_3)$ such as

$$J\left(u_1^{*}\left(t\right),u_2^{*}\left(t\right),u_3^{*}\left(t\right)\right)=min\left\{J\left(u_1\left(t\right),u_2\left(t\right),u_3\left(t\right)\right)\mid \left(u_1\left(t\right),u_2\left(t\right),u_3\left(t\right)\right)\in U\right\},$$

the control set U equals that

$$U=\left\{u_1\left(t\right),u_2\left(t\right),u_3\left(t\right)\mid u_i\left(t\right)(i=1,2,3)\mathrm{can}\mathrm{be}\mathrm{measured},u_i\left(t\right)\in (0,1),t\in [0,T]\right\}.$$

Through the Filippov–Cesari Existence Theorem and the Mangasarian Theorem [49], a special method $u_1^{*}\left(t\right),u_2^{*}\left(t\right),u_3^{*}\left(t\right)$ which is the optimal control exists. After that, we can find this method in Pontryagin’s Minimum Principle [50].

Theorem 4. 

Suppose $u_1^{*}$, $u_2^{*}$ and $u_3^{*}$ are the optimal control variables, and correspondingly $S_1^{*}$, $I_1^{*}$, $Q_1^{*}$, $R_1^{*}$, $S_2^{*}$, $I_2^{*}$, $Q_2^{*}$, $R_2^{*}$ are the control system’s (8) optimal state variables, then there exists adjoint variable $\lambda \left(t\right)=({\lambda }_1\left(t\right),{\lambda }_2\left(t\right),{\lambda }_3\left(t\right),{\lambda }_4\left(t\right),{\lambda }_5\left(t\right),{\lambda }_6\left(t\right),{\lambda }_7\left(t\right),{\lambda }_8\left(t\right))\in {\mathbb{R}}_8$ that satisfies the following adjoint equations:

$$\begin{array}{cc}\hfill \frac{d{\lambda }_1\left(t\right)}{dt}=& {\lambda }_1\left(t\right)\left[\beta \left(1-u_1\right)I_1+\mu \right]+\left({\lambda }_5\left(t\right)-{\lambda }_7\left(t\right){\theta }_e\right)\left(1-{\theta }_d\right)\gamma \alpha \left(1-u_2\right)I_1\hfill \\ \hfill & -{\chi }_{[0,T-\tau ]}\left(t\right)\beta \left(1-u_1\right)I_1(t-\tau )e^{-\mu \tau }{\lambda }_2(t+\tau )+\left({\lambda }_1\left(t\right)-{\lambda }_5\left(t\right)\right)\alpha \left(1-u_2\right)\hfill \\ \hfill & -{\chi }_{[0,T-\tau ]}\left(t\right)\left(1-{\theta }_d\right)\left(1-{\theta }_e\right)\gamma \alpha \left(1-u_2\right)I_1(t-\tau )e^{-\mu \tau }{\lambda }_6(t+\tau ),\hfill \\ \hfill \frac{d{\lambda }_2\left(t\right)}{dt}=& -A_1(1-{\theta }_d)(1-{\theta }_e)\alpha +{\lambda }_1\left(t\right)\beta \left(1-u_1\right)S_1+{\lambda }_2\left(t\right)\left(b+d\left(1+u_3\right)+\alpha \left(1-u_2\right)\right)\hfill \\ \hfill & -{\lambda }_4\left(t\right)(1+d{u}_3)+{\lambda }_5\left(t\right)\left(1-{\theta }_d\right)\gamma \alpha (1-\left.u_2\right)S_1-{\lambda }_7\left(t\right){\theta }_e\left(1-{\theta }_d\right)\gamma \alpha \left(1-u_2\right)S_1\hfill \\ \hfill & -{\lambda }_6\left(t\right)\left(1-{\theta }_d\right)\left(1-{\theta }_e\right)\alpha \left(1-u_2\right)-{\lambda }_7\left(t\right){\theta }_e\left(1-{\theta }_d\right)\alpha (1-u_2)\hfill \\ \hfill & -{\chi }_{[0,T-\tau ]}\left(t\right)\beta \left(1-u_1\right)S_1(t-\tau )e^{-\mu \tau }{\lambda }_2(t+\tau )-{\lambda }_3\left(t\right){\theta }_d\alpha (1-u_2)\hfill \\ \hfill & -{\chi }_{[0,T-\tau ]}\left(t\right)(1-{\theta }_d)(1-{\theta }_e)\gamma \alpha (1-\left.u_2\right)S_1(t-\tau )e^{-\mu \tau }{\lambda }_6(t+\tau ),\hfill \\ \hfill \frac{d{\lambda }_3\left(t\right)}{dt}=& {\lambda }_3\left(t\right)(c+f)-{\lambda }_4\left(t\right)f,\hfill \\ \hfill \frac{d{\lambda }_4\left(t\right)}{dt}=& \left({\lambda }_4\left(t\right)-{\lambda }_8\left(t\right)\right)\alpha \left(1-u_2\right)+{\lambda }_4\left(t\right)\mu ,\hfill \\ \hfill \frac{d{\lambda }_5\left(t\right)}{dt}=& {\lambda }_5\left(t\right)\left[\beta \left(1-u_1\right)I_2+\mu \right]+\left({\lambda }_1\left(t\right)-{\lambda }_3\left(t\right){\theta }_e\right)\left(1-{\theta }_d\right)\gamma \alpha \left(1-u_2\right)I_2\hfill \\ \hfill & -{\chi }_{[0,T-\tau ]}\left(t\right)\beta \left(1-u_1\right)I_2(t-\tau )e^{-\mu \tau }{\lambda }_6(t+\tau )+\left({\lambda }_5\left(t\right)-{\lambda }_1\left(t\right)\right)\alpha \left(1-u_2\right)\hfill \\ \hfill & -{\chi }_{[0,T-\tau ]}\left(t\right)\left(1-{\theta }_d\right)\left(1-{\theta }_e\right)\gamma \alpha \left(1-u_2\right)I_2(t-\tau )e^{-\mu \tau }{\lambda }_2(t+\tau ),\hfill \\ \hfill \frac{d{\lambda }_6\left(t\right)}{dt}=& -A_2(1-{\theta }_d)(1-{\theta }_e)\alpha +{\lambda }_5\left(t\right)\beta \left(1-u_1\right)S_2+{\lambda }_6\left(t\right)\left(b+d\left(1+u_3\right)+\alpha \left(1-u_2\right)\right)\hfill \\ & -{\lambda }_8\left(t\right)(1+d{u}_3)+{\lambda }_1\left(t\right)\left(1-{\theta }_d\right)\gamma \alpha (1-\left.u_2\right)S_2-{\lambda }_3\left(t\right){\theta }_e\left(1-{\theta }_d\right)\gamma \alpha \left(1-u_2\right)S_2\hfill \\ & -{\lambda }_2\left(t\right)\left(1-{\theta }_d\right)\left(1-{\theta }_e\right)\alpha \left(1-u_2\right)-{\lambda }_3\left(t\right){\theta }_e\left(1-{\theta }_d\right)\alpha (1-\left.u_2\right)\hfill \\ \hfill & -{\chi }_{[0,T-\tau ]}\left(t\right)\beta \left(1-u_1\right)S_2(t-\tau )e^{-\mu \tau }{\lambda }_6(t+\tau )-{\lambda }_7\left(t\right){\theta }_d\alpha (1-u_2)\hfill \\ \hfill & -{\chi }_{[0,T-\tau ]}\left(t\right)\left[1-{\theta }_d\right]\left[1-{\theta }_e\right]\gamma \alpha (1-\left.u_2\right)S_2(t-\tau )e^{-\mu \tau }{\lambda }_2(t+\tau ),\hfill \\ \hfill \frac{d_{\lambda 7}\left(t\right)}{dt}=& {\lambda }_7\left(t\right)(c+f)-{\lambda }_8\left(t\right)f,\hfill \\ \hfill \frac{d{\lambda }_8\left(t\right)}{dt}=& \left({\lambda }_8\left(t\right)-{\lambda }_4\left(t\right)\right)\alpha \left(1-u_2\right)+{\lambda }_8\left(t\right)\mu .\hfill \end{array}$$

Given the transversality conditions ${\lambda }_i\left(T\right)=0(i=1,2,3,4,5,6,7,8)$. Moreover, we provide the corresponding parameters below:

$$\begin{array}{cc}& u_1^{*}\left(t\right)=min\left\{1,max\left(0,{\overline{u}}_1\left(t\right)\right)\right\},\hfill \\ & u_2^{*}\left(t\right)=min\left\{1,max\left(0,{\overline{u}}_2\left(t\right)\right)\right\},\hfill \\ & u_3^{*}\left(t\right)=min\left\{1,max\left(0,{\overline{u}}_3\left(t\right)\right)\right\}.\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\overline{u}}_1\left(t\right)& =\frac{{\lambda }_2\left(t\right)\beta S_1(t-\tau )I_1(t-\tau )e^{-\mu \tau }+{\lambda }_6\left(t\right)\beta S_2(t-\tau )I_2(t-\tau )e^{-\mu \tau }-{\lambda }_5\left(t\right)\beta S_2{I}_2-{\lambda }_1\left(t\right)\beta S_1{I}_1}{B_1},\hfill \\ \hfill {\overline{u}}_2\left(t\right)& =\left\{{\lambda }_2\left(t\right)\left[\left(1-{\theta }_d\right)\left(1-{\theta }_e\right)\gamma \alpha S_2(t-\tau )I_2(t-\tau )e^{-\mu \tau }+\left(1-{\theta }_d\right)\left(1-{\theta }_e\right)\alpha I_2-\alpha I_1\right]\right.\hfill \\ \hfill & +{\lambda }_3\left(t\right)\left[{\theta }_e\left(1-{\theta }_d\right)\gamma \alpha S_2{I}_2+{\theta }_d\alpha I_1+{\theta }_e\left(1-{\theta }_d\right)\alpha I_2\right]\hfill \\ \hfill & +{\lambda }_6\left(t\right)\left[\left(1-{\theta }_d\right)\left(1-{\theta }_e\right)\gamma \alpha S_1(t-\tau )I_1(t-\tau )e^{-\mu \tau }+\left(1-{\theta }_d\right)\left(1-{\theta }_e\right)\alpha I_1-\alpha I_2\right]\hfill \\ \hfill & +{\lambda }_7\left(t\right)\left[{\theta }_e\left(1-{\theta }_d\right)\gamma \alpha S_1{I}_1+{\theta }_d\alpha I_2+{\theta }_e\left(1-{\theta }_d\right)\alpha I_1\right]-\left({\lambda }_8\left(t\right)-{\lambda }_4\left(t\right)\right)\left(\alpha R_1-\alpha R_2\right)\hfill \\ \hfill & \left.-{\lambda }_1\left(t\right)\left(\alpha S_1-\alpha S_2+\left(1-{\theta }_d\right)\gamma \alpha S_2{I}_2\right)-{\lambda }_5\left(t\right)\left(\alpha S_2-\alpha S_1+\left(1-{\theta }_d\right)\gamma \alpha S_1{I}_1\right)\right\}/B_2,\hfill \\ \hfill {\overline{u}}_3\left(t\right)& =\frac{\left({\lambda }_2\left(t\right)-{\lambda }_4\left(t\right)\right)d{I}_1+\left({\lambda }_6\left(t\right)-{\lambda }_8\left(t\right)\right)d{I}_2}{B_3}.\hfill \\ \hfill \end{array}$$

Proof. 

The control problem of the Hamilton function H is defined as

$$\begin{array}{cc}\hfill H\left(t\right)=& A_1(1-{\theta }_d)(1-{\theta }_e)\alpha I_1+A_2(1-{\theta }_d)(1-{\theta }_e)\alpha I_2+\frac{1}{2}\left(B_1{u}_1^2\left(t\right)+B_2{u}_2^2\left(t\right)+B_3{u}_3^2\left(t\right)\right)\hfill \\ \hfill & +{\lambda }_1\left(t\right)\left[\mathsf{\Lambda }-\beta \left(1-u_1\right)S_1{I}_1-\alpha \left(1-u_2\right)S_1+\alpha \left(1-u_2\right)S_2-\left(1-{\theta }_d\right)\gamma \alpha \left(1-u_2\right)S_2{I}_2-\mu S_1\right]\hfill \\ \hfill & +{\lambda }_2\left(t\right)\left[\beta \left(1-u_1\right)S_1(t-\tau )I_1(t-\tau )e^{-\mu \tau }+\left(1-{\theta }_d\right)\left(1-{\theta }_e\right)\gamma \alpha \left(1-u_2\right)S_2(t-\tau )I_2(t-\tau )e^{-\mu \tau }\right.\hfill \\ \hfill & +\left.\left(1-{\theta }_d\right)\left(1-{\theta }_e\right)\alpha \left(1-u_2\right)I_2-\left(b+d\left(1+u_3\right)+\alpha \left(1-u_2\right)\right)I_1\right]\hfill \\ \hfill & +{\lambda }_3\left(t\right)\left[{\theta }_e\left(1-{\theta }_d\right)\gamma \alpha \left(1-u_2\right)S_2{I}_2+{\theta }_d\alpha \left(1-u_2\right)I_1+{\theta }_e\left(1-{\theta }_d\right)\alpha \left(1-u_2\right)I_2-(c+f)Q_1\right]\hfill \\ \hfill & +{\lambda }_4\left[f{Q}_1+d\left(1+u_3\right)I_1-\alpha \left(1-u_2\right)R_1+\alpha \left(1-u_2\right)R_2-\mu R_1\right]\hfill \\ \hfill & +{\lambda }_5\left(t\right)\left[\mathsf{\Lambda }-\beta \left(1-u_1\right)S_2{I}_2-\alpha \left(1-u_2\right)S_2+\alpha \left(1-u_2\right)S_1-\left(1-{\theta }_d\right)\gamma \alpha \left(1-u_2\right)S_1{I}_1-\mu S_2\right]\hfill \\ \hfill & +{\lambda }_6\left(t\right)\left[\beta \left(1-u_1\right)S_2(t-\tau )I_2(t-\tau )e^{-\mu \tau }+\left(1-{\theta }_d\right)\left(1-{\theta }_e\right)\gamma \alpha \left(1-u_2\right)S_1(t-\tau )I_1(t-\tau )e^{-\mu \tau }\right.\hfill \\ \hfill & +\left.\left(1-{\theta }_d\right)\left(1-{\theta }_e\right)\alpha \left(1-u_2\right)I_1-\left(b+d\left(1+u_3\right)+\alpha \left(1-u_2\right)\right)I_2\right]\hfill \\ \hfill & +{\lambda }_7\left(t\right)\left[{\theta }_e\left(1-{\theta }_d\right)\gamma \alpha \left(1-u_2\right)S_1{I}_1+{\theta }_d\alpha \left(1-u_2\right)I_2+{\theta }_e\left(1-{\theta }_d\right)\alpha \left(1-u_2\right)I_1-(c+f)Q_2\right]\hfill \\ \hfill & +{\lambda }_8\left[f{Q}_2+d\left(1+u_3\right)I_2-\alpha \left(1-u_2\right)R_2+\alpha \left(1-u_2\right)R_1-\mu R_2\right],\hfill \end{array}$$

where $\lambda \left(t\right)={({\lambda }_1\left(t\right),{\lambda }_2\left(t\right),{\lambda }_3\left(t\right),{\lambda }_4\left(t\right),{\lambda }_5\left(t\right),{\lambda }_6\left(t\right),{\lambda }_7\left(t\right),{\lambda }_8\left(t\right))}^T$ is the adjoint variable satisfying

$\frac{d{\lambda }_1\left(t\right)}{dt}=-\frac{\partial H}{\partial S_1}\left(t\right)-{\chi }_{[0,T-\tau ]}\left(t\right)\frac{\partial H}{\partial S_1(t-\tau )}(t+\tau ),$

$\frac{d{\lambda }_2\left(t\right)}{dt}=-\frac{\partial H}{\partial I_1}\left(t\right)-{\chi }_{[0,T-\tau ]}\left(t\right)\frac{\partial H}{\partial I_1(t-\tau )}(t+\tau ),$

$\frac{d{\lambda }_3\left(t\right)}{dt}=-\frac{\partial H}{\partial Q_1}\left(t\right),$

$\frac{d{\lambda }_4\left(t\right)}{dt}=-\frac{\partial H}{\partial R_1}\left(t\right),$

$\frac{d{\lambda }_5\left(t\right)}{dt}=-\frac{\partial H}{\partial S_2}\left(t\right)-{\chi }_{[0,T-\tau ]}\left(t\right)\frac{\partial H}{\partial S_2(t-\tau )}(t+\tau ),$

$\frac{d{\lambda }_6\left(t\right)}{dt}=-\frac{\partial H}{\partial I_2}\left(t\right)-{\chi }_{[0,T-\tau ]}\left(t\right)\frac{\partial H}{\partial I_2(t-\tau )}(t+\tau ),$

$\frac{d{\lambda }_7\left(t\right)}{dt}=-\frac{\partial H}{\partial Q_2}\left(t\right),$

$\frac{d{\lambda }_8\left(t\right)}{dt}=-\frac{\partial H}{\partial R_2}\left(t\right).$

On $[0,T-\tau ]$, ${\chi }_{[T-\tau ]}\left(t\right)$ are indicator functions under the transversality conditions ${\lambda }_i\left(T\right)=0(i=1,2,3,4,5,6,7,8)$, which fullfill

$${\chi }_{[T-\tau ]}\left(t\right)=\left\{\begin{array}{cc}1,\hfill & t\in [0,T-\tau ],\hfill \\ 0,\hfill & otherwise.\hfill \end{array}\right.$$

The optimality equation

$$\frac{\partial H}{\partial u_1}\left(t\right)=\frac{\partial H}{\partial u_2}\left(t\right)=\frac{\partial H}{\partial u_3}\left(t\right)=0.$$

From the optimality equation,

$$\begin{array}{cc}\hfill B_1{u}_1& +{\lambda }_1\left(t\right)\beta S_1{I}_1-{\lambda }_2\left(t\right)\beta S_1(t-\tau )I_1(t-\tau )e^{-\mu \tau }+{\lambda }_5\left(t\right)\beta S_2{I}_2-{\lambda }_6\left(t\right)\beta S_2(t-\tau )I_2(t-\tau )e^{-\mu \tau }=0,\hfill \\ \hfill B_2{u}_2& +{\lambda }_1\left(t\right)\left(\alpha S_1-\right.\left.\alpha S_2+\left(1-{\theta }_d\right)\gamma \alpha S_2{I}_2\right)\hfill \\ \hfill & -{\lambda }_2\left(t\right)\left[\left(1-{\theta }_d\right)\left(1-{\theta }_e\right)\gamma \alpha S_2(t-\tau )I_2(t-\tau )e^{-\mu \tau }+\left(1-{\theta }_d\right)\left(1-{\theta }_e\right)\alpha I_2-\alpha I_1\right]\hfill \\ \hfill & -{\lambda }_3\left(t\right)\left[{\theta }_e\left(1-{\theta }_d\right)\gamma \alpha S_2{I}_2+{\theta }_d\alpha I_1+{\theta }_e\left(1-{\theta }_d\right)\alpha I_2\right]+{\lambda }_5\left(t\right)\left(\alpha S_2-\alpha S_1+\left(1-{\theta }_d\right)\gamma \alpha S_1{I}_1\right)\hfill \\ \hfill & -{\lambda }_6\left(t\right)\left[\left(1-{\theta }_d\right)\left(1-{\theta }_e\right)\gamma \alpha S_1(t-\tau )I_1(t-\tau )e^{-\mu \tau }+\left(1-{\theta }_d\right)\left(1-{\theta }_e\right)\alpha I_1-\alpha I_2\right]\hfill \\ \hfill & -{\lambda }_7\left(t\right)\left[{\theta }_e\left(1-{\theta }_d\right)\gamma \alpha S_1{I}_1+{\theta }_d\alpha I_2+{\theta }_e\left(1-{\theta }_d\right)\alpha I_1\right]-\left({\lambda }_8\left(t\right)-{\lambda }_4\left(t\right)\right)\left(\alpha R_1-\alpha R_2\right)=0,\hfill \\ \hfill B_3{u}_3& -\left({\lambda }_2\left(t\right)-{\lambda }_4\left(t\right)\right)d{I}_1-\left({\lambda }_6\left(t\right)-{\lambda }_8\left(t\right)\right)d{I}_2=0.\hfill \end{array}$$

Thus we have,

$$\overline{u}\left(t\right)=\left({\overline{u}}_1\left(t\right),{\overline{u}}_2\left(t\right),{\overline{u}}_3\left(t\right)\right).$$

Given that $u^{*}\left(t\right)$ is part of U, we see

$$\begin{array}{cc}\hfill & u_1^{*}\left(t\right)=min\left\{1,max\left(0,{\overline{u}}_1\left(t\right)\right)\right\},\hfill \\ \hfill & u_2^{*}\left(t\right)=min\left\{1,max\left(0,{\overline{u}}_2\left(t\right)\right)\right\},\hfill \\ \hfill & u_3^{*}\left(t\right)=min\left\{1,max\left(0,{\overline{u}}_3\left(t\right)\right)\right\}.\hfill \end{array}$$

(8)

□

We present some numerical outcomes that belong to the optimality system. In order to determine the optimal control and the state system, we use the scheme proposed in [47,50], which is based on the forward-backward sweep scheme. Furthermore, many other algorithms describe approximation methods for obtaining optimal control [51,52]. We choose the figure from Table 2 as the parameter values in the optimal control. The weights in the goal function $A_1=A_2=10$, $B_1=20$, $B_2=80$, $B_3=120$, we will explain these values in next section. The results in Figure 5 suggest that optimal control effectively controls the development and deterioration of the disease. In fact, the intensity of various measures implemented by the government is different, such as wearing masks. Combined with the prevention and control situation in the United States at that time[53,54,55], we assume that the upper limit of these three measures is $0.11$. Figure 6 shows the control intensity of each measure in the control system over time, and we can see the magnitude of these three measures when they are taken at the same time. $u_2\left(t\right)$ (controlling population migration) can end about 20 days earlier than $u_1\left(t\right)$ (wearing masks) and $u_3\left(t\right)$ (take medication and actively pursue treatment). This means that we can reduce the control intensity of $u_2\left(t\right)$ (controlling population migration) sometime in advance when all three measures are implemented. Implementing various strategies has a certain cost, and we discuss the most cost-effective strategies in the next section.

5. Cost-Effectiveness Analysis

We use three ways, which are: infection averted ratio (IAR) [56], average cost-effectiveness ratio (ACER) [57], and incremental cost-effectiveness ratio (ICER) [57] to complete the cost-effectiveness analysis. Below are the interpretations of the three ways:

(1) IAR can be denoted as

$$\mathrm{IAR}=\frac{\mathrm{Number}\mathrm{of}\mathrm{people}\mathrm{avoiding}\mathrm{infection}}{\mathrm{Number}\mathrm{of}\mathrm{people}\mathrm{who}\mathrm{have}\mathrm{recovered}\mathrm{after}\mathrm{infected}}.$$

We regard the control way whose IAR is the highest as the economical choice [57].

(2) ACER is told as

$$\mathrm{ACER}=\frac{\mathrm{The}\mathrm{overall}\cos \mathrm{t}\mathrm{resulting}\mathrm{from}\mathrm{the}\mathrm{carrying}\mathrm{out}\mathrm{of}\mathrm{a}\mathrm{specific}\mathrm{intervention}\mathrm{strategy}}{\mathrm{Numeber}\mathrm{of}\mathrm{people}\mathrm{who}\mathrm{have}\mathrm{avioding}\mathrm{infection}\mathrm{after}\mathrm{performing}\mathrm{strategy}}.$$

We assess the money which is used to perform a particular intervention strategy is

$$c\left(u\right)=\frac{1}{2}{\int }_0^T\sum _{i=1}^3{B}_i{u}_i^{2}dt.$$

(3) ICER assess the variations in cost and health advantages of two distinct intervention strategies under the same limited resourses. Viewing strategies m and n as two viral control intervention strategies, then ICER is expressed as

$$\mathrm{ICER}=\frac{\mathrm{Total}\cos \mathrm{t}\mathrm{variations}\mathrm{in}\mathrm{strategies}m\mathrm{and}n}{\mathrm{Control}\mathrm{benefits}\mathrm{variations}\mathrm{in}\mathrm{strategies}m\mathrm{and}n}.$$

The infection that is prevented with the controls during the time period T is expressed as

$$Infectionaverted=T\left(I_1\left(0\right)+I_2\left(0\right)\right)-{\int }_0^T{I}_1{ }^{*}\left(t\right)+I_2{ }^{*}\left(t\right)dt.$$

Here, $B_i,(i=1,2,3)$ is related to the per capita unit cost of the measures $u_i^{*},i=1,2,3$. Let us take the example of Washington. Washington mandates wearing masks in public settings [53]. However, due to cost and availability, individuals in low-income environments are likelier to use reusable cloth masks instead of disposable medical masks. There is still uncertainty in mask-wearing practices, such as the duration of wearing a mask [58]. Therefore, we can reasonably assume that the daily cost of a mask for each person is approximately $0.06$ dollars, which amounts to approximately $21.9$ dollars per year.

Due to early stay-at-home orders in the United States, population mobility and migration decreased [54], and local governments also provided subsidies to the public during the pandemic [55]. We assume the average annual subsidy per person is approximately 80 dollars. As the prices of different drugs for treating and preventing COVID-19 vary [59,60], we assume that the annual cost for each person to purchase drugs for COVID-19 treatment is approximately 120 dollars.

Because of the specific numerical values given by the fourth column of Table 3, it is clear that Strategy 1 (wearing a mask) is the most cost-effective measure.

According to the ACER values in column 5 of Table 3, Strategy 1 is the lowest, so the most cost-effective Strategy is Strategy 1 (wearing masks).

Next, we calculate the ICER value, where we rank the number of people we want to avoid infection in increasing order.

The ICER is computed as

$$\begin{array}{cc}\hfill & ICER\left(3\right)=\frac{2.188\times {10}^7-0}{1.823\times {10}^5-0}=120.021,\hfill \\ \hfill & ICER\left(2\right)=\frac{4.921\times {10}^7-2.188\times {10}^7}{6.152\times {10}^5-1.823\times {10}^5}=63.132,\hfill \\ \hfill & ICER\left(6\right)=\frac{7.109\times {10}^7-4.921\times {10}^7}{6.157\times {10}^5-6.152\times {10}^5}=599.452,\hfill \\ \hfill & ICER\left(1\right)=\frac{1.515\times {10}^7-7.109\times {10}^7}{6.921\times {10}^5-6.517\times {10}^5}=-1385.650,\hfill \\ \hfill & ICER\left(5\right)=\frac{3.703\times {10}^7-1.515\times {10}^7}{7.080\times {10}^5-6.921\times {10}^5}=1376.100,\hfill \\ \hfill & ICER\left(4\right)=\frac{6.436\times {10}^7-3.703\times {10}^7}{8.301\times {10}^5-7.080\times {10}^5}=223.833,\hfill \\ \hfill & ICER\left(7\right)=\frac{8.624\times {10}^7-6.436\times {10}^7}{8.787\times {10}^5-8.301\times {10}^5}=450.210.\hfill \end{array}$$

The results are presented in Table 3. Compared with the ICER value of Strategy 6, Strategy 5 is significantly superior. Therefore, Strategy 5 is removed from the control list. Similarly, in sequence, we can eliminate strategies 6, 7, 4, 3, and 2. Detailed calculation process can be found in .

Therefore, to sum up the analysis, we can know that when one considers the following control measures: $u_1$ (wearing masks ), $u_2$ (policy interventions and city lockdowns), $u_3$ (taking medication and actively seeking treatment). $u_1$ (wearing masks) is the least incremental cost-effective, providing the optimal cost over all other strategies on a large scale.

6. Discussion

Considering that different genotypes of viruses have different incubative periods, we select different viruses circulating in different periods of COVID-19. The original SARS-CoV-2, Delta, and Omicron incubation periods are 5.3 days, 4 days, and 2.5 days [42,61]. We only discuss the influence of the length of the incubative period on the number of patients, so other parameters in the model are determined, and Figure 7 is obtained. To observe the extent of infection for different incubative periods, we have extracted data on the total quantity of infected individuals in 300 days with different incubative periods in two regions from Figure 7, as shown in Figure 8. We obtain conclusions similar to those of [62]. As the disease incubative period becomes shorter, the disease tends to end earlier, leading to a larger scale of infection.

To verify the effectiveness of transportation-related infection and entry-exit screenings, we take the below four reproduction numbers into consideration,

$${\mathcal{R}}_0^0=\frac{\beta S^0{e}^{-\mu \tau }}{b+d},{\mathcal{R}}_0^{\gamma }=\frac{(\beta +\gamma \alpha )S^0{e}^{-\mu \tau }}{b+d},{\mathcal{R}}_0^{\gamma {\theta }_d}=\frac{(\beta +(1-{\theta }_d)\gamma \alpha )S^0{e}^{-\mu \tau }}{b+d+{\theta }_d\alpha },$$

$${\mathcal{R}}_0^{\gamma {\theta }_e}=\frac{(\beta +(1-{\theta }_e)\gamma \alpha )S^0{e}^{-\mu \tau }}{b+d+{\theta }_e\alpha },{\mathcal{R}}_0=\frac{(\beta +\left(1-{\theta }_e\right)\left(1-{\theta }_d\right)\gamma \alpha )S^0{e}^{-\mu \tau }}{b+d+{\theta }_d\alpha +{\theta }_e\left(1-{\theta }_d\right)\alpha }.$$

The ${\mathcal{R}}_0^0$ is the basic reproduction number for each region when no individuals are moving between regions, and ${\mathcal{R}}_0^{\gamma }$ is the basic reproduction number when there is transportation-related infection. ${\mathcal{R}}_0^{\gamma {\theta }_e}$ is the value of the basic reproduction number for transportation-related infection and entry screening, while ${\mathcal{R}}_0^{\gamma {\theta }_d}$ is the one for transportation-related infection and exit screening, the ${\mathcal{R}}_0$ is the basic reproduction number for model (1) with transportation-related infection and entry-exit screenings. In order to consider the effects of ${\theta }_d$ and ${\theta }_e$ on the disease, we take the basic reproduction number ${\mathcal{R}}_0^0$ with no population movement between the two cities as the base. We fix all the parameters except ${\theta }_d$ and ${\theta }_e$ and observe the change in the basic reproduction number. We obtain Figure 9 of ${\mathcal{R}}^{*}$ with respect to ${\theta }_d$ and ${\theta }_e$. The results in Figure 9 indicate that transportation-related infections increase the basic reproduction number and enhance disease transmission. Entry or exit screening can inhibit disease transmission by reducing the basic reproduction number.

7. Conclusion

We have established a time-delayed SIQR epidemiological model for infectious diseases with transportation-related infections between two regions. Although many diseases, such as SARS and influenza, can be studied using this model, we focus on using COVID-19 as a specific example to consider our model. Despite the end of the COVID-19 pandemic, it is still essential to remain vigilant and understand how to prevent such emerging diseases in real-time.

We obtain some useful conclusions. Model (1) has a unique global asymptotically stable disease-free equilibrium when ${\mathcal{R}}_0<1$ and a unique endemic equilibrium exists when ${\mathcal{R}}_0>1$. We use COVID-19 data from Washington, USA, between June 10, 2020, and September 10, 2020, to fit the model and obtain a set of optimal parameters. We conduct a sensitivity analysis on the basic reproduction number (${\mathcal{R}}_0$) and obtain that $\beta$, $\alpha$, and d are more sensitive parameters compared to other parameters. Based on these sensitive parameters, we propose three control measures: (1) Wearing masks; (2) Policy interventions and city lockdowns; (3) Taking medication and actively pursuing treatment. Considering infections between regions through travel, we formulate an optimal control problem and find an optimal solution to reduce the number of undetected patients traveling between the two regions. We use Pontryagin’s minimum principle to obtain the optimal solution to achieve our goal. Subsequently, we obtain some numerical results using an optimal control approach based on the forward-backward sweep scheme. The results indicate that the control strategies help reduce infections. We can reduce the control intensity of $u_2\left(t\right)$ (controlling population migration) sometime in advance when all three measures ($u_1\left(t\right)$, $u_2\left(t\right)$, $u_3\left(t\right)$) are implemented. Additionally, due to a decrease in the number of infected individuals, there is a relative reduction in isolating patients in hospitals, which aids in reducing the required medical resources for significant emergent diseases. We compare the three control strategies using three methods: Infection avoidance ratio (IAR), Average cost-effectiveness ratio (ACER), and Incremental cost-effectiveness Ratio (ICER). Strategy 1 (wearing masks) is the most cost-effective measure. Therefore, in the face of some major emerging infectious diseases, such as SARS, COVID-19, and other diseases. We can all encourage and require people to wear masks.

Finally, we discuss the effects of varying the disease’s incubative period and entry-exit screening on the disease. We observe that a shorter incubative period result in a shorter disease duration but leads to a larger scale of infection, and the peaks are reduced. This has a certain guiding role in the face of many major emerging infectious diseases between two regions. We can formulate some effective prevention and control measures according to the length of the disease’s incubative period, which can reduce the large-scale spread of the disease to a certain extent. Transportation-related infections increase the basic reproduction number (${\mathcal{R}}_0$) and enhance disease spread. Entry-exit screening have the ability to curb disease spread by lowering the basic reproduction number. Indeed, our model also has many factors that have not been considered. For instance, when the total populations of the two regions in the model are different, and the migration rates between these regions also vary, what is the impact on the spread of this disease? There are many uncertain factors in disease transmission, so what influence do they have on the disease? The risk of reinfection in recovered patients, differences in disease detection coverage in different regions, and transmission in n regions are all questions that deserve to be explored in the future.