Mathematical analysis and optimal control of an epidemic model with vaccination and different infectivity

1.   Introduction

Infectious diseases have always been an important part of human history, where various societies and countries have all witnessed devastating epidemics that resulted in millions of casualties. A notable example is the Spanish flu, which originated in Spain in 1918 and claimed the lives of around 40 million people. The COVID-19 pandemic, which began at the close of 2019, was declared a major global public health concern by the World Health Organization (WHO). Due to the virus's rapid global transmission, it has reached over 200 countries. As of July 22, 2023, there have been approximately 760 million confirmed cases of the virus worldwide, with approximately 6.9 million reported deaths. This has posed a significant threat to the well-being of humans, giving rise to major public health challenges and economic crises ^[1,2]. Furthermore, it has placed additional strain on global healthcare systems ^[3].

Vaccinations stand as the most effective approach for curtailing the transmission of diseases within a population, thus representing a long-term strategy in disease prevention. The introduction of vaccines has played a pivotal role in effectively managing several infectious diseases, including smallpox, cholera, and plague. Notably, the oral polio vaccine, often referred to as the sugar pill, has been instrumental in halting the spread of indigenous wild poliovirus, thereby safeguarding thousands of children from physical disabilities. The field of mathematical epidemiology has witnessed significant and remarkable advancements in modeling vaccination strategies. Lhimn et al. ^[4] developed a two-group Susceptible-Infected-Recovered (SIR) mathematical model to assess the effectiveness of split-dose vaccines in disease control. Their findings indicated that when vaccine efficacy reaches a certain threshold, vaccines exhibit a positive impact on both reducing the size of epidemics and delay the peak outbreak timing. Gao et al. ^[5] employed a two-sex mathematical model to assess the impact of vaccinations on Human Papilloma Virus (HPV) infection. Their research demonstrated that implementing a variable vaccination strategy could substantially contribute to reducing HPV prevalence.

Along with the successful development of the COVID-19 vaccine, many authors have assessed its impact on the transmission of novel coronaviruses. Thompson et al. ^[6] proposed a detailed mathematical model to study the impact of various preventive and control measures and vaccinations on the spread of COVID-19. Harizi et al. ^[7] used data from Canada to discuss the impact of universal vaccinations on the transmission of COVID-19. Diagne et al. ^[8] developed a mathematical model to discuss the effect of the vaccination rate and vaccine inefficiency on the spread of epidemic. Paul et al. ^[9] used a six-compartment Susceptible-Vaccinated-Exposed-Infected-Hospitalized-Recovered (SVEIHR) model to analyze the effect of vaccines on the spread of COVID-19, and showed that vaccination can help to reduce social stress, which, in turn helps to strengthen the immune system. Rocha et al. ^[10] used an expanded version of the Susceptible-Infected-Recovered-Death (SIRD) model, including the vaccination effect, to study the impact of vaccines on COVID-19 transmission and concluded that vaccination results in a 5% relative reduction in the total number of deaths. Makhoul et al. ^[11] evaluated the effect of vaccines against infection and disease and concluded that an increased vaccination rate plays an important role in controlling the spread of the disease. Furthermore, several studies have used mathematical models to explore optimized vaccine distribution ^[12,13]. Brody et al. showed that the best strategy is to prioritize for older adults and high-exposure populations ^[14]. Li et al. ^[15] illustrated that universal booster injections are a cost-effective strategy in fully vaccinated persons aged $\geq$ 65 years.

Optimal control theory encompasses several crucial principles for guiding the management of epidemic spread through practical control measures. Numerous researchers have applied this theory to explore the optimal strategy for epidemic models. Alrabaiah et al. ^[16] developed an optimized control model with three control variables and found that an integrated control strategy, including isolation, treatment and vaccination, is the most effective in eradicating Hepatitis B virus (HBV). Odionyenma et al. ^[17] employed optimal control theory to conclude that a control strategy combining treatment and prevention is the most cost-effective approach in controlling Chlamydia trachomatis, as per a cost-effectiveness analysis.

However, upon reviewing the aforementioned literature, it becomes apparent that most works tend to focus on either vaccine failure rates or vaccine-induced immune decline while neglecting the possibility of asymptomatic infections progressing to symptomatic ones. Furthermore, there is limited research on optimal comprehensive control strategies that consider multiple factors. In light of these considerations, our study concentrates on evaluating the effects of reducing social distancing, increasing vaccination rates, and enhancing hospitalization rates on an epidemic spread within the population. We establish a transmission model that accounts for vaccination, distinct infection states, and varying infectivity. Through an exploration of its dynamic behavior and numerical simulations, we propose prevention and control strategies for infectious diseases with similar characteristics.

The structure of this paper is as follows. In Section 2, we mainly make some basic assumptions and introduce the model. In Section 3, we discuss the stability of the disease-free equilibrium, and then conduct numerical simulations to discuss the effects of vaccine-related parameters on the control reproduction number and the spread of an epidemic. In Section 4, we perform sensitivity analyses of a control reproduction number, vaccinated individuals ($V$), and symptomatic infections ($I$) to derive three parameters that influence most on-disease transmission. In Section 5, the corresponding optimal control problem is considered, and several different control combinations are compared in terms of the control effectiveness and the cost effectiveness. This paper ends with the conclusion and discussion.

2.   Model formulation

Assuming that the vaccine, which is available and is partially effective on preventing infection, we divide the total population ($N$) into susceptible ($S$), vaccinated ($V$), exposed ($E$), asymptomatic ($A$), symptomatic ($I$), hospitalized ($H$), and recovered ($R$). Based on the divided compartments, we will formulate a transmission epidemic model.

It is assumed that susceptible individuals are recruited into the population at rate $\Lambda$, and receive vaccinations at rate $\varphi$. Symptomatic and asymptomatic infections can transmit a disease at rates $\beta\varepsilon_I$ and $\beta\varepsilon_A$, respectively, where $\beta$ is the basic transmission rate, and $\varepsilon_I$ and $\varepsilon_A$ are adjusted factors for the transmission ability of symptomatic and asymptomatic infections, respectively. Although the hospitalized cases are infectious, we assume that they will not reinfect others due to the increased isolation conditions. We assume that the vaccine may not offer complete protection and vaccinated individuals can become infected with a rate $\sigma\beta$, where $\sigma \in [0, 1]$ is the modified factor and denote $\overline{\sigma} = 1-\sigma$. The vaccinated individuals can return to be susceptible after the vaccine-induced immunity wanes at a rate $\theta$. After contact with infectious individuals, they will enter into an exposed compartment and experience an exposed period ($1/\alpha$). Due to individual heterogeneity, exposed individuals with a fraction $p$ develop to be asymptomatically infectious while others with a fraction $1-p$ become symptomatically infectious. Furthermore, asymptomatic infections progress to symptomatic infections at a rate $\delta$. The asymptomatic and symptomatic infections can naturally recover at rates $r_A$ and $r_I$, respectively. We assume that asymptomatic infections do not need hospitalization and cannot cause an additional death. Alternatively, symptomatic infections require hospitalization at a rate $k$, and has an additional death at a rate $d_I$. The hospitalized individuals can recover and may die at rates $r_H$ and $d_H$, respectively. All recovered individuals enter the recovered compartment and have no secondary infections. Moreover, individuals in all compartments have a natural death rate $d$. The detailed schematic diagram is shown in Figure 1. Table 1 summarizes the model parameters. Based on the aforementioned assumptions and Figure 1, we formulate the following transmission epidemic model in the form of a nonlinear differential equations given by the following:

with initial conditions $\bar{X_0} = (S(0), V(0), E(0), A(0), I(0), H(0), R(0))$, where $S(0)$ is positive and the other initial values are non-negative. For simplicity, we introduce the notations in model (2.1) as follows:

It should be noted that $H$ and $R$ are independent of the other equations in model (2.1). Therefore, we will focus on the following subsystem:

Denote the initial value of system (2.2) as $X(0) = (S(0), V(0), E(0), A(0), I(0))$.

3.   Model analysis

In this section, we will focus on the mathematical analysis of model (2.2) and the associated simulations to explore the influences of vaccine-related parameters on the spread of an epidemic.

3.1.   Positive invariance

Clearly, system (2.2) with the non-negative initial conditions is wellposed. It is easy to check that the solution of system (2.2) is always unique and non-negative in the region $\Omega = \{(S, V, E, A, I)\in\mathbb{R}^5_+: 0 < N \leq \Lambda/d \}$ for all time $t\geq 0$.

Lemma 3.1. The region $\Omega$ is positively invariant of model (2.2).

It follows from model (2.2) that we can add all the equations and observe that the total population satisfies $dN/dt = \Lambda-dN-(d_II+d_HH)\leq \Lambda-dN$, which implies that $N \leq \Lambda / d$.

3.2.   Two reproduction numbers

Clearly, model (2.2) has a unique disease-free equilibrium $E_0 = (S_0, V_0, 0, 0, 0)$, where $S_0 = \Lambda \phi_2/dD_1$, $V_0 = \Lambda \varphi/dD_1$ and $D_1 = \phi_2+\varphi$.

By the next generation matrix approach ^[30,31], we obtain the following two matrices:

Furthermore, linearizing system (2.2) at $E_0$ produces the following two sensitivity matrixes of $\mathcal{F}$ and $\mathcal{V}$:

where $D_2 = \phi_2+\sigma \varphi$. The $\Re_V$ is obtained by the spectral radius of the next generation matrix $(FV^{-1})$ with the following expression:

It can be seen that $D_1$ and $D_2$ are vaccine-related parameters. If we do not take the vaccination strategy into account, that is, $\varphi = \sigma = \theta = 0$, we have $D_1 = D_2 = d$. In this case, we can obtain the corresponding basic reproduction number, denoted by $\Re_0$, as follows:

Clearly, $\Re_0$ consists of three parts, which implies that there are three transmission routes. Specifically, the first item denotes the number of new, asymptomatic infections during his/her duration period. $\beta \varepsilon_A$ is the transmission rate of asymptomatic infections. The probability of survival and transfer from exposed individuals to asymptomatic infections is $p \alpha /\phi_3$, whereas $1/\phi_4$ refers to the average duration of transmission in asymptomatic infections. The second and third items denote the number of new, symptomatic infections during his/her average duration period. There are two ways. $\beta \varepsilon_I$ is the infection rate of symptomatic infections. First, survival and transmission from exposed individuals to asymptomatic infections is $p \alpha / \phi_3$, and then transmission from asymptomatic infections to symptomatic infections is $\delta/\phi_4$. Alternatively, survival and transmission from exposed individuals to symptomatic infections is $(1-p)\alpha /\phi_3$. $1/\phi_5$ is the mean infection duration of symptomatic infections.

3.3.   The stability of $E_0$

Theorem 3.1. The disease-free equilibrium $E_0$ is globally asymptotically stable if $\Re_{V} < 1$ and $\sigma = 1$.

Proof. Step1. Local stability.

The Jacobian of model (2.2) at $E_0$ is as follows:

The associated characteristic equation of $E_0$ is given by $C(\lambda) = (\lambda +d)(\lambda +D_2)f(\lambda) = 0$, where $f(\lambda) = \lambda^3+\Psi_1\lambda^2+\Psi_2 \lambda+\Psi_3$ with

Clearly, $\lambda_1 = -d$ and $\lambda_2 = -D_2$ are two negative roots of $C(\lambda) = 0$, and the other three roots are determined by $f(\lambda) = 0$. Now, we focus on $f(\lambda) = 0$. From $\Re_V < 1$, one can obtain the following:

Thus, $\Psi_2 > 0$. Since the only negative item $\phi_3\phi_4\phi_5$ can be counterbalanced by one of items in $\Psi_1\Psi_2$, $\Psi_1\Psi_2-\Psi_3 > 0$ holds.

By Hurwitz criteria, one has that all eigenvalues of $f(\lambda) = 0$ have negative real parts. Moreover, $C(\lambda) = 0$ has five roots with negative real parts. Hence, $E_0$ is locally asymptotically stable if $\Re_V < 1$.

Step2. Globally stability.

Define the Lyapunov function as follows:

where

It is easily seen that $\mathcal{L}$ is non-negative in $\Omega$. The time differentiation of $\mathcal{L}$ along with the solution of model (2.2) gives the following:

Obviously, $\Re_V = \Re_0$ if $\sigma = 1$. Hence, one has

if $\Re_V < 1$. Moreover $d\mathcal{L}/dt = 0$ only if $A = I = 0$. By LaSalle's invariance principle, $E_0$ is globally asymptotically stable if $\Re_V < 1$ and $\sigma = 1$. □

3.4.   The influences of vaccine on the spread of epidemic

In this section, based on the parameter values in Table 1, we perform several numerical simulations to explore how the vaccine influences the control reproduction number and the infection size of the epidemic.

First, Figure 2 shows the influence of vaccines on the control reproduction number. From Figure 2, one can observe the control reproduction number is positively correlated with the immune waning rate, and negatively correlated with the effectiveness rate of the vaccine. Figure 2(a) implies that solely increasing vaccination rate cannot reduce the value of $\Re_V$ to below 1. A clear observation is that the vaccine effectiveness rate is more sensitive than the vaccination rate, as shown in Figure 2(b). Concretely, when the value of the vaccination rate is fixed, the value of the control production number decreases linearly with the increasing effectiveness of the vaccine. This suggests that promoting the development of more effective and less wanning vaccines is a more important way to cope with the epidemic.

Next, Figure 3 shows how the vaccination rate influences the finial size of vaccinated individuals and symptomatic infections. In Figure 3, we set $\varphi = 0.038$ as the baseline value. We observe the changes in behavior by increasing $\varphi$ by 30%, 60% and 90%, which are summarized in Table 2. From Figure 3, we observe that the finial size of vaccinated individuals increases as the vaccination rate is enhanced, while the size of symptomatic infections does the opposite. Specially, a clear comparison from Table 2 shows that a 30% increase in the vaccination rate leads to a 17.57% increase in the peak value of vaccinated individuals and a 15.73% decrease in the peak value of symptomatic infections; however, if the vaccination rate increases by 90%, the peak value of vaccinated individuals will increase by 47.2% and the peak value of symptomatic infections will decrease by 41.39%. Therefore, additional measures, such as extensively advertising and more encouraging policies, are taken to increase the vaccination rate, which will effectively reduce the peak value of infections and contribute to alleviating the epidemic.

Finally, Figure 4 shows the influences of the immune waning rate on the finial size of epidemic. To do so, we set $\theta = 0.0057$ as the baseline value, and observe the changes by reducing $\theta$ by 30%, 60% and 90%. Figure 4 shows that the effects on the peak value of vaccinated individuals and symptomatic infections are not so significant. The additional simulation values are summarized in Table 2, which shows that if the immune waning rate is reduced by 90%, the peak value of vaccinated individuals will increase by 4.77% and the peak value of symptomatic infections will decrease by 4.87%. Comparing Figure 3(b) and Figure 4(b), one can conclude that improving the vaccination rate is a more conductive and effective measure to control the epidemic.

4.   Sensitivity analysis

4.1.   Sensitivity analysis of $\Re_V$

A sensitivity analysis shows the importance of parameters in relation to disease transmission. This information is essential not only for experimental design, but also for the reduction of complex non-linear models. Two types of sensitivity analyses are used, namely local and global sensitivity analyses. The local sensitivity analysis technique examines the local response of the output by changing one input parameter at a time and keeping the other parameters at their central values ^[32]. The global sensitivity analysis perturbs the input parameters of the model on a large scale to quantify the overall impact of the model inputs on the model output ^[33,34]. In this section, we employ a sensitivity analysis to explore how parameters influence the control reproduction number. To do so, we only select eight key parameters, $\beta$, $\sigma$, $\theta$, $\alpha$, $r_I$, $r_A$, $\varphi$ and $k$, which can be interpreted by some actual measurements.

4.1.1.   Local sensitivity analysis

The local sensitivity analysis can be quantified by the sensitivity index, which is defined as the forward sensitivity index ^[33] calculated by $\mathcal{\gamma}_\xi^{\mathcal{R}_V} = \frac{\partial \mathcal{R}_V}{\partial \xi}\frac{\xi}{\mathcal{R}_V}$.

The analytical expressions of each parameter sensitivity index for $\Re_{V}$ are as follows:

with $D_3 = p\varepsilon_A\phi_5+(1-p) \varepsilon_I\phi_4+p\varepsilon_I\delta$.

It should be noted that some sensitivity indicators depend on one or more parameters, while partial sensitivity indicators are constant and do not depend on any parameter value. Based on the parameter values in Table 1, the calculated indicators are shown in Table 3 and the absolute sensitivity indicators are shown in Figure 5(a).

From Table 3, the parameters with positive sensitivity indicators are $\beta$, $\sigma$, $\theta$, and $\alpha$, which means that they have a positive effect on $\Re_V$, that is, the value of the control reproduction number will increase/decrease as parameters $\beta$, $\sigma$, $\theta$, and $\alpha$ increase/decrease. Furthermore, the sensitivity indicators show that the transmission rate $\beta$ is the most sensitive positive parameter. The parameters with negative sensitivity indicators are $r_I$, $r_A$, $\varphi$ and $k$, which means that they have a negative effect on $\Re_V$, that is, the increase of these parameters will decrease the value of $\Re_V$, and vice versa. Sensitivity indicators show that the hospitalization rate $k$ is the most sensitive negative parameter.

4.1.2.   Global sensitivity analysis

Compared to the local sensitivity analysis, the global sensitivity analysis extends the range of model input parameters. A global sensitivity analysis is performed by calculating the partial rank correction coefficient (PRCC) ^[35,36] to sample by the Latin hypercube sampling (LHS) method.

The PRCC of the control reproduction number is described in Figure 5(b). One can observe that the parameters $\beta$, $\sigma$, $\theta$, and $\alpha$ are positively correlated with $\Re_V$, and the parameters $k$, $\varphi$, $r_A$ and $r_I$ are negatively correlated with $\Re_V$. We can further find that the transmission rate $\beta$ is the most sensitive positive indicator, and the hospitalization rate $k$ and the vaccination rate $\varphi$ are the most sensitive negative indicators. Therefore, in order to effectively reduce the value of the control reproduction number, it is necessary to carry out control measures to deduce the transmission rate, and increase the hospitalization and vaccination rates, which is consistent with the results of the local sensitivity analysis.

4.2.   Sensitivity analysis of $V$ and $I$

To assess the importance of parameters on the entire duration of the epidemic, we focus on the state variables $V$ and $I$. Based on the parameter values in Table 1, we chose a time range from 0 to 100 for the PRCCs, which is plotted in Figure 6, where the gray region indicates that there is no significant difference from zero.

From Figure 6(a), one can observe that the parameters fall into three categories. The parameters $\varphi$ and $\beta$ belong to the first category, where $\varphi$ is always positively correlated with $V$ and $\beta$ is always negatively correlated with $V$. The parameter $\theta$ belongs to the second category, where $\theta$ is negatively correlated with $V$ at the beginning, then enters the grey region and leaves it again as time changes, and finally remains negatively correlated with $V$ and eventually stabilizes. Any additional parameters fall into the third category, where they start in the gray region, either increase or decrease over time, then leave the gray region and eventually stabilize.

From Figure 6(b), the parameters are also divided into three categories. The first category includes $\beta$, $\alpha$, $r_I$ and $k$, where $\beta$ is always positively correlated with $I$ and parameters $r_I$, $\alpha$ and $k$ are always negatively correlated with $I$. The second category contains $\sigma$, $\theta$ and $r_A$, which start in the gray region, either increase or decrease over time, and then leave the gray region. The third category is the parameter $\varphi$, which is initially negatively correlated with $I$ and moves into the grey region over time and is uncorrelated with $I$.

5.   Model with optimal control

From the above analysis, we can observe that the three most important parameters affecting the transmission threshold are the transmission rate $\beta$, the vaccination rate $\varphi$ and the hospitalization rate $k$. Thus, an optimal control problem to mitigate the COVID-19 epidemic is proposed by introducing three time-varying control measures denoted by $u_i(t)$ $(i = 1, 2, 3)$, where $u_1(t)$ denotes reducing social distance to decrease the transmission ability, $u_2(t)$ denotes enhancing the vaccination rate, $u_3(t)$ denotes enhancing the hospitalization rate. Thus, the transmission rate is reduced to $(1 - u_1 (t))\beta$, the vaccination rate is increased to $(1 + u_2(t))\varphi$ and the hospitalization rate is adjusted as $(1+u_3(t))k$. Therefore, model (2.1) with an integrated control strategy is given by the following system:

where $\bar{\phi_5} = r_I+d_I+d$. Denote three control measures as $U(t) = (u_1(t), u_2(t), u_3(t))$ and the control set as $\Theta = {\{u_{i}(\cdot)\in (L^{\infty}[0, T], \mathbb{R})}|0 \leq u_i(t) \leq 1, i = 1, 2, 3\}.$ Thus, we define the objective functional as follows:

where the integrand function is given by the following:

The objective of optimal control is to find the optimal control $U^*$ such that the number of infections (asymptomatic and symptomatic) is minimized with the smallest cost of the control measures. Here, $A_1$ and $A_2$ represent the weights of asymptomatic and symptomatic infections, respectively. $B_i$ ($i = 1, 2, 3$) denote the weight coefficients standing for the cost associated with the control variables $u_i(t)$ ($i = 1, 2, 3$), respectively.

In the following, we will analyze the existence of optimal control and its characteristics by Pontryagin's Maximum Principle ^[38,39].

5.1.   Existence of optimal control

Theorem 5.1. There exists an optimal control $U^* = (u_1^*, u_2^*, u_3^*)\in \Theta$ such that $J(U^*) = \min \limits_{U \in \Theta}J(U)$.

Proof. Note that following: (i) the state variables and control variables are non-negative; (ii) the control set $\Theta$ is closed and convex; (iii) the optimal system is bounded, which implies the compactness of the optimal control; (iv) the integrand of the objective functional $J(U)$ is convex on $\Theta$; and (v) there exist constants $a_1, a_2 > 0$ and $q > 1$ such that $J(U)$ satisfies the following:

Therefore, it follows from the results in ^[37] that system (5.1) exists an optimal control $U^* = (u_1^*, u_2^*, u_3^*) \in \Theta$ such that $J(U^*) = \min \limits_{U \in \Theta}J(U)$. □

5.2.   Characterization of the optimal control

To find the characteristic express of optimal control, we denote $X = (S, V, E, A, I, H, R)^T$ and $\lambda = (\lambda_1, ..., \lambda_7)$ and then define the Hamiltonian function as follows:

Here, $\lambda_i$ is called to be adjoint functions. When taking the state system (5.1) with Hamiltonian function (5.2) together, the adjoint system is derived by $d\lambda_i/dt = -\partial\lambda_i/\partial X_i$ ($i = 1, 2..., 7$), that is

The transversality conditions are as follows:

Combining system (5.1), adjoint system (5.3) and transversality conditions (5.4) forms an optimal control problem. By Pontryagin's Maximum Principle, one can have the following theorem.

Theorem 5.2. For any optimal control $U^* \in \Theta$ and solution of system (5.1), the optimal control solution of optimal control problem can be obtained as follows:

where

Proof. The optimal control solution can be obtained by solving the following equations:

Thus, one has the following:

which, together with the upper and the lower bounds of $u_i(t)$ ($i = 1, 2, 3$), derive the characteristic expression of the optimal control solution. □

5.3.   Numerical simulations

All parameter values are shown in Table 1 and the control period is set as 80 days. Furthermore, we assume that the values of the weight coefficients are $A_1 = 20$, $A_2 = 50$, $B_1 = 10$, $B_2 = 20$ and $B_3 = 30$. The initial state variables are chosen as $S(0) = 500,000$, $V(0) = 80,000$, $E(0) = 140$, $A(0) = 40$, $I(0) = 28$, $H(0) = 19$ and $R(0) = 0$.

5.3.1.   Constant control

To better visualise the dependence of peak value and peak timing of vaccinated individuals and symptomatic infections on three constant controls, we designed four different constant values for $u_i$ ($i = 1, 2, 3$): (1) no control $u_1 = u_2 = u_3 = 0$; (2) weak control $u_1 = u_2 = u_3 = 0.1$; (3) moderate control $u_1 = u_2 = u_3 = 0.2$; and (4) strong control $u_1 = u_2 = u_3 = 0.3$. Figure 7 illustrates the effect of the different control intensities on the peak value and peak timing of vaccinated individuals and symptomatic infections, where the numerical results summarized in Table 4. As can be seen in Figure 7 and Table 4, the peak value of vaccinated individuals significantly increases with increasing control intensity, whereas the peak value of symptomatic infections significantly decreases. Additionally, it can be seen that the peak timing of vaccinated individuals and symptomatic infections gradually increases as the control intensity increases. Specifically, for the peak value of vaccinated individuals (symptomatic infections), there is an increase (decrease) about 13.2% (26.1%) for weak control, 28.8% (68.9%) for moderate control, and 47.5% (71.9%) for strong control. In addition, with the increase of control intensity, the peak timing of vaccinated individuals and symptomatic infections are delayed.

5.3.2.   Optimal control

In this section, we combine the forward-backward sweep method ^[40] and the Runge-Kutta scheme to conduct a numerical simulation of the control problem. To investigate the effectiveness of control strategies, we designed four different strategies:

● Strategy A: $(u_1, u_2)$.

● Strategy B: $(u_1, u_3)$.

● Strategy C: $(u_1, u_2, u_3)$.

● Strategy D: Black group (No control).

Strategy A includes the control combination of reducing social distance and improving the vaccination rate. Figure 8(a) shows the optimal curve of Strategy A, where both $u_1$ and $u_2$ are at their maximum values at the beginning, and $u_1$ and $u_2$ decreases from day 64 and day 24, respectively. As shown in Table 5, a comparison between Strategy A and D shows that Strategy A can significantly reduce the peak value of symptomatic infections by 99.91% and shorten the peak timing of symptomatic infections from 20 days to 4 days.

Strategy B combines reducing social distance and improving the hospitalization rate. It follows from Figure 8(b) that both $u_1$ and $u_3$ display the maximum intensity for about 64 days and 38 days, respectively, and then gradually and slowly decrease till they hit 0. Using data from Table 5 to compare Strategy B and D, the results show that Strategy B significantly reduces the peak value of symptomatic infections by 99.93% and shortens the peak timing to 2.4 days.

Strategy C contains three control measures. In the beginning, $u_1$ is at its maximum value for 56 days and then starts to decrease, $u_2$ is at its maximum value only for 18 days and then begins to decrease slowly, and $u_3$ remains at its maximum value for 37 days and then begins to decrease. From Table 5, one can see that Strategy C has the best effect on reducing the peak value and delaying the peak timing of symptomatic infections.

The comparison between the control effects of Strategies A–C is shown in Figure 8(d). Compared to both Strategies A and B, Strategy B has a better effect on reducing the peak value of symptomatic infections. This implies that improving the hospitalization rate is more important than enhancing the vaccination rate. If the effectiveness of the vaccine can be greatly raised and the vaccine immunity waning can be enormously lowered, the situations may undergo completely different changes. Compared to both Strategies A and C, Strategy C adds the control measure of improving the hospitalization rate control measure and has a lower peak value of symptomatic infections. Compared to both Strategies B and C, Strategy C adds the control measure of improving the vaccination rate, in which one can see the peak value of symptomatic infections is reduced but not significantly changed. This shows that increasing the hospitalization rate has a more remarkable impaction on curbing the spread of the epidemic. The more effective vaccine should be quickly developed to improve the role of epidemic prevention. Thus, Strategy C is the most effective strategy to reduce the number of symptomatic infections.

5.3.3.   Cost-effectiveness analysis

The impact of an intervention should not only refer to its effectiveness, but also the corresponding control cost. Hence, to further quantify the benefits of different strategies, we introduce the incremental cost-effectiveness ratio (ICER) ^[41,42], which is defined as $ICER(Y) = \Delta_1 /\Delta_2,$ where $\Delta_1$ = Cost of Strategy Y $-$ Cost of Strategy X, and $\Delta_2$ = Number of infections averted under Strategy Y $-$ Number of infections averted under Strategy X.

Furthermore, the cost of the strategy can be calculated by $\int_0^{t_f} (B_1 u_1 + B_2 u_2+B_3u_3)dt$, and the number of symptomatic infections averted is $\int_0^{t_f} (I-\tilde{I})dt$, where $\tilde{I}$ is the optimal solution under the corresponding strategy. Based on the numerical simulation results, the difference in the number of symptomatic infections averted and the cost are shown in Tables 6 and 7.

A clear observation from Table 6 is that if ICER(A) $<$ ICER(B), then Strategy B is less effective than Strategy A. Moreover, Strategy B has a higher cost and a lower benefit compared to Strategy A. Therefore, Strategy B is excluded.

The next step is a comparison of the ICER values for Strategy A and Strategy C. From Table 7, one can have ICER(A) $<$ ICER(C). Therefore, Strategy C is excluded and Strategy A is more effective.

From the aforementioned comparisons, we can conclude that Strategy A is the most cost-effective among the three strategies.

6.   Conclusions and Discussion

6.1.   Conclusions

To investigate the impacts of vaccinations and varying infectivity across different infection statuses, we developed a mathematical epidemic model and conducted a comprehensive mathematical analysis of its dynamic behaviors and optimal control. This involved establishing the well-posedness of the solution, deriving the expression for the control reproduction number, and demonstrating the existence and global stability of the disease-free equilibrium.

Furthermore, we conducted two types of sensitivity analyses to identify which parameters have the most significant influence on the threshold and the spread of infections. This information serves as the theoretical foundation for devising corresponding control strategies. Notably, the top three parameters found to be most influential are the transmission rate, the vaccination rate, and the hospitalization rate. Consequently, we proposed and rigorously analyzed a time-varying optimal control problem, thus providing specific expressions for optimal control strategies. These strategies hold the potential to assist public health efforts in effectively preventing and controlling epidemics.

Our findings indicated that the control reproduction number plays a pivotal role in determining the global stability of the disease-free equilibrium. The results of the sensitivity analyses are valuable for guiding the control strategy implementation. When controls remain constant, higher control intensities lead to larger final sizes of vaccinated individuals, thus suggesting that a robust control strategy offers the most effective means of reducing the peak value and delaying the peak timing of infections. In the case of time-varying controls, we compared various optimal control combinations in terms of minimizing control costs. The results underscore the effectiveness of an integrated control strategy in curtailing epidemics, with increased vaccine protection rates significantly hastening the suppression of infections. Additionally, we conducted a cost-effectiveness analysis to identify the most economically efficient strategy, which involves a combination of reducing social distance and enhancing vaccination rates, which is distinct from the integrated optimal control strategy. This discovery guides devising strategies that align with limited resources and economic considerations.

6.2.   Discussion

In this section, we will discuss the influence of several key factors. As for model (2.1), it incorporates the natural births and deaths (NBD), while ignoring natural-immunity waning (NIW). If either ignoring NBD or considering NIW, one must discover what changes will occur in the long-term dynamic behavior. To understand the influences of such factors, we set the following four cases to perform the numeral simulations:

● Case 1: Model (2.1)

● Case 2: Model (2.1) without NBD

● Case 3: Model (2.1) without NIW

● Case 4: Model (2.1) without NIW and with NBD

First, we compare Case 1 and Case 2 to explore the influence of NBD on the long-term dynamics. Based on the parameter values in Table 1, we can observe from Figure 9 that NBD takes a slight influence in the time-horizon of an epidemic. To be specific, there are some little changes in the peak value of each compartment. When considering NBD, the peak values of each compartment (except S) are higher than those not considered, which is consistent with the reality for the long-time behavior.

Then, we compare Case 1 and Case 3 to understand how NIW affects the dynamical behavior. Based on the previous work in ^[46], we set the NIW rate as 0.005. Combined with the parameter values in Table 1, we plotted the time series diagrams of each subpopulation. From Figure 10, one can observe that NIW can obviously influence the long-time behavior of susceptible and vaccinated populations, but only has a slight influence on the other populations. It is undeniable that considering NIW within the model will cause non-monotonicity in the proposed system, which will make the analysis on the global dynamical behavior difficult. Our future work will include how to deeply analyze the impact of NIW on the global dynamical behavior.

Finally, we compare Case 1 and Case 4 to investigate the joint influence of NBD and NIW. It follows from Figure 11 that there is an obvious influence on the dynamic behavior of each subpopulation. In the absence of NBD, NIW has the most obvious impact on susceptible and vaccinated individuals, whose numbers significantly increase in the later stage of the epidemic; however, there was a slight influence on the number of other populations (i.e., exposed, asymptomatic, symptomatic, and hospitalized). Furthermore, comparing Figure 10 and Figure 11, one can observe that NBD still has little influence on disease transmission, even in the presence of NIW. This implies that NBD can only slightly affect the long-term dynamic behavior, and perhaps this phenomenon is related with either the proposed model or the parameter values.

Dynamical systems have consistently served as valuable tools in the advancement of life sciences and infectious disease research ^[43,44,45]. Nonetheless, several complex issues warrant deeper consideration. (ⅰ) The formulation of suitable Lyapunov functions and the application of LaSalle's invariant principle demonstrate the global stability of the endemic equilibrium when the control reproduction number exceeds one. (ⅱ) The impact of virus mutations leading to either multiple or super-infections and the modeling of vaccines with varying protection rates against different virus variants. (ⅲ) In the case of diseases like COVID-19, where vaccines may render a significant proportion of individuals as asymptomatic carriers, it is essential to develop a more detailed description and model of this phenomenon. (ⅳ) Acknowledging that both vaccinated and unvaccinated individuals experience natural immune waning after a certain period of recovery from infection, we must explore how secondary infections among the population influence vaccine scheduling, resource allocation, and control strategies. These considerations represent potential directions for future research and analyses.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

We are grateful to the editors and referees for their careful reading and valuable comments which have led to a great improvement of this paper. The authors thank Zhen Jin from Shanxi University for his useful suggestions, and Yiwen Tao from Zhengzhou University for improving English writing. This work was supported by the National Natural Science Foundation of China (Nos. 12126349, 12361102, 11901326, 12071382), Natural Science Foundation of Guizhou Province (No. QianKeHe Jichu-ZK[2021] Yiban002), Guizhou Provincial Department of Education Innovation Group Project (No. Qian Jiao He KY Zi[2019]067).

Conflict of interest

The authors declare there is no conflict of interest.