Next Article in Journal
A Rényi-Type Limit Theorem on Random Sums and the Accuracy of Likelihood-Based Classification of Random Sequences with Application to Genomics
Next Article in Special Issue
Two-Age-Structured COVID-19 Epidemic Model: Estimation of Virulence Parameters through New Data Incorporation
Previous Article in Journal
The Roots of the Reliability Polynomials of Circular Consecutive-k-out-of-n:F Systems
Previous Article in Special Issue
Cumulative Incidence Functions for Competing Risks Survival Data from Subjects with COVID-19
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Assessing the Impact of Time-Varying Optimal Vaccination and Non-Pharmaceutical Interventions on the Dynamics and Control of COVID-19: A Computational Epidemic Modeling Approach

1
School of Mathematics and Data Sciences, Changji University, Changji 831100, China
2
Department of Mathematics, Abdul Wali Khan University Mardan, Mardan 23200, Pakistan
3
Department of Mathematics, University of Peshawar, Peshawar 25120, Pakistan
4
Department of Quantitative Analysis, College of Business Administration, King Saud University, P.O. Box 71115, Riyadh 11587, Saudi Arabia
5
Faculty of Engineering, Future University in Egypt, New Cairo 11835, Egypt
*
Author to whom correspondence should be addressed.
Mathematics 2023, 11(20), 4253; https://doi.org/10.3390/math11204253
Submission received: 30 August 2023 / Revised: 3 October 2023 / Accepted: 4 October 2023 / Published: 11 October 2023

Abstract

:
Vaccination strategies remain one of the most effective and feasible preventive measures in combating infectious diseases, particularly during the COVID-19 pandemic. With the passage of time, continuous long-term lockdowns became impractical, and the effectiveness of contact-tracing procedures significantly declined as the number of cases increased. This paper presents a mathematical assessment of the dynamics and prevention of COVID-19, taking into account the constant and time-varying optimal COVID-19 vaccine with multiple doses. We attempt to develop a mathematical model by incorporating compartments with individuals receiving primary, secondary, and booster shots of the COVID-19 vaccine in a basic epidemic model. Initially, the model is rigorously studied in terms of qualitative analysis. The stability analysis and mathematical results are presented to demonstrate that the model is asymptotically stable both locally and globally at the COVID-19-free equilibrium state. We also investigate the impact of multiple vaccinations on the COVID-19 model’s results, revealing that the infection risk can be reduced by administrating the booster vaccine dose to those individuals who already received their first vaccine doses. The existence of backward bifurcation phenomena is studied. A sensitivity analysis is carried out to determine the most sensitive parameter on the disease incidence. Furthermore, we developed a control model by introducing time-varying controls to suggest the optimal strategy for disease minimization. These controls are isolation, multiple vaccine efficacy, and reduction in the probability that different vaccine doses do not develop antibodies against the original virus. The existence and numerical solution to the COVID-19 control problem are presented. A detailed simulation is illustrated demonstrating the population-level impact of the constant and time-varying optimal controls on disease eradication. Using the novel concept of human awareness and several vaccination doses, the elimination of COVID-19 infections could be significantly enhanced.

1. Introduction

The world’s economy, communities, and public health have been profoundly impacted by the emergence of the novel coronavirus, SARS-CoV-2, which led to the COVID-19 pandemic. The first case of infection with this virus appeared in December 2019 in Wuhan, Hubei province of China [1]. This disease quickly spread outside of China only a few weeks after it first appeared. Initial reports of the disease outside China came from Japan and Thailand [2]. The COVID-19 pandemic has had a significant impact on society in different aspects. From a global health perspective, the virus has led to millions of worldwide cases of infection and loss of life. It has strained healthcare systems, particularly in regions with limited resources and capacity. This global disease also has far-reaching impacts on the economy, with disruptions in worldwide supply chains, business closures, and job losses. Socially, COVID-19 has led to school closures, travel restrictions, and changes in social behaviors, affecting individuals’ mental health and well-being [3]. Extensive scientific investigation has been carried out to comprehend the dynamics and transmission of the SARS-CoV-2 virus [4]. There are many different indications and symptoms that COVID-19 might exhibit. On the other hand, some patients may remain symptom-free, and some may experience mild and then severe disease symptoms. Fever, coughing, sore throat, shortness of breath, feelings of exhaustion or low energy, muscle discomfort, and body pains are among the main signs and symptoms of COVID-19 [5,6]. The presence and intensity of symptoms can differ from person to person, and some people may not experience any symptoms at all or only experience minor ones. Different countries employ various techniques to stop the spread of infections, and the majority of countries adhere to similar rules such as social distancing, isolation, and self-quarantine [7].
Vaccination remains one of the effective interventions against severe infectious diseases [8,9,10]. The careful application of non-pharmaceutical therapies and vaccination programs work together to reduce the spread of COVID-19 globally. To prevent the disease, multiple dosages of vaccines are given. The majority of vaccinations are given to a person as primary (first time), secondary (second time), and booster shots (third time). The first course of vaccination doses administered to people who have never had any vaccination doses before is referred to as primary vaccination. The subsequent dosage given after the original immunization is referred to second-time vaccine and is also referred to as a second dose. The third-time vaccine is an additional dosage administered to people who have finished their first vaccine course to strengthen and prolong their immunological response [7].
In order to investigate the complex dynamics of infectious disease, different methodologies have been developed. The implementation of mathematical models is a valuable tool that has been utilized successfully to present different aspects of infectious diseases. Usually, these models include classical (ordinary and partial) derivatives [11,12,13], stochastic derivatives [14,15] and fractional derivatives [16,17]. In particular, to better explore the dynamic aspects of COVID-19, several epidemic models have been developed [18,19,20]. The impact of face mask use by the general public to curtail the COVID-19 pandemic was studied in [21]. Augusto et al. studied the changing behavior of the COVID-19 model [22]. The global impact of the first year of COVID-19 vaccination programs was studied in [23]. The impact of vaccination on two variants of COVID-19, alpha and delta, was studied in [24]. Ngonghala et al. [25] considered the omicron and delta variants of COVID-19 in the presence of multiple vaccinations. Their study revealed that treatment leads to a reduction in hospitalization rates, and the potential for COVID-19 elimination is increased when investments in control resources are directed toward promoting mask usage and vaccine intervention. In [26], the authors analyzed the mitigation of the pandemic via double-dose vaccination using an epidemic modeling approach. The outcomes of their study indicated that primary and secondary vaccination alone is not adequate for infection reduction, and thus it is essential to provide the booster shot (third time) of the vaccine for a better eradication of the infection. Recently, a similar fractional and fractal-fractional study with an exponential-type kernel analyzing the dynamics of COVID-19 under vaccination was presented in [27].
In this paper, we develop a mathematical model based on the SEVIHR type of compartmental model and incorporate three vaccine compartments for the first, second, and booster shots of COVID-19 vaccines. This study stands out from previous literature by considering multiple vaccination compartments in relation to the rates of antibody production against the original virus. We also used the assumption that vaccinated individuals might become infected with the virus if the antibodies have not been developed even after vaccinating to account for breakthrough infections. This enabled us to illustrate the reproductive number causing a significant widespread issue of the disease. The outcomes of these simulations aid those who are skeptical of vaccination in making thoughtful decisions. The section-wise description of the present work is as follows: The formulation of the COVID-19 model is presented in Section 2. Section 3 covers the basic qualitative properties and results of the proposed model. Section 4 presents the role of basic reproduction numbers, vaccination coverage, and bifurcation analysis. A simulation of the proposed model with constant control measures is given in Section 5. The sensitivity analysis is performed in Section 6. The formulation of the optimal control problem is given in Section 6. A simulation of the optimal control problem estimating the optimal solution is presented in Section 8. Finally, Section 9 presents the conclusion.

2. Modeling the Dynamics of COVID-19 with Multiple Vaccine Doses

We divide the total population into the following compartments: 1. susceptible individuals S; 2. exposed individuals E; 3. individuals with first-time vaccination V 1 ; 4. individuals with second-time vaccination V 2 ; 5. individuals with booster shots V 3 ; 6. infected individuals I; 7. hospitalized individuals I H ; and 8. R represents the recovered individuals. Therefore, the entire population can be expressed as
N ( t ) = S ( t ) + E ( t ) + V 1 ( t ) + V 2 ( t ) + V 3 ( t ) + I ( t ) + I H ( t ) + R ( t ) .
The group of susceptible people is generated as a result of the birth rate θ , which reduces due to the transmission to the vaccinated class V 1 upon receiving the first vaccine dose at a rate ξ 1 . Further, this class experiences a decrease after becoming infected at the contact rate  α . All population groups experience natural mortality at a rate of μ . Thus, we obtain the following differential equation.
S t = θ α I N S ζ 1 + μ S ,
Individuals in the exposed class are generated as a result of effective contacts between individuals in the infectious class I with those in classes S, V 1 , V 2 , and V 3 at the contact rates α . The symbols δ 1 , δ 2 , and δ 3 are the respective probabilities that vaccine recipients in the V 1 , V 2 and V 3 groups do not develop antibodies to the original viruses after 28 days of inoculation. Therefore, the individuals susceptible to the original viruses considered in the model are S, δ 1 V 1 , δ 2 V 2 and δ 3 V 3 . The transmission rate of the exposed individuals to the infectious class is denoted by σ and as a result, we derive the following equation for this class:
E t = α I S N + α δ 1 I V 1 N + α δ 2 I V 2 N + α δ 3 I V 3 N σ + μ E ,
Individuals in the susceptible class move to the class of first-time vaccinated individuals after receiving the initial vaccination at a rate of ζ 1 . Subsequently, the individuals in this group reduce due to the contact rate α with infectious individuals and the administration of the second dose at the rate ζ 2 . We obtain the following mathematical form for the dynamics of the first-time vaccinated class.
V 1 t = ζ 1 S α δ 1 I V 1 N ζ 2 + μ V 1 ,
The population of second-time-vaccinated individuals is initially formed by administering a second dose to individuals in V 1 class at a rate ζ 2 . This group experiences a decrease due to the contact rate α with infectious individuals and the administration of booster shots at the rate ζ 3 . We obtain the following equation for the dynamics of the second-dose-vaccinated class.
V 2 t = ζ 2 V 1 α δ 2 I V 2 N ζ 3 + μ V 2 .
The class of third-time-vaccinated individuals is initiated through by administrating the booster shot to individuals in the V 2 class at a rate of ζ 3 . This population experiences a decrease due to several factors, such as the contact rate α with infected people and the natural mortality rate μ . Thus, we obtain the following differential equation.
V 3 t = ζ 3 V 2 α δ 3 I V 3 N μ V 3 .
The exposed individuals become infected and move to I class at the rate σ . The population in this class declined because of the hospitalization at a rate η , and the COVID-19 and natural deaths at rates denoted by γ and μ , respectively. Thus, we derive the following expression.
I t = σ E η + γ + μ I .
The fraction 0 < b < 1 of the η I moves to R because of their natural immunity, while the remaining ( 1 b ) are hospitalized. The individuals in the I H class convalesce after proper treatment and join the recovered class at the recovery rate d. Thus, we obtain the following equations for the dynamics of hospitalized and recovered individuals.
I H t = η 1 b I d + γ + μ I H , R t = η b I + d I H μ R .
In the result of the above discussion, the compartmental model for COVID-19 transmission is summarized as follows:
S t = θ α I N S ζ 1 + μ S , E t = α I S N + α δ 1 I V 1 N + α δ 2 I V 2 N + α δ 3 I V 3 N σ + μ E , V 1 t = ζ 1 S α δ 1 I V 1 N ζ 2 + μ V 1 , V 2 t = ζ 2 V 1 α δ 2 I V 2 N ζ 3 + μ V 2 , V 3 t = ζ 3 V 2 α δ 3 I V 3 N μ V 3 , I t = σ E γ + η + μ I , I H t = η 1 b I μ + γ + d I H , R t = η b I + d I H μ R .
subject to non-negative initial conditions, S ( 0 ) = S 0 , E ( 0 ) = E 0 , V 1 ( 0 ) = V 1 , V 2 ( 0 ) = V 2 , V 3 ( 0 ) = V 3 , I ( 0 ) = I 0 , I H ( 0 ) = I H 0 , R ( 0 ) = R 0 . Let
λ = ( α S + α δ 1 V 1 + α δ 2 V 2 + α δ 3 V 3 ) N ,
and
q 1 = ( ζ 1 + μ ) , q 2 = ( σ + μ ) , q 3 = ( ζ 2 + μ ) , q 4 = ( ζ 3 + μ ) , q 5 = ( η + γ + μ ) , q 6 = ( d + γ + μ ) .
Then, (1) becomes
S t = θ α I N S q 1 S , E t = λ I q 2 E , V 1 t = ζ 1 S α δ 1 I N V 1 q 3 V 1 , V 2 t = ζ 2 V 1 α δ 2 I N V 2 q 4 V 2 , V 3 t = ζ 3 V 2 α δ 3 I N V 3 μ V 3 , I t = σ E q 5 I , I H t = η 1 b I q 6 I H , R t = η b I + d I H μ R .
The transmission between various compartments are shown in the Figure 1.

3. Qualitative Analysis of the Model

This section presents necessary mathematical aspects of the proposed COVID-19 transmission model having multiple vaccine compartments. We proceed as follows:

3.1. Positivity and Boundedness

3.1.1. Positivity

Theorem 1. 
Given S 0 0 , E 0 0 , V 1 0 0 , V 2 0 0 , V 3 0 0 , I 0 0 , I H 0 0 , R 0 0 . Then, the solution ( S ( t ) , E ( t ) , V 1 ( t ) , V 2 ( t ) , V 3 ( t ) , I ( t ) , I H ( t ) , R ( t ) ) of model (2) are positive for all t > 0 .
Proof. 
First equation of the model (2) gives
S t = θ α I N S q 1 S , S t α ¯ q 1 S , α ¯ ( t ) = α I N .
The integrating factor is given by exp 0 t α ¯ s d s + q 1 t , multiplying the inequality (3) by the integrating factor, we obtain,
S t exp ( 0 t α ¯ s d s + q 1 t ) 0 .
The solution of (4) implies
S t S 0 exp ( ( 0 t α ¯ s d s + q 1 t ) ) > 0 .
In similar pattern, it can be shown that E 0 , V 1 0 , V 2 0 , V 3 0 , I 0 , I H 0 , R 0 in the model (2). This implies that ( S ( t ) , E ( t ) , V 1 ( t ) , V 2 ( t ) , V 3 ( t ) , I ( t ) , I H ( t ) , R ( t ) ) are all non-negative for non-negative initial conditions.    □

3.1.2. Boundedness

Theorem 2. 
The solution ( S ( t ) , E ( t ) , V 1 ( t ) , V 2 ( t ) , V 3 ( t ) , I ( t ) , I H ( t ) , R ( t ) ) of the model (2) are bounded.
Proof. 
Adding all the equations of the proposed model (2), we have
N t = θ μ N γ I + I H , N t θ μ N .
Solving the inequality in (5), we have the following steps.
d e μ t N t θ e μ t d t , e μ t N t θ μ e μ t + C , N t e μ t θ μ e μ t + C .
Using the initial condition N ( t ) = N ( 0 ) , at t = 0 implies that
N t e μ t θ μ e μ t + N 0 θ μ , e μ t N 0 + θ μ 1 e μ t , lim t N t θ μ .
Hence, N ( t ) is bounded by θ μ , and by using the comparison theorem [28], we deduce that N ( t ) θ μ , if N ( 0 ) θ μ . Therefore, θ μ remains the upper bound of the region Ω = S t , E t , V 1 t , V 2 t , V 3 t , I t , I H t , R t : S t + E t + V 1 t + V 2 t + V 3 t + I t + I H t + R t θ μ . . Thus, the region Ω is also positively invariant [29].    □

3.2. Equilibria and the Threshold Parameter of the Model

The COVID-19 compartmental model (2) exhibits two equilibria, namely, the COVID-free equilibrium ( CFE ) , denoted as P 0 . This equilibrium can be given by the following expressions:
P 0 = S 0 , E 0 , V 1 0 , V 2 0 , V 3 0 , I 0 , I H 0 , R 0 = θ q 1 , 0 , θ ζ 1 q 1 q 3 , θ ζ 1 ζ 2 q 1 q 3 q 4 , θ ζ 1 ζ 2 ζ 3 μ q 1 q 2 q 3 , 0 , 0 , 0 .
Taking the next generation approach into account [30], the threshold number R 0 is formulated as follows:
R 0 = α σ q 2 q 5 N 0 S 0 + δ 1 V 1 0 + δ 2 V 2 0 + δ 3 V 3 0 .

3.3. Local Stability at COVID-Free Equilibrium Point

Theorem 3. 
The CFE point of the proposed model (2) is locally asymptotically stable (LAS) if R 0 < 1 and unstable when R 0 > 1 .
Proof. 
To establish the desired outcome, we need to demonstrate that the Jacobian of the linearized system at the CFE has negative eigenvalues. The matrix J ( P 0 ) can be evaluated as
J P 0 = q 1 0 0 0 0 α μ q 1 0 0 0 q 2 0 0 0 α μ q 1 + α μ δ 1 ζ 1 q 1 q 3 + α μ δ 2 ζ 1 ζ 2 q 1 q 3 q 4 + α δ 3 ζ 1 ζ 2 ζ 3 q 1 q 3 q 4 0 0 ζ 1 0 q 3 0 0 α μ δ 1 q 1 q 3 0 0 0 0 ζ 2 q 4 0 α μ δ 2 ζ 1 ζ 2 q 1 q 3 q 4 0 0 0 0 0 ζ 3 μ α δ 3 ζ 1 ζ 2 ζ 3 q 1 q 3 q 4 0 0 0 σ 0 0 0 q 5 0 0 0 0 0 0 0 1 b η q 6 0 0 0 0 0 0 b η d μ .
The eigenvalues of J ( P 0 ) are μ , μ , q 1 , q 3 , q 4 , q 6 , and the roots of the 2nd degree equation c 0 λ 2 + c 1 λ + c 2 = 0 . Here, the coefficients are defined in terms of the basic reproduction number R 0 as
c 0 = 1 , c 1 = q 2 + q 5 , c 2 = q 2 q 5 1 R 0 .
Since, c 0 , c 1 > 0 and c 2 > 0 if R 0 < 1 . Therefore, the well-known Ruth–Hurwitz criteria demonstrate that the system (2) is LAS.    □

3.4. Global Asymptotic Stability of (GAS) for Special Case

The GAS of the model at CFE point is provided for a special case. Consider the special case of the model (2) with γ = δ 1 = δ 2 = δ 3 = 0 . The assumption regarding the subsequent parameters is made to facilitate mathematical analysis. The feasible region for the special case of the model (2) is constructed as
Ω * = S , E , V 1 , V 2 , V 3 , I , I H , R R + * : S S 0 , V 1 V 1 0 , V 2 V 2 0 , V 3 V 3 0 .
It can be shown that Ω * is a positive invariant set and attracts the solution of system (2) with respect to the special case. Moreover, the basic reproductive number of the reduced model is as follows:
R ^ 0 = R 0 γ = δ 1 = δ 2 = δ 3 = 0 = α σ S 0 q 1 q 2 q 5 * N 0 , w h e r e q 5 * = η + μ , q 6 * = μ + d .
Theorem 4. 
The CFE of the special case of the model (2) with γ = δ 1 = δ 2 = δ 3 = 0 is GAS in Ω * whenever, R ^ 0 1 .
Proof. 
L t = A 1 S S 0 2 2 S 0 + A 2 E + A 3 I , L = A 1 S S 0 S 0 S + A 2 E + A 3 I , L = A 1 S S 0 S 0 α I N S S 0 q 1 S S 0 α I N S 0 + A 2 α I S N q 2 E + A 3 σ E q 5 * I A 1 α S S 0 2 S 0 I + A 2 α I S 0 N 0 q 2 E + A 3 σ E q 5 * I , = A 1 α S S 0 2 S 0 I + σ A 3 q 2 A 2 E + α S 0 N 0 A 2 q 5 * A 3 I , = A 1 α S S 0 2 S 0 I + σ A 3 q 2 A 2 E + q 5 * A 3 α S 0 q 5 * A 3 N 0 A 2 1 I , L S S 0 2 I + q 2 R ^ 0 1 E ,
where,
A 1 = S 0 α , A 2 = σ A 3 q 2 , A 3 = 1 .
Thus, L ( t ) 0 if R ^ 0 1 further L ( t ) = 0 if E = I = I H = 0 . So, the number of infected individuals becomes zero as t . Using E = I = I H = 0 in the above model, we have S θ q 1 , V 1 θ ζ 1 q 1 q 3 , V 2 θ ζ 1 ζ 2 q 1 q 3 q 4 , V 3 θ ζ 1 ζ 2 ζ 3 μ q 1 q 3 q 4 and R 0 as t . Thus, using Lyapunov stability theorem, every solution of the given model with non-negative initial condition approaches to P 0 as t in Ω . Thus, it follows that the system (2) is GAS.    □

4. Interpretation of Vaccination Coverages Based on R 0

We provide the critical rate of vaccine coverage that could lead to eradication of the infection, denoted as R 0 ( ζ 1 , ζ 2 , ζ 3 ) = R 0 e . In a scenario with no vaccine, i.e., when ζ 1 = 0 , ζ 2 = 0 , and ζ 3 = 0 , then R 0 is reduced to
R 0 e = R 0 0 , 0 , 0 = α σ q 2 q 5 .
After some rearrangement R 0 ( ζ 1 ) and with ζ 2 = ζ 3 = 0 can be written as
R 0 ( ζ 1 ) = α σ ζ 1 δ 1 + μ q 2 q 5 ζ 1 + μ .
Thus,
R 0 = lim ζ 1 R 0 ζ 1 = α σ q 2 q 5 δ 1 = R 0 e δ 1 .
The partial derivative with respect to ζ 1 leads to the following form:
R 0 ζ 1 = R 0 e μ 1 δ 1 q 1 2 < 0 .
Therefore, R 0 e δ 1 R 0 ( ζ 1 ) R 0 e , and hence, R 0 e < 1 implies R 0 ( ζ 1 ) < 1 . Further, if R 0 e > 1 , it is important to note that
R 0 < 1 R 0 δ 1 < 1 δ 1 > δ 1 * = 1 R 0 e .
This interpretation suggests that when the value of δ 1 is low, and if the effective reproductive number R 0 e > 1 , the disease might not spread extensively provided that the first-dose vaccination coverage is substantial. Similarly, if ζ 1 , ζ 2 0 ,
R 0 ζ 1 , ζ 2 = R 0 e μ 2 + μ δ 1 ζ 1 + μ + δ 2 ζ 1 ζ 2 ζ 1 + μ ζ 2 + μ .
After some rearrangement, we have
R 0 ζ 1 , ζ 2 = 1 μ q 3 R 0 ζ 1 , 0 + R 0 e μ + δ 2 ζ 1 ζ 2 q 1 q 3 .
Thus,
R 0 ζ 1 , = lim ζ 2 R 0 ζ 1 , ζ 2 = R 0 e μ + δ 2 ζ 1 q 1 .
The partial differentiation with respect to ζ 2 yields
R 0 ζ 2 = R 0 e δ 1 δ 2 ζ 1 μ q 1 q 2 q 3 2 q 5 < 0 .
Therefore, R 0 e ( μ + δ 2 ζ 1 ) q 1 R 0 ( ζ 1 , ζ 2 ) R 0 e , and hence R 0 e < 1 implies R 0 ( ζ 1 , ζ 2 ) < 1 . Further if R 0 e > 1 , it is important to note that
R 0 ζ 1 , < 1 R 0 e μ + δ 2 ζ 1 q 1 < 1 δ 2 > δ 2 * = q 1 R 0 e μ ζ 1 .
The above interpretation implies that when the value of δ 2 is minimal and if the effective reproductive number R 0 e > 1 , there is a possibility of infection elimination if the second vaccination coverage is substantial. Finally, in a scenario if not only first and second vaccine doses are administrated, but additionally, a booster short is provided, i.e., when ζ 1 , ζ 2 , ζ 3 0 then,
R 0 ζ 1 , ζ 2 , ζ 3 = μ R 0 ζ 1 , ζ 2 , 0 q 4 + R 0 e μ 2 + μ δ 1 ζ 1 + μ + δ 3 ζ 1 ζ 2 ζ 3 q 1 q 3 q 4 .
Thus,
R 0 ζ 1 , ζ 2 , = lim ζ 3 R 0 ζ 1 , ζ 2 , ζ 3 = R 0 e μ 2 + μ δ 1 ζ 1 + μ + δ 3 ζ 1 ζ 2 q 1 q 3 .
Taking partial derivative with respect to ζ 3 yields
R 0 ζ 3 = α σ μ δ 2 δ 3 ζ 1 ζ 2 q 1 q 2 q 3 q 4 2 q 5 < 0 .
Therefore, R 0 e μ 2 + μ δ 1 ζ 1 + μ + δ 3 ζ 1 ζ 2 q 1 q 3 R 0 ζ 1 , ζ 2 , ζ 3 R 0 e , and hence R 0 e < 1 implies R 0 ( ζ 1 , ζ 2 , ζ 3 ) < 1 . Further, if R 0 e > 1 it is important to note that
R 0 ζ 1 , ζ 2 , < 1 R 0 e μ 2 + μ δ 1 ζ 1 + μ + δ 3 ζ 1 ζ 2 q 1 q 3 < 1 δ 3 > δ 3 * = q 1 q 3 R 0 e ζ 1 μ μ + δ 1 ζ 1 + ζ 2 ζ 1 .
This implies that when δ 3 is minimal and the effective reproductive number R 0 e > 1 , the disease has the potential to be eliminated through the administration of a booster shot to a person who had already received the initial dose.

Backward Bifurcation Analysis

In this subsection, we discuss the existence of backward bifurcation of the system (2). To analyze this, we follow the center manifold theory. If we take α as bifurcation parameter then at R 0 = 1 we have,
α * = q 1 q 2 q 3 q 4 q 5 σ μ q 3 q 4 + ζ 1 μ q 4 δ 1 + ζ 2 μ δ 2 + δ 3 ζ 3 .
The variables in the system (2) are subjected to the following variations so that S = x 1 , E = x 2 , V 1 = x 3 , V 2 = x 4 , V 3 = x 5 , I = x 6 , I H = x 7 , and R = x 8 . Further, using the vector notation x = ( x 1 , x 2 , x 3 , x 4 , x 5 , x 6 , x 7 , x 8 ) t r . COVID-19 (2) can be written equivalently as d x d t = f , where f = ( f 1 , f 2 , f 3 , f 4 , f 5 , f 6 , f 7 , f 8 ) t r as shown below:
x 1 t = θ α x 6 N x 1 q 1 x 1 , x 2 t = α x 6 N x 1 + α δ 1 x 6 N x 3 + α δ 2 x 6 N x 4 + α δ 3 x 6 N x 5 q 2 x 2 , x 3 t = ζ 1 x 1 α δ 1 x 6 N x 3 q 3 x 3 , x 4 t = ζ 2 x 3 α δ 2 x 6 N x 4 q 4 x 4 , x 5 t = ζ 3 x 4 α δ 3 x 6 N x 5 μ x 5 , x 6 t = σ x 2 q 5 x 6 , x 7 t = η 1 b x 6 q 6 x 7 , x 8 t = η b x 6 + d x 7 μ x 8 ,
where N = i = 1 8 x i . The Jacobian matrix evaluated at CFE point P 0 with α * is
J P 0 = q 1 0 0 0 0 α * μ q 1 0 0 0 q 2 0 0 0 α * μ q 1 + α * μ δ 1 ζ 1 q 1 q 3 + α * μ δ 2 ζ 1 ζ 2 q 1 q 3 q 4 + α * δ 3 ζ 1 ζ 2 ζ 3 q 1 q 3 q 4 0 0 ζ 1 0 q 3 0 0 α * μ δ 1 ζ 1 q 1 q 3 0 0 0 0 ζ 2 q 4 0 α * μ δ 2 ζ 1 ζ 2 q 1 q 3 q 4 0 0 0 0 0 ζ 3 μ α * δ 3 ζ 1 ζ 2 ζ 3 q 1 q 3 q 4 0 0 0 σ 0 0 0 q 5 0 0 0 0 0 0 0 η 1 b q 6 0 0 0 0 0 0 b η d μ ,
The Jacobian matrix has a simple eigenvalue calculated at α * . The right and left eigenvectors are denoted by W = ( w 1 , w 2 , w 3 , w 4 , w 5 , w 6 , w 7 , w 8 ) and V = ( v 1 , v 2 , v 3 , v 4 , v 5 , v 6 , v 7 , v 8 ) , respectively, where
w 1 = μ q 2 q 3 q 4 q 5 w 6 σ q 1 μ q 3 q 4 + ζ 1 μ δ 1 q 4 + ζ 2 μ δ 2 + δ 3 ζ 3 , w 2 = q 5 w 6 σ , w 3 = w 6 μ q 2 q 5 δ 1 ζ 1 q 4 σ q 3 μ q 3 q 4 + ζ 1 μ δ 1 q 4 + ζ 2 μ δ 2 + δ 3 ζ 3 + μ q 2 q 5 ζ 1 q 3 q 4 σ q 1 q 3 μ q 3 q 4 + ζ 1 μ δ 1 q 4 + ζ 2 μ δ 2 + δ 3 ζ 3 , w 4 = μ q 2 q 5 w 6 ζ 1 ζ 2 δ 1 q 1 q 4 + q 2 μ + δ 2 q 1 + ζ 3 σ q 1 q 3 q 4 μ 2 q 4 + μ δ 1 ζ 1 q 4 + ζ 2 μ μ + δ 2 ζ 1 + μ + δ 3 ζ 1 ζ 3 , w 5 = q 2 q 5 w 6 ζ 1 ζ 2 ζ 3 μ δ 1 q 1 q 4 + q 3 μ δ 2 q 1 + μ + δ 3 q 1 q 4 μ σ q 1 q 3 q 4 μ 2 q 4 + μ δ 1 ζ 1 q 4 + ζ 2 μ μ + δ 2 ζ 1 + μ + δ 3 ζ 1 ζ 3 , w 6 > 0 , w 7 = η 1 b w 6 q 6 , w 8 = d η + b γ η + b η μ μ q 6 w 6 ,
and,
v 1 = 0 , v 2 > 0 , v 3 = v 4 = v 5 = 0 , v 6 = q 2 v 2 σ , v 7 = v 8 = 0 .
Furthermore, bifurcation coefficients a ¯ , and b ¯ of the proposed model (2) evaluated at ( P 0 , α * ) are calculated as
a ¯ = k , i , j = 1 8 v k w i w j 2 f k P 0 , α * x i x j , b ¯ = k , i = 1 8 v k w i 2 f k P 0 , α * x i β .
So we have
a ¯ = 2 α * μ v 2 w 6 θ q 1 q 3 q 4 q 3 q 4 ζ 1 w 1 + μ w 2 + w 3 + w 4 + w 5 + w 6 + w 7 + w 8 q 1 w 3 δ 1 + w 4 δ 2 + w 5 δ 3 + w 1 + w 2 + w 3 + w 4 + w 5 + w 6 + w 7 + w 8 ζ 1 μ q 4 δ 1 + ζ 2 μ δ 2 + ζ 3 δ 3 < 0 , b ¯ = v 2 w 6 μ q 3 q 4 + ζ 1 μ q 4 δ 1 + ζ 2 μ δ 2 + δ 3 ζ 3 q 1 q 3 q 4 > 0 .
The presence of the backward bifurcation in the system depends on the sign of b ¯ , as indicated in [31]. Specifically, if a ¯ > 0 and b ¯ > 0 , the system will experience a backward bifurcation. In the context of backward bifurcation, the endemic equilibrium and the CFE coexist with one being stable and the other being unstable when R 0 < 1 . For the COVID-19 model (2), the biological significance of backward bifurcation lies in the fact that the condition R 0 < 1 is essential, but not solely adequate for mitigating the infection from the community. The outcome depends on the initial population size.

5. Estimating the Time Series Solution

This section focuses on the numerical simulation of the epidemic model (2) to analyze the dynamics of various state variables for R 0 > 1 and R 0 < 1 . The model without optimal controls is numerically solved using the well-known Runge–Kutta iterative scheme of fourth order. Simulation is conducted in Matlab version R2022b for the 0–700 days. Initially, the model is simulated with the parameter values given in Table 1 such that R 0 < 1 , and the resulting plots are depicted in Figure 2 and Figure 3. The figures show that the population curves converge to the COVID-19-free equilibrium state. Similarly, in Figure 4 and Figure 5, a few values of the model’s parameters provided in Table 1 are changed in a way that R 0 > 1 . These graphical interpretations show that the solution trajectories converge to an endemic state. Moreover, we analyzed the dynamics of the model with different initial values of state variables for R 0 > 1 and R 0 < 1 as shown in the Figure 6 and Figure 7, respectively.
In the next section, the application of optimal control theory for the mitigation of pandemic is carried out. Before establishing effective optimal controls, we perform the sensitivity analysis to investigate the most influential factors on the transmission of disease.

6. Sensitivity Analysis

Sensitivity theory is a powerful tool that provides insights into the model’s behavior and helps in better understanding the influence of the variations in the input parameters on the respective output of the model. This technique is used in various fields, including economics, engineering, and biological sciences, to assess the potential factors of the problem under consideration. The analysis examines how small changes in the input values of the model affect the corresponding output. In the epidemiological modeling approach, sensitivity analysis provides valuable insights into identifying the parameters that play a substantial role in influencing disease transmission and control. In this study, we specifically performed sensitivity analysis of some crucial parameters. We employ the well-known normalized parametric scheme based on the forward sensitivity index of the model’s parameters, as described by Chitnis et al. [34]. A positive (or negative) index illustrates that the parameter has a direct (or inverse) effect on R 0 .
Definition 1. 
The normalized sensitivity index, which is used to quantify the relative change in R 0 concerning changes in model’s various parameters, is defined as
Y x = x R 0 × R 0 x .
Applying Equation (7), the Table 2 shows the normalized sensitivity index assessing the proportional change of R 0 with respect to the model parameters. The parameters α , σ , μ , δ 1 , δ 2 , and δ 3 directly effect R 0 . This shows that the value of R 0 will grow with an increase in the parameters above (or decrease respectively). The parameters γ , η , ζ 1 , ζ 2 , ζ 3 have an inverse relationship with R 0 . The graphical representation of sensitivity indices are depicted in the plot (Figure 8).
The significant insights gained from the sensitivity analysis can be helpful in various applications such as
  • Identifying potential parameters: This helps identify the input parameter(s) that play a substantial role in influencing disease dynamics.
  • This helps understand how uncertainties in input parameters can propagate to uncertainties in the model’s predictions.
  • Optimization: In optimization problems, the sensitivity indices of model parameters can guide the setting of appropriate optimal interventions.
  • Decision-making: Understanding the sensitivity of the model to various inputs can assist decision-makers in making informed choices.

7. Optimal Control Analysis of COVID-19 Model

In this section, we introduce an optimal control strategy for the model (2). The aim of optimal control is to minimize the COVID-19 disease by providing the following control strategies.
  • u 1 ( t ) : The first control is used to reduce the number of effective contacts between infected and susceptible individuals.
  • u 2 ( t ) : The second control is used to enhance the first vaccine dose efficacy.
  • u 3 ( t ) : The third control is used to reduce the possibility that after 28 days of vaccination, those who received the first dose do not develop immunity to the original virus.
  • u 4 ( t ) : This control is used to reduce the possibility that after 28 days of vaccination, those who received the second dose do not develop immunity to the original virus.
  • u 5 ( t ) : The fifth control is used to reduce the possibility that after 28 days of vaccination, those who received the third dose do not develop immunity to the original virus.
  • u 6 ( t ) : The sixth control is used to enhance the second-time vaccination rate.
  • u 7 ( t ) : This control is used to enhance the third-time vaccination rate.
By employing the controls mentioned above, this section formulates an optimal control problem that elucidates how time-dependent control strategies contribute to disease eradication. The control model is established in (8). Based on the sensitivity index, the desired controls are selected. Hence, the resulting control model is structured as follows:
S t = θ α I N S ( 1 u 1 ( t ) ) u 2 ( t ) + μ S , E t = α I S N ( 1 u 1 ( t ) ) + α δ 1 I V 1 N ( 1 u 3 ( t ) ) + α δ 2 I V 2 N ( 1 u 4 ( t ) ) + α δ 3 I V 3 N ( 1 u 5 ( t ) ) σ + μ E , V 1 t = u 2 ( t ) S α δ 1 I V 1 N ( 1 u 3 ( t ) ) u 6 ( t ) + μ V 1 , V 2 t = u 6 ( t ) V 1 α δ 2 I V 2 N ( 1 u 4 ( t ) ) u 7 ( t ) + μ V 2 , V 3 t = u 7 ( t ) V 2 α δ 3 I V 3 N ( 1 u 5 ( t ) ) μ V 3 , I t = σ E η + μ + γ I , H t = η 1 b I γ + μ + d I H , R t = η b I + d I H μ R .
with the same initial conditions given in (2). The respective objective functional is described as
J u 1 , u 2 , u 3 , u 4 , u 5 , u 6 , u 7 = 0 T A 1 E + A 2 V 1 + A 3 V 2 + A 4 V 3 + A 5 u 1 2 2 + A 6 u 2 2 2 + A 7 u 3 2 2 + A 8 u 4 2 2 + A 9 u 5 2 2 + A 10 u 6 2 2 + A 11 u 7 2 2 d t ,
where, A 1 , A 2 , A 3 , A 4 are the balancing constants associated with the suggested variables of the objective function, A 5 , A 6 , A 7 , A 8 are the the cost factors while T represents the final time. We used the quadratic objective functional due to the non-linearity of the intervention considered for the mitigation of the pandemic. For a more comprehensive understanding, please refer to the work and associated references [35,36,37]. Our primary goal is to identify the optimal controls
u ¯ i t f o r i = 1 , 2 , 3 , 4 , 5 , 6 , 7 ,
so that,
J u ¯ 1 , u ¯ 2 , u ¯ 3 , u ¯ 4 , u ¯ 5 , u ¯ 6 , u ¯ 7 = min Ξ J u 1 , u 2 , u 3 , u 4 , u 5 , u 6 , u 7 .
The corresponding control set is given by
Ξ = u 1 , u 2 , u 3 , u 4 , u 5 , u 6 , u 7 : [ 0 , T ] 0 , 1 u 1 , u 2 , u 3 , u 4 , u 5 , u 6 , u 7 i s a L e b e s g u e m e a s u r a b l e .
The Lagrangian and Hamiltonian for the provided control system (8) are shown as L and H , respectively, and are given as follows:
L = A 1 E + A 2 V 1 + A 3 V 2 + A 4 V 3 + 1 2 A 5 u 1 2 + A 6 u 2 2 + A 7 u 3 2 + A 8 u 4 2 + A 9 u 5 2 + A 10 u 6 2 + A 11 u 7 2 ,
and
H = A 1 E + A 2 V 1 + A 3 V 2 + A 4 V 3 + 1 2 A 5 u 1 2 + A 6 u 2 2 + A 7 u 3 2 + A 8 u 4 2 + A 9 u 5 2 + A 10 u 6 2 + A 11 u 7 2 + λ 1 θ α I S N 1 u 1 u 2 + μ S + λ 2 α I N S 1 u 1 + V 1 δ 1 1 u 3 + V 2 δ 2 1 u 4 + V 3 δ 3 1 u 5 σ + μ E + λ 3 u 2 S α I V 1 δ 1 N 1 u 3 u 6 + μ V 1 + λ 4 V 1 u 6 I V 2 δ 2 α N 1 u 4 u 7 + μ V 2 + λ 5 u 7 V 3 α I V 3 δ 3 N 1 u 5 μ V 3 + λ 6 σ E μ + η + γ I + λ 7 η 1 b I μ + d + γ I H + λ 8 η b I + d I H μ R ,
where, λ m for m = 1 , 2 , 3 , , 8 represents the adjoint variables.

Solution of Optimal Control Problem

In this section, the solution to the optimal control COVID-19 as outlined in (8) is established. To achieve this, famous Pontryagin’s principle [38,39] is employed. The desired optimal solution is denoted by u ¯ 1 , u ¯ 2 , u ¯ 3 , u ¯ 4 , u ¯ 5 , u ¯ 6 , u ¯ 7 . Furthermore, the following are the corresponding necessary optimality conditions used in the solution procedure stated as
d z d t = λ m H t , u ¯ i , λ m , u H t , u ¯ i , λ m = 0 , d λ m t d t = z λ m t , u ¯ i , λ m .
The criterias mentioned in (13) and the following theorem has been utilized to derive the solution of the optimal system.
Theorem 5. 
The controls u ¯ 1 , u ¯ 2 , u ¯ 3 , u ¯ 4 , u ¯ 5 , u ¯ 6 , u ¯ 7 and the solution S ¯ , E ¯ , V ¯ 1 , V ¯ 2 , V ¯ 3 , I ¯ , I H ¯ , R ¯ of the control system (8) minimizing the objective functional in the problem, then there exist adjoint variables (co-state variables) λ m , m = 1 , 2 , , 8 . Further, the transversality conditions λ m T = 0 , m = 1 , 2 , 3 8 , such that
λ 1 t = μ λ 1 + u 2 λ 1 λ 3 + λ 1 λ 2 α I ¯ N ¯ S ¯ N ¯ 2 1 u 1 + λ 2 λ 3 1 u 3 α V ¯ 1 δ 1 I ¯ N ¯ 2 + λ 2 λ 4 1 u 4 α V ¯ 2 δ 2 I ¯ N ¯ 2 + λ 2 λ 5 1 u 5 α V ¯ 3 δ 3 I ¯ N ¯ 2 , λ 2 t = A 1 + μ λ 2 + σ λ 2 λ 6 + λ 2 λ 1 1 u 1 α I ¯ S ¯ N ¯ 2 + λ 2 λ 3 1 u 3 α V ¯ 1 δ 1 I ¯ N ¯ 2 + λ 2 λ 4 1 u 4 α V ¯ 2 δ 2 I ¯ N ¯ 2 + λ 2 λ 5 1 u 5 α V ¯ 3 δ 3 I ¯ N ¯ 2 , λ 3 t = A 2 + μ λ 3 + λ 3 λ 4 u 6 + λ 2 λ 1 1 u 1 α I ¯ S ¯ N ¯ 2 + λ 3 λ 2 1 u 3 α δ 1 I ¯ N ¯ V ¯ 1 N ¯ 2 + λ 2 λ 4 1 u 4 α V ¯ 2 δ 2 I ¯ N ¯ 2 + λ 2 λ 5 1 u 5 α V ¯ 3 δ 3 I ¯ N ¯ 2 , λ 4 t = A 3 + μ λ 4 + u 7 λ 4 λ 5 + λ 2 λ 1 1 u 1 α I ¯ S ¯ N ¯ 2 + λ 2 λ 3 1 u 3 α V ¯ 1 δ 1 I ¯ N ¯ 2 + λ 4 λ 2 1 u 4 α δ 2 I ¯ N ¯ V ¯ 2 N ¯ 2 + λ 2 λ 5 1 u 5 α V ¯ 3 δ 3 I ¯ N ¯ 2 , λ 5 t = A 4 + μ λ 5 + λ 2 λ 1 1 u 1 α I ¯ S ¯ N ¯ 2 + λ 2 λ 3 1 u 3 α V ¯ 1 δ 1 I ¯ N ¯ 2 + λ 2 λ 4 1 u 4 α V ¯ 2 δ 2 I ¯ N ¯ 2 + λ 5 λ 2 1 u 5 α δ 3 I ¯ N ¯ V ¯ 3 N ¯ 2 , λ 6 t = λ 1 λ 2 1 u 1 α S ¯ N ¯ I ¯ N ¯ 2 + λ 3 λ 2 1 u 3 α δ 1 V ¯ 1 N ¯ I ¯ N ¯ 2 + λ 4 λ 2 1 u 4 α δ 2 V ¯ 2 N ¯ I ¯ N ¯ 2 + λ 5 λ 2 1 u 5 α δ 3 V ¯ 3 N ¯ I ¯ N ¯ 2 + μ + γ λ 6 + η λ 6 λ 7 + b η λ 7 λ 8 , λ 7 t = d λ 7 λ 8 + μ + γ λ 7 + λ 2 λ 1 1 u 1 α I ¯ S ¯ N ¯ 2 + λ 2 λ 3 1 u 3 α V ¯ 1 δ 1 I ¯ N ¯ 2 + λ 2 λ 4 1 u 4 α V ¯ 2 δ 2 I ¯ N ¯ 2 + λ 2 λ 5 1 u 5 α V ¯ 3 δ 3 I ¯ N ¯ 2 , λ 8 t = μ λ 8 + λ 2 λ 1 1 u 1 α I ¯ S ¯ N ¯ 2 + λ 2 λ 3 1 u 3 α V ¯ 1 δ 1 I ¯ N ¯ 2 + λ 2 λ 4 1 u 4 α V ¯ 2 δ 2 I ¯ N ¯ 2 + λ 2 λ 5 1 u 5 α V ¯ 3 δ 3 I ¯ N ¯ 2 .
Furthermore, the associated optimal controls u ¯ 1 , u ¯ 2 , u ¯ 3 , u ¯ 4 , u ¯ 5 , u ¯ 6 , and u ¯ 7 are given by
u ¯ 1 = I ¯ S ¯ α λ 2 λ 1 A 5 N , u ¯ 2 = S ¯ λ 1 λ 3 A 6 , u ¯ 3 = I ¯ V ¯ 1 α δ 1 λ 2 λ 3 A 7 N , u ¯ 4 = I ¯ V ¯ 2 α δ 2 λ 2 λ 4 A 8 N , u ¯ 5 = I ¯ V ¯ 3 α δ 3 λ 2 λ 5 A 9 N , u ¯ 6 = V ¯ 1 λ 3 λ 4 A 10 , u ¯ 7 = V ¯ 2 λ 4 λ 5 A 11 .
Proof. 
By using the condition stated in (13), the transversality conditions and results given in (14) are obtained for the Hamiltonian function given in (12) using S = S ¯ , E = E ¯ , V 1 = V ¯ 1 , V 2 = V ¯ 2 , V 3 = V ¯ 3 , I = I ¯ , I H = I ¯ H , and R = R ¯ . Moreover, using the condition H ( t , u ¯ i , λ m ) u i = 0 given in (13), the optimal controls u ¯ 1 , u ¯ 2 , u ¯ 3 , u ¯ 4 , u ¯ 5 , u ¯ 6 , and u ¯ 7 shown in (15) are derived. □

8. Estimating the Optimal Solution

This section aims to analyze the significant impact of the proposed time-varying controls on disease incidence and potential mitigation. For this purpose, the COVID-19 control model is simulated with the aforementioned time-varying control variables to display their impact on the disease dynamics. An effective iterative approach known as RK4 is used to carry out the simulation process. The weights and balancing constants are chosen as A 1 = 10 , A 2 = 0.1 , A 3 = 0.1 , A 4 = 0.01 , A 5 = 50 , A 6 = 10 , A 7 = 30 , A 8 = 100 , A 9 = 50 , A 10 = 20 , and A 11 = 100 , while the simulation’s parameters are taken from Table 1. It is important to note that the numerical values of weighed and balanced constants are taken for the sake of simulation. The red dashed curves illustrate the dynamical behavior of various populations under varied control measures (implementing all the control measures at the time), while the black solid curves exhibit the changing behavior with constant controls. The dynamics of susceptible, primary vaccinated, exposed, and secondary vaccinated individuals with and without time-varying controls are analyzed in Figure 9, while a similar analysis of individuals with a booster shots, hospitalized, infected, and recovered individuals are presented in Figure 10. Figure 11 and Figure 12 demonstrate the corresponding control profiles. When the suggested optimal controls are actively utilized, the population in the susceptible class reduces while the population in all vaccinated classes increases significantly and reaches its maximum level with time. However, the population in the exposed and infected classes dramatically decreased and vanished after 60 and 90 days, respectively, with the implementation of time-varying controls. Because of averting infection, the populations in the hospitalized and recovered classes also reduced significantly, as shown in Figure 10c,d. The time-varying personal protection control u 1 is utilized at the maximum level from the initial time to approximately 220 days and decreases gradually until the end of the considered time level, as shown in Figure 11a. The time-varying control for first-time vaccine efficacy enhancement u 2 is initially at the maximum level until the first 100 days and then immediately reduces until the end, as can be seen in Figure 11b. Figure 11c,d and Figure 12a illustrate the intensity of controls u 3 , u 4 , and u 5 evaluating the reduction in the possibility that first, second, and third vaccine doses do not develop immunity to the original viruses, respectively. It can be observed that all of these controls are initially maintained at the maximum level and then gradually reduced until the end of the considered time interval. The implementation level of controls used for the enhancement of second- and third-time vaccinations u 6 and u 7 are shown in Figure 12b,c, respectively. These controls are implemented at the maximum level from the start to the end of the time level under consideration.

9. Conclusions

Vaccination remains one of the most effective preventive interventions and is a critical component in mitigating disease outbreaks. In this study, we developed a novel mathematical model to assess the impact of administering multiple constant and time-varying vaccines, including the first and second doses and booster shots, on infection incidence and persistence. Initially, the necessary mathematical assessment of the model is presented and basic reproduction is calculated. To curb infection, we determined the critical vaccination convergence rate, which is dependent on the reproduction number. Moreover, using backward bifurcation analysis, we concluded that bifurcation depends on the sign of a ¯ . To measure each parameter’s relative influence on disease spread, the sensitivity indices of the proposed model parameters for the basic reproduction number are tabulated using a normalized approach. Furthermore, we reconstructed the model using optimal control theory to identify the best control strategy for minimizing infection. We introduced seven time-varying controls, where u 1 ( t ) is the personal protection control used for the reduction in effective contacts of infected individuals with susceptible, u 2 ( t ) , u 6 ( t ) , u 7 ( t ) controls are used for the 1st, 2nd, and 3rd vaccine dose efficacy enhancement. The controls u 3 ( t ) , u 4 ( t ) , u 5 ( t ) are used to reduce the possibility that after 28 days of vaccination, those who received the first, second, and third doses do not develop immunity to the original virus. Using the well-known Pontryagin’s maximum principle, the necessary optimal conditions are determined. The implementation of the suggested time-varying optimal controls has been found to play a crucial role in effectively minimizing the risk of disease transmission. Moreover, by effectively adjusting vaccination strategies, such as the timing and frequency of doses, the model shows that the spread of the disease can be significantly controlled. These time-varying controls allow for a more adaptive and responsive approach, enabling healthcare authorities to tailor vaccination campaigns on the basis of evolving epidemiological conditions, the emergence of new variants, and the overall vaccination coverage in the population. This flexibility empowers the public health system to optimize vaccination efforts and better control the spread of the infection, ultimately leading to a reduced incidence of the disease and improved public health outcomes.

Author Contributions

Y.L.: Conceptualization, formal analysis, writing—review and editing, S. and L.Z.: Software, methodology, visualization, validation, writing—review & editing, E.A.A.I. and F.A.A.: Conceptualization, supervision, project administration, formal analysis, writing—review and editing, funding acquisition, A.M.H.: Supervision, investigation, funding acquisition, resources. All authors have read and agreed to the published version of the manuscript.

Funding

This project is funded by King Saud University, Riyadh, Saudi Arabia.

Data Availability Statement

Not applicable.

Acknowledgments

Researchers Supporting Project number (RSPD2023R1060), King Saud University, Riyadh, Saudi Arabia.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Andersen, K.G.; Rambaut, A.; Lipkin, W.I.; Holmes, E.C.; Garry, R.F. The proximal origin of SARS-CoV-2. Nat. Med. 2020, 26, 450–452. [Google Scholar] [CrossRef] [PubMed]
  2. World Health Organization. Novel Coronavirus (2019-nCoV): Situation Report; World Health Organization: Geneva, Switzerland, 2020; p. 11. [Google Scholar]
  3. Pfefferbaum, B.; North, C.S. Mental health and the COVID-19 pandemic. N. Engl. J. Med. 2020, 383, 510–512. [Google Scholar] [CrossRef] [PubMed]
  4. Tu, Y.F.; Chien, C.S.; Yarmishyn, A.A.; Lin, Y.Y.; Luo, Y.H.; Lin, Y.T.; Lai, W.Y.; Yang, D.M.; Chou, S.J.; Yang, Y.P.; et al. A review of SARS-CoV-2 and the ongoing clinical trials. Int. J. Mol. Sci. 2020, 21, 2657. [Google Scholar] [CrossRef]
  5. Guan, W.J.; Ni, Z.Y.; Hu, Y.; Liang, W.H.; Ou, C.Q.; He, J.X.; Liu, L.; Shan, H.; Lei, C.L.; Hui, D.S.; et al. Clinical characteristics of coronavirus disease 2019 in China. N. Engl. J. Med. 2020, 382, 1708–1720. [Google Scholar] [CrossRef]
  6. Lauer, S.A.; Grantz, K.H.; Bi, Q.; Jones, F.K.; Zheng, Q.; Meredith, H.R.; Azman, A.S.; Reich, N.G.; Lessler, J. The incubation period of coronavirus disease 2019 (COVID-19) from publicly reported confirmed cases: Estimation and application. Ann. Intern. Med. 2020, 172, 577–582. [Google Scholar] [CrossRef]
  7. World Health Organization. Non-Pharmaceutical Public Health Measures for Mitigating the Risk and Impact of Epidemic and Pandemic Influenza: Annex: Report of Systematic Literature Reviews (No. WHO/WHE/IHM/GIP/2019.1); World Health Organization: Geneva, Switzerland, 2019. [Google Scholar]
  8. Polack, F.P. Safety and Efficacy of the BNT162b2 mRNA COVID-19 Vaccine. N. Engl. J. Med. 2020, 383, 2603–2615. [Google Scholar] [CrossRef] [PubMed]
  9. Al-arydah, M. Mathematical modeling and optimal control for COVID-19 with population behavior. Math. Meth. Appl. Sci. 2023, 1–15. [Google Scholar] [CrossRef]
  10. Al-arydah, M.; Berhe, H.; Dib, K.; Madhu, K. Mathematical modeling of the spread of the coronavirus under strict social restrictions. Math. Meth. Appl. Sci. 2021, 1–11. [Google Scholar] [CrossRef]
  11. Aatif, A.; Ullah, S.; Khan, M.A. The impact of vaccination on the modeling of COVID-19 dynamics: A fractional order model. Nonlinear Dyn. 2022, 110, 3921–3940. [Google Scholar]
  12. Zarin, R. Numerical study of a nonlinear COVID-19 pandemic model by finite difference and meshless methods. Partial. Differ. Equations Appl. Math. 2022, 6, 100460. [Google Scholar] [CrossRef]
  13. Alshehri, A.; Ullah, S. A numerical study of COVID-19 epidemic model with vaccination and diffusion. Math. Biosci. Eng. 2023, 20, 4643–4672. [Google Scholar] [CrossRef]
  14. Ali, I.; Khan, S.U. Dynamics and simulations of stochastic COVID-19 epidemic model using Legendre spectral collocation method. AIMS Math. 2023, 8, 4220–4236. [Google Scholar] [CrossRef]
  15. Din, A.; Amine, S.; Allali, A. A stochastically perturbed co-infection epidemic model for COVID-19 and hepatitis B virus. Nonlinear Dyn. 2023, 111, 1921–1945. [Google Scholar] [CrossRef]
  16. Rahat, Z.; Khan, A.; Yusuf, A.; Sayed, A.-K.; Inc, M. Analysis of fractional COVID-19 epidemic model under Caputo operator. Math. Methods Appl. Sci. 2023, 46, 7944–7964. [Google Scholar]
  17. Ravichandran, C.; Logeswari, K.; Khan, A.; Abdeljawad, T.; Goamez-Aguilar, J.F. An epidemiological model for computer virus with Atangana-Baleanu fractional derivative. Results Phys. 2023, 51, 106601. [Google Scholar] [CrossRef]
  18. Lou, J.; Zheng, H.; Zhao, S.; Cao, L.; Wong, E.L.; Chen, Z.; Chan, R.W.; Chong, M.K.; Zee, B.C.; Chan, P.K.; et al. Quantifying the effect of government interventions and virus mutations on transmission advantage during COVID-19 pandemic. J. Infect. Public Health 2022, 15, 338–342. [Google Scholar] [CrossRef]
  19. Wang, Y.; Wang, P.; Zhang, S.; Pan, H. Uncertainty modeling of a modified SEIR epidemic model for COVID-19. Biology 2022, 11, 1157. [Google Scholar] [CrossRef]
  20. Liu, P.; Huang, X.; Zarin, R.; Cui, T.; Din, A. Modeling and numerical analysis of a fractional order model for dual variants of SARS-CoV-2. Alex. Eng. J. 2023, 65, 427–442. [Google Scholar] [CrossRef]
  21. Eikenberry, S.E.; Mancuso, M.; Iboi, E.; Phan, T.; Eikenberry, K.; Kuang, Y.; Kostelich, E.; Gumel, A.B. To mask or not to mask: Modeling the potential for face mask use by the general public to curtail the COVID-19 pandemic. Infect. Dis. Model. 2020, 5, 293–308. [Google Scholar] [CrossRef] [PubMed]
  22. Agusto, F.B.; Erovenko, I.V.; Fulk, A.; Abu-Saymeh, Q.; Romero-Alvarez, D.; Ponce, J.; Sindi, S.; Ortega, O.; Saint Onge, J.M.; Peterson, A.T. To isolate or not to isolate: The impact of changing behavior on COVID-19 transmission. BMC Public Health 2022, 22, 138. [Google Scholar]
  23. Watson, O.J.; Barnsley, G.; Toor, J.; Hogan, A.B.; Winskill, P.; Ghani, A.C. Global impact of the first year of COVID-19 vaccination: A mathematical modelling study. Lancet Infect. Dis. 2022, 22, 1293–1302. [Google Scholar] [CrossRef] [PubMed]
  24. Eyre, D.W.; Taylor, D.; Purver, M.; Chapman, D.; Fowler, T.; Pouwels, K.B.; Walker, A.S.; Peto, T.E. Effect of COVID-19 vaccination on transmission of alpha and delta variants. N. Engl. J. Med. 2022, 386, 744–756. [Google Scholar] [CrossRef] [PubMed]
  25. Ngonghala, C.N.; Taboe, H.B.; Safdar, S.; Gumel, A.B. Unraveling the dynamics of the Omicron and Delta variants of the 2019 coronavirus in the presence of vaccination, mask usage, and antiviral treatment. Appl. Math. Model. 2023, 114, 447–465. [Google Scholar] [CrossRef] [PubMed]
  26. Peter, O.J.; Panigoro, H.S.; Abidemi, A.; Ojo, M.M.; Oguntolu, F.A. Mathematical model of COVID-19 pandemic with double dose vaccination. Acta Biotheor. 2023, 71, 9. [Google Scholar] [CrossRef] [PubMed]
  27. Wang, Y.; Ullah, S.S.; Khan, I.U.; AlQahtani, S.A.; Hassan, A.M. Numerical assessment of multiple vaccinations to mitigate the transmission of COVID-19 via a new epidemiological modeling approach. Results Phys. 2023, 52, 106889. [Google Scholar] [CrossRef]
  28. Lakshmikantham, V.; Leela, S.; Martynyuk, A.A. Stability Analysis of Nonlinear Systems; M. Dekker: New York, NY, USA, 1989; pp. 249–275. [Google Scholar]
  29. LaSalle, J.P.; Lefschetz, S. The Stability of Dynamical Systems (SIAM, Philadelphia, 1976); Zhonghuai Wu Yueyang Vocational Technical College Yueyang: Hunan, China, 1976. [Google Scholar]
  30. Van den Driessche, P.; Watmough, J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Math. Biosci. 2002, 180, 29–48. [Google Scholar] [CrossRef] [PubMed]
  31. Castillo-Chavez, C.; Song, B. Dynamical models of tuberculosis and their applications. Math. Biosci. Eng. 2004, 1, 361–404. [Google Scholar] [CrossRef] [PubMed]
  32. Akinwande, N.I.; Ashezua, T.T.; Gweryina, R.I.; Somma, S.A.; Oguntolu, F.A.; Usman, A.; Abdurrahman, O.N.; Kaduna, F.S.; Adajime, T.P.; Kuta, F.A.; et al. Mathematical model of COVID-19 transmission dynamics incorporating booster vaccine program and environmental contamination. Heliyon 2022, 8, e11513. [Google Scholar] [CrossRef]
  33. Kim, Y.R.; Choi, Y.J.; Min, Y. A model of COVID-19 pandemic with vaccines and mutant viruses. PLoS ONE 2022, 17, e0275851. [Google Scholar] [CrossRef]
  34. Chitnis, N.; Hyman, J.M.; Cushing, J.M. Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model. Bull. Math. Biol. 2008, 70, 1272–1296. [Google Scholar] [CrossRef]
  35. Saif, U.; Khan, M.A. Modeling the impact of non-pharmaceutical interventions on the dynamics of novel coronavirus with optimal control analysis with a case study. Chaos Solitons Fractals 2020, 139, 110075. [Google Scholar]
  36. Agusto, F.B.; Khan, M.A. Optimal control strategies for dengue transmission in Pakistan. Math. Biosci. 2018, 305, 102–121. [Google Scholar] [CrossRef]
  37. Saif, U.; Khan, M.A.; Gmez-Aguilar, J.F. Mathematical formulation of hepatitis B virus with optimal control analysis. Optim. Control. Appl. Methods 2019, 40, 529–544. [Google Scholar]
  38. Pontryagin, L.S. Mathematical Theory of Optimal Processes; CRC Press: Boca Raton, FL, USA, 1987. [Google Scholar]
  39. Fleming, W.H.; Rishel, R.W. Deterministic and Stochastic Optimal Control; Springer Science and Business Media: Berlin/Heidelberg, Germany, 2012; Volume 1. [Google Scholar]
Figure 1. Flow-Chart of the COVID-19 model with vaccinations (2).
Figure 1. Flow-Chart of the COVID-19 model with vaccinations (2).
Mathematics 11 04253 g001
Figure 2. Graphical results of the proposed model for exposed, infected, hospitalized and recovered compartments when R 0 < 1 .
Figure 2. Graphical results of the proposed model for exposed, infected, hospitalized and recovered compartments when R 0 < 1 .
Mathematics 11 04253 g002
Figure 3. Simulation of the proposed model for susceptible and vaccinated compartments using R 0 < 1 .
Figure 3. Simulation of the proposed model for susceptible and vaccinated compartments using R 0 < 1 .
Mathematics 11 04253 g003
Figure 4. Simulation of the proposed model for individuals with booster short, exposed, infected, hospitalized and recovered compartments when R 0 > 1 .
Figure 4. Simulation of the proposed model for individuals with booster short, exposed, infected, hospitalized and recovered compartments when R 0 > 1 .
Mathematics 11 04253 g004
Figure 5. Simulation of the proposed model for susceptible and primary and secondary vaccinated compartments when R 0 > 1 .
Figure 5. Simulation of the proposed model for susceptible and primary and secondary vaccinated compartments when R 0 > 1 .
Mathematics 11 04253 g005
Figure 6. Simulation of the proposed COVID-19 model with different initial values of the state variables and R 0 > 1 , where (a) susceptible, (b) exposed, (c) recovered, (d) hospitalized, (e) infected, (f) shows the individuals with booster short. The curves with different colors show the dynamics by choosing different values of the corresponding state variables.
Figure 6. Simulation of the proposed COVID-19 model with different initial values of the state variables and R 0 > 1 , where (a) susceptible, (b) exposed, (c) recovered, (d) hospitalized, (e) infected, (f) shows the individuals with booster short. The curves with different colors show the dynamics by choosing different values of the corresponding state variables.
Mathematics 11 04253 g006
Figure 7. Dynamics of the COVID-19 model with different initial values of the state variables and R 0 < 1 , where (a) secondary vaccinated, (b) individuals with booster short, (c) exposed, (d) infected, (e) hospitalized, (f) R] recovered individuals. The curves with different colors show the dynamics by choosing different values of the corresponding state variables.
Figure 7. Dynamics of the COVID-19 model with different initial values of the state variables and R 0 < 1 , where (a) secondary vaccinated, (b) individuals with booster short, (c) exposed, (d) infected, (e) hospitalized, (f) R] recovered individuals. The curves with different colors show the dynamics by choosing different values of the corresponding state variables.
Mathematics 11 04253 g007
Figure 8. Bar graph of sensitivity indices of the COVID-19 model.
Figure 8. Bar graph of sensitivity indices of the COVID-19 model.
Mathematics 11 04253 g008
Figure 9. Simulation of (a) susceptible, (b) exposed, (c) primary and (d) secondary vaccinated individuals in the model (8) with optimal and with constant controls. The constant values of first, second and booster COVID-19 vaccine controls are considered as 0.7, 0.6, and 0.29, respectively.
Figure 9. Simulation of (a) susceptible, (b) exposed, (c) primary and (d) secondary vaccinated individuals in the model (8) with optimal and with constant controls. The constant values of first, second and booster COVID-19 vaccine controls are considered as 0.7, 0.6, and 0.29, respectively.
Mathematics 11 04253 g009
Figure 10. Simulation of people with (a) booster short, (b) infected, (c) hospitalized and (d) recovered compartment in the model (8) with optimal and with constant controls. The constant values of first, second and booster COVID-19 vaccine controls are considered as 0.7, 0.6, and 0.29, respectively.
Figure 10. Simulation of people with (a) booster short, (b) infected, (c) hospitalized and (d) recovered compartment in the model (8) with optimal and with constant controls. The constant values of first, second and booster COVID-19 vaccine controls are considered as 0.7, 0.6, and 0.29, respectively.
Mathematics 11 04253 g010
Figure 11. The corresponding optimal control profiles with controls (a) u 1 ( t ) , (b) u 2 ( t ) , (c) u 3 ( t ) and (d) u 4 ( t ) .
Figure 11. The corresponding optimal control profiles with controls (a) u 1 ( t ) , (b) u 2 ( t ) , (c) u 3 ( t ) and (d) u 4 ( t ) .
Mathematics 11 04253 g011
Figure 12. The corresponding optimal control profiles with controls (a) u 5 ( t ) , (b) u 6 ( t ) and (c) u 7 ( t ) .
Figure 12. The corresponding optimal control profiles with controls (a) u 5 ( t ) , (b) u 6 ( t ) and (c) u 7 ( t ) .
Mathematics 11 04253 g012
Table 1. Description with numerical values of the model’s parameters.
Table 1. Description with numerical values of the model’s parameters.
ParameterMeaningValueReference
θ humans’ recruitment rate7,828,143[32]
μ natural death rate0.011380[32]
ζ 1 rate of first COVID vaccine dose 0.710 [32]
ζ 2 rate of second COVID vaccine dose0.650Assumed
ζ 3 rate of booster shot 0.290 [32]
δ 1 possibility that after 28 days of vaccinations, those who received the first vaccine dose has not produced immunity to the original virus 4.7 % [33]
δ 2 possibility that after 28 days of vaccinations, those who received the 2nd vaccine dose has not produced immunity to the original virus 0.2 % [33]
δ 3 possibility that after 28 days of vaccinations, those who received the booster shot has not produced immunity to the original virus 0.1 % [33]
σ transmission rate of exposed to infectious class 1 / 5.2 [26]
brecovery rate of infectious people 0.76210 [26]
η flow rate of infectious individuals 0.0720 [32]
α effective contacts rate 0.750 Assumed
γ mortality rate due to infection 0.000010 Assumed
drecovery rate of hospitalization individuals 0.01070 Assumed
Table 2. Sensitivity index of the selected model’s parameters.
Table 2. Sensitivity index of the selected model’s parameters.
ParameterIndex
Y α 1
Y σ 0.0558698
Y μ 0.73501
Y δ 1 0.0452883
Y δ 2 0.0041564
Y δ 3 0.0529594
Y γ −0.000119918
Y η 0.863413
Y ζ 1 −0.881821
Y ζ 2 −0.0435263
Y ζ 3 −0.00199973
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Li, Y.; Samreen; Zada, L.; Ismail, E.A.A.; Awwad, F.A.; Hassan, A.M. Assessing the Impact of Time-Varying Optimal Vaccination and Non-Pharmaceutical Interventions on the Dynamics and Control of COVID-19: A Computational Epidemic Modeling Approach. Mathematics 2023, 11, 4253. https://doi.org/10.3390/math11204253

AMA Style

Li Y, Samreen, Zada L, Ismail EAA, Awwad FA, Hassan AM. Assessing the Impact of Time-Varying Optimal Vaccination and Non-Pharmaceutical Interventions on the Dynamics and Control of COVID-19: A Computational Epidemic Modeling Approach. Mathematics. 2023; 11(20):4253. https://doi.org/10.3390/math11204253

Chicago/Turabian Style

Li, Yan, Samreen, Laique Zada, Emad A. A. Ismail, Fuad A. Awwad, and Ahmed M. Hassan. 2023. "Assessing the Impact of Time-Varying Optimal Vaccination and Non-Pharmaceutical Interventions on the Dynamics and Control of COVID-19: A Computational Epidemic Modeling Approach" Mathematics 11, no. 20: 4253. https://doi.org/10.3390/math11204253

APA Style

Li, Y., Samreen, Zada, L., Ismail, E. A. A., Awwad, F. A., & Hassan, A. M. (2023). Assessing the Impact of Time-Varying Optimal Vaccination and Non-Pharmaceutical Interventions on the Dynamics and Control of COVID-19: A Computational Epidemic Modeling Approach. Mathematics, 11(20), 4253. https://doi.org/10.3390/math11204253

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop