Mathematical modelling of COVID-19 dynamics using SVEAIQHR model

Ambalarajan Venkatesh; Mallela Ankamma Rao; Murugadoss Prakash Raj; Karuppusamy Arun Kumar; D. K. K. Vamsi

doi:10.1515/cmb-2023-0112

Article Open Access

Mathematical modelling of COVID-19 dynamics using SVEAIQHR model

Ambalarajan Venkatesh , Mallela Ankamma Rao , Murugadoss Prakash Raj , Karuppusamy Arun Kumar and D. K. K. Vamsi

Published/Copyright: February 29, 2024

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information

From the journal Computational and Mathematical Biophysics Volume 12 Issue 1

Abstract

In this study, we formulate an eight-compartment mathematical model with vaccination as one of the compartments to analyze the dynamics of COVID-19 transmission. We examine the model’s qualitative properties, such as positivity and boundedness of solutions, and stability analysis of the illness-free equilibrium with respect to the basic reproduction number. We estimate ten significant parameters and also compute the magnitude of the basic reproduction number for India by fitting the proposed model to daily confirmed and cumulative confirmed COVID-19 cases in India. Sensitivity analysis with respect to basic reproduction number is conducted, and the main parameters that impact the widespread of disease are determined. We further extend this model to an optimal control problem by including four non-pharmaceutical and pharmaceutical intervention measures as control functions. Our numerical results show that the four control strategy has greater impact than the three control strategies, two control strategies, and single control strategies on reducing the dynamics of COVID-19 transmission.

Keywords: COVID-19; epidemic model; vaccination; bifurcation; sensitivity analysis; optimal control problem

MSC 2010: 92D30; 34C23; 49J15

1 Introduction

Currently, most countries are affected by COVID-19 pandemic and are experiencing enormous cases of infection and fatalities. Up to August 16, 2023, globally 769,806,130 infected cases and 6,955,497 death cases were reported from this COVID-19 pandemic [40]. India was facing COVID-19 pandemic with massive infected and death cases. As of August 16, 2023, India has reported 44,996,599 COVID-19 cases, with 531,925 illness fatalities, determining India as the third country among those countries affected by the pandemic.

Due to different dangerous variants of COVID-19 pandemic, the researchers, biologists, and medical professionals are constantly working to develop powerful vaccines and remedies for minimization of COVID-19 infection. The available literature demonstrates that numerous research articles with various viewpoints have been written and published in relation to minimizing COVID-19 infection. Mathematical models are the best possible way to determine the transmission dynamics of infectious and evaluate the effective of intervention measures in reducing the prevalence of diseases. Numerous mathematical models with vaccination [2,4,14,25,34] to reduce the COVID-19 disease burden have already been performed. To calculate the magnitude of the COVID-19 epidemic in Wuhan, Imai et al. [13] used computational modeling with an attention on the human-to-human transmission. Their results reveals that the preventative measures must effectively stop well over 60% of disease transmission. In order to lessen the impact of COVID-19, Ngonghala et al. [26] emphasized in his study the importance of maintaining social distance. Ndaïrou et al. [24] presented a deterministic model for the COVID-19 disease transmission, emphasizing on the spreadability of super-spreaders among people. To reduce the spread of COVID-19, Mandal et al. [21] developed a mathematical model that incorporates a quarantine compartment and government intervention measures. To study the effects of confined, hospitalized and latent compartment persons, Prathumwan et al. [30] produced a mathematical model for assessing the dynamics of COVID-19 transmission. Okuonghae and Omame [27] investigated the impact of control measures including face masks, social isolation, and detecting cases in reducing the dynamics of COVID-19. The role of face mask use, hospitalization, and quarantine in minimizing the COVID-19 pandemic in India were studied by Srivastav et al. [35]. Khajanchi et al. [16] proposed a model that categorizes infection into nine stages, and it was found that combination pharmaceutical and non-pharmaceutical preventive measures is more effective than individual measure. Rai et al. [31] constructed a mathematical model to evaluate the effectiveness of social media advertising in battling the COVID-19 outbreak in India using non-pharmaceutical intervention (NPI) measures. Acuña-Zegarra et al. [1] created a mixed constraint optimal control problem to define vaccination schedules. The outcomes imply that the main factor in COVID-19 mitigation would be the response to immunity is caused by vaccines and reinfection periods. Kurmi and Chouhan [17] created a multi-compartment model using vaccination as a control to evaluate the effect of immunization on the spread of COVID-19. Based on studies [10], the combination of therapy and vaccine measures is more effective in minimizing the COVID spread and improving the recovery rate of infected individuals. Ankamma Rao and Venkatesh [3] formulated the 9-compartment mathematical model including quarantine, hospitalization and insusceptible by using NPIs to mitigate the spread of the COVID disease. By employing the pharmaceutical interventions, Venkatesh et al. [37] developed a multi-strain mathematical model with 16 compartments to reduce the COVID-19 infection and improve recovery. Venkatesh and Ankamma Rao [38] constructed an epidemic model that comprises effective and ineffective vaccination to analyze the dissemination of COVID-19. Furthermore, an optimal control model was developed to evaluate the vaccination and treatment measures in mitigating the spread of COVID-19. To assess the effectiveness of vaccine and treatment approaches in reducing the COVID-19 burden, Chhetri et al. [7] proposed a mathematical model that includes age-related COVID-19 transmission dynamics.

Inspired by all of the above work, we developed a novel deterministic compartmental model that includes the hospitalized, vaccinated, and quarantined classes. The main objective of this study is to analyze the impacts of non-pharmaceutical and pharmaceutical intervention techniques, which comprises effort of awareness initiatives, vaccination for susceptible individuals, and treatment for infected individuals on the controlling of COVID-19 spread.

The rest of the article is organized as follows: Section 2 contains the model formulation; in Section 3, the positivity and boundedness of solutions are verified and stability analysis of illness- free equilibrium (IFE) with respect to R 0 is performed; in Section 4, the model is fitted for COVID-19 data in India and estimated parameter values, sensitive analysis of R 0 is performed to identify the significant of parameters, and the effect of vaccination rate on infected populations is verified; furthermore, in Section 5, an optimal control problem is developed using four control variables and numerical simulation for combinations of a single, two, and three control strategies are provided; and finally, a thorough discussion of our results is presented in Section 6.

2 Model formulation

In this study , we formulate an epidemiological model with eight compartments for the transmission dynamics COVID-19 disease in India. The total human population N ( t ) has been divided into eight distinct compartments, namely susceptible S ( t ) , vaccinated V ( t ) , exposed E ( t ) , asymptomatic infected A ( t ) , symptomatic infected I ( t ) , quarantined Q ( t ) , hospitalized H ( t ) , and recovered R ( t ) . Then,

N ( t ) = S ( t ) + V ( t ) + E ( t ) + I ( t ) + A ( t ) + H ( t ) + R ( t ) .

In order to structure the system of differential equation and obtain the desired model, we assume the following assumptions.

Let Λ be birth or recruitment rate of individuals in the susceptible population, and let δ be the natural death rate. Suppose that ζ and ρ are adjustment factors for asymptomatic and symptomatic infected populations. Assume that the disease transmission rate is β and vaccination rate of susceptible individuals is α . Some people who received vaccinations were re-exposed to the pathogen as a result of vaccine inefficiency or vaccine failure. Let η be transfer rate of vaccinated individuals to exposed individuals. A portion σ of exposed individuals becomes asymptomatic infected, and the remaining portion ( 1 − σ ) becomes symptomatic infected at the rate ψ . The symptomatic individuals with mild symptoms are quarantined/isolated at the rate π i , and those with extremely severe symptoms are hospitalized at the rate κ i . Due to their severe illness, some symptomatic infected people die at a rate of δ i . Similarly, the asymptomatic individuals are quarantined at the rate π a and some of the asymptomatic infected individuals may recover at a rate γ a , while others might die at the rate δ a before showing symptoms. It may happen for individuals under quarantine to experience specific difficulties and be hospitalized at a rate of κ q . Some of the hospitalized individuals recover at the rate γ h , and other critically ill individuals have a chance of not recovering even after being hospitalized, so their mortality rate is δ h .

A schematic representation of above assumptions is shown in Figure 1, and the nonlinear differential equations of the proposed model is given as follows:

(1) d S d t = Λ − β ( ζ A + ρ I ) S N − ( α + δ ) S d V d t = α S − β η ( ζ A + ρ I ) V N − δ V d E d t = β ( ζ A + ρ I ) S N + β η ( ζ A + ρ I ) V N − ( ψ + δ ) E d A d t = σ ψ E − ( π a + γ a + δ a + δ ) A d I d t = ( 1 − σ ) ψ E − ( π i + κ i + δ i + δ ) I d Q d t = π a A + π i I − ( κ q + γ q + δ q + δ ) Q d H d t = κ i I + κ q Q − ( γ h + δ h + δ ) H d R d t = γ a A + γ q Q + γ h H − δ R ,

with initial conditions as follows:

(2) S ( 0 ) ≥ 0 , V ( 0 ) ≥ 0 , E ( 0 ) ≥ 0 , A ( 0 ) ≥ 0 , I ( 0 ) ≥ 0 , Q ( 0 ) ≥ 0 , H ( 0 ) ≥ 0 , and D ( 0 ) ≥ 0 .

Table 1 presents a description of all the parameters that are used in system (1).

Figure 1

Schematic diagram of S V E A I Q H R model.

Table 1

The explanation of parameters of S V E A I Q H R model

Parameter	Description
Λ	Recruitment rate of individuals in the susceptible class
( 1 ‒ σ )	Fraction of exposed individuals move to symptomatic infected class
ψ	Incubation rate of exposed individuals
α	Vaccination rate of susceptible individuals
η	Transition rate for vaccinated individuals become exposed
ζ	Modifying factor for asymptomatic infected individuals
ρ	Modifying factor for symptomatic infected individuals
β	Transmission rate of virus
π a	Quarantine rate of asymptomatic infected individuals
π i	Quarantine rate of symptomatic infected individuals
κ q	Hospitalization rate of symptomatic infected individuals
κ i	Hospitalization rate of quarantined individuals
γ a	Recovery rate of asymptomatic infected individuals
γ q	Recovery rate of quarantined individuals
γ h	Recovery rate of hospitalized individuals
δ a	Mortality rate of asymptomatic infected individuals
δ i	Mortality rate of symptomatic infected individuals
δ q	Mortality rate of quarantined individuals
δ h	Mortality rate of hospitalized individuals
δ	Natural mortality rate of human individuals

3 S V E A I Q H R model analysis

3.1 Positivity and boundedness

Theorem 3.1

The solutions of system (1) with initial conditions (2) remain positive over the region in R + 8 .

Proof

The positivity of the solutions of system (1) is established using a method described in [7, 8].

In order to prove the positivity of system (1), we consider any solution emerging from the nonnegative region R + 8 which remains positive for every t > 0 .

To do this, we demonstrate that the vector field points on each plane is bounded by nonnegative region R + 8 .

Using the equations of system (1), we obtain

d S d t S = 0 = Λ ≥ 0 , d V d t V = 0 = α S ≥ 0 , d E d t E = 0 = β ( ζ A + ρ I ) S N + η V N ≥ 0 , d A d t A = 0 = σ ψ E ≥ 0 , d I d t I = 0 = ( 1 − σ ) ψ E ≥ 0 , d Q d t Q = 0 = π a A + π i I ≥ 0 , d H d t H = 0 = κ i I + κ q Q ≥ 0 , and d R d t R = 0 = γ a A + γ q Q + γ h H ≥ 0 .

Therefore, on the bounding planes indicated by S = 0 , V = 0 , A = 0 , I = 0 , Q = 0 , H = 0 , and R = 0 of the nonnegative region R + 8 , all of the above rates are nonnegative.

Hence, it can be concluded that all solutions of system (1) will remain positive for any time t > 0 , given that the initial conditions (2) are positive.□

Theorem 3.2

The solutions of system (1) with initial conditions (2) are uniformly bounded in the region Φ = { ( S , V , E , A , I , Q , H , R ) ∈ R + 8 : 0 ≤ S + V + E + A + I + Q + H + R ≤ Λ δ } .

Proof

Let N = S + V + E + A + I + Q + H + R

By adding all equations of system (1), we obtain

d N d t = Λ − ( δ a A + δ i I + δ q Q + δ h H ) − δ N ≤ Λ − δ N

⇒ d N d t + δ N ≤ Λ

This leads to

N ( t ) ≤ Λ δ + N ( 0 ) − Λ δ e − δ t

Considering t → ∞ , we have

lim t → ∞ sup N ( t ) ≤ Λ δ

Thus, we obtain N ( t ) ≤ Λ δ , and consequently, S + V + E + A + I + Q + H + R ≤ Λ δ .

Therefore, every solution trajectory ( S , V , E , A , I , Q , H , R ) of system (1) with the initial conditions (2) is uniformly bounded in region Φ .□

3.2 Basic reproduction number R 0

In order to attain the IFE E 0 of system (1), for E = A = I = Q = H = 0 , all equations of equation (1) are set to zero and by solving, we obtain

E 0 = ( S 0 , V 0 , E 0 , A 0 , I 0 , Q 0 , H 0 , R 0 ) = Λ ( α + δ ) , α Λ δ ( α + δ ) , 0 , 0 , 0 , 0 , 0 , 0 .

One of the most significant factor to evaluate the potential of disease spread within a population is the basic reproduction number R 0 . The basic reproduction number is the expected number of secondary infections caused by single primary infected individual over his infectious period in an otherwise susceptible population. We compute the R 0 value by using the next-generation matrix approach [36]. Now for the new infection and transition rates, we obtain the matrices ℱ and V , respectively, as follows:

ℱ = β ( ζ A + ρ I ) S N + β η ( ζ A + ρ I ) V N 0 0 and V = ( ψ + δ ) E − σ ψ E + ( π a + γ a + δ a + δ ) A − ( 1 − σ ) ψ E + ( π i + κ i + δ i + δ ) I .

At E 0 , the variational matrix F = d ℱ ∕ d Z and V = d V ∕ d Z , where Z = [ E , A , I ] will be determined as follows:

F = 0 β Λ ζ ( δ + η α ) δ ( α + δ ) β Λ ρ ( δ + η α ) δ ( α + δ ) 0 0 0 0 0 0 and V = ( ψ + δ ) 0 0 − σ ψ ( π a + γ a + δ a + δ ) 0 − ( 1 − σ ) ψ 0 ( π i + κ i + δ i + δ ) .

Since the basic reproduction number is the largest eigenvalue of the matrix F V − 1 , we have the basic reproduction number for the proposed model as follows:

R 0 = β Λ ρ ( δ + η α ) σ ψ δ ( α + δ ) ( π a + γ a + δ a + δ ) ( ψ + δ ) + β Λ ζ ( δ + η α ) ( 1 − σ ) ψ δ ( α + δ ) ( π i + κ i + δ i + δ ) ( ψ + δ )

3.3 Stability analysis of infection free equilibrium

Theorem 3.3

If R 0 < 1 , then the IFE E 0 is locally asymptotically stable; otherwise, unstable for R 0 > 1 .

Proof

The Jacobian matrix for system (1) at IFE E 0 is given as follows:

J E 0 = − ( α + δ ) 0 0 − β ζ S 0 − β ρ S 0 0 0 0 α − δ 0 − β ζ η V 0 − β ρ η V 0 0 0 0 0 0 − ( ψ + δ ) β ζ ( S 0 + η V 0 ) β ρ ( S 0 + η V 0 ) 0 0 0 0 0 σ ψ − ( π a + γ a + δ a + δ ) 0 0 0 0 0 0 ( 1 − σ ) ψ 0 − ( π i + κ i + δ i + δ ) 0 0 0 0 0 0 π a π i − ( κ q + γ q + δ q + δ ) 0 0 0 0 0 0 κ i κ q − ( γ h + δ h + δ ) 0 0 0 0 γ a 0 γ q γ h − δ .

Let us consider ( α + δ ) = K 1 , ( ψ + δ ) = K 2 , ( π a + γ a + δ a + δ ) = K 3 , ( π i + κ i + δ i + δ ) = K 4 , ( κ q + γ q + δ q + δ ) = K 5 , and ( γ h + δ h + δ ) = K 6 in the above matrix.

By substituting S 0 = Λ ( α + δ ) and V 0 = α Λ δ ( α + δ ) , the above matrix becomes

J E 0 = − K 1 0 0 − β ζ Λ α + δ − β ρ Λ α + δ 0 0 0 α − δ 0 − β ζ η α Λ δ ( α + δ ) − β ρ η α Λ δ ( α + δ ) 0 0 0 0 0 − K 2 β Λ ζ ( δ + η α ) δ ( α + δ ) β Λ ρ ( δ + η α ) δ ( α + δ ) 0 0 0 0 0 σ ψ − K 3 0 0 0 0 0 0 ( 1 − σ ) ψ 0 − K 4 0 0 0 0 0 0 π a π i − K 5 0 0 0 0 0 0 κ i κ q − K 6 0 0 0 0 γ a 0 γ q γ h − δ .

Then, the characteristic equation of the matrix J E 0 is ∣ J E 0 − λ I ∣ = 0

⇒ ( λ + δ ) ( λ + δ ) ( λ + K 1 ) ( λ + K 5 ) ( λ + K 6 ) ( λ 3 + c 1 λ 2 + c 2 λ + c 3 ) = 0 , where

c 1 = K 2 + K 3 + K 4 , c 2 = K 2 K 3 + K 3 K 4 + K 2 K 4 − β Λ ( δ + η α ) δ ( α + δ ) ( ζ σ ψ + ρ σ ( 1 − ψ ) ) = β Λ ( δ + η α ) δ ( α + δ ) ζ σ ψ K 4 K 3 + ρ σ ( 1 − ψ ) K 3 K 4 + K 3 K 4 + K 2 ( K 3 + K 4 ) ( 1 − R 0 ) c 3 = K 2 K 3 K 4 ( 1 − R 0 ) .

Hence, − δ , − δ , − K 1 , − K 5 , and − K 6 are first five eigenvalues of Jacobian matrix J E 0 , and the remaining three eigenvalues obtained from zeros of cubic equation ( λ 3 + c 1 λ 2 + c 2 λ + c 3 ) = 0 .

Liénard-Chipart test states [11,19] that the IFE E 0 is locally asymptotically stable if and only if c 1 > 0 , c 2 > 0 , c 3 > 0 , and c 1 c 2 > c 3 .

Clearly, c 1 > 0 , and if R 0 < 1 , then c 2 > 0 and c 3 > 0

c 1 c 2 − c 3 = ( K 2 + K 3 + K 4 ) β Λ ( δ + η α ) δ ( α + δ ) ζ σ ψ K 4 K 3 + ρ σ ( 1 − ψ ) K 3 K 4 + ( K 2 + K 3 + K 4 ) K 3 K 4 + K 2 [ ( K 2 + K 4 ) ( K 3 + K 4 ) + K 3 2 ] ( 1 − R 0 ) .

Therefore, c 1 c 2 > c 3 for R 0 < 1 .

Thus, the two Liénard-Chipart test conditions are satisfied for R 0 < 1 .

Hence, IFE E 0 is locally asymptotically stable for R 0 < 1 .

Since we obtain c 3 > 0 at R 0 > 1 , according to Descarte’s rule of signs, there must be at least one positive eigenvalue.

Therefore, E 0 is locally asymptotically unstable for R 0 > 1 .□

Theorem 3.4

The IFE E 0 = ( X * , 0 ) of system (1) is globally asymptotically stable for R 0 < 1 .

Proof

In system (1) the disinfection compartments are S , V , and R and infection compartments are E , A , I , Q , and H .

Then system (1) can be written as d X d t = M ( X , Y ) and d Y d t = W ( X , Y ) , W ( X , 0 ) = 0 , where X = ( S , V , R ) ∈ R + 3 and Y = ( E , A , I , Q , H ) ∈ R + 5

Using the Castillo-Chavez technique [5], the global stability of IFE at E 0 = ( X * , 0 ) for R 0 < 1 will be determined if the following two conditions are satisfied:

d X d t = M ( X , 0 ) where X * was globally asymptotically stable.
W ( X , Y ) = U Y − W ¯ ( X , Y ) , where W ¯ ( X , Y ) ≥ 0 for all ( X , Y ) ∈ Δ , U = D Y W ( X * , 0 ) is a Metzler Matrix having the positive nondiagonal entries, and Δ represents the feasible region.

For system (1), M ( X , Y ) and W ( X , Y ) are defined as follows:

M ( X , Y ) = Λ − β ( ζ A + ρ I ) S N − ( α + δ ) S α S − β η ( ζ A + ρ I ) V N − δ V γ a A + γ q Q + γ h H − δ R

and

W ( X , Y ) = β ( ζ A + ρ I ) S N − β η ( ζ A + ρ I ) V N − ( ψ + δ ) E σ ψ E − ( π a + γ a + δ a + δ ) A ( 1 − σ ) ψ E − ( π i + κ i + δ i + δ ) I π a A + π i I − ( κ q + γ q + δ q + δ ) Q κ i I + κ q Q − ( γ h + δ h + δ ) H .

From the above expressions, it can be easily shown that W ( X , 0 ) = 0 and

M ( X , 0 ) = Λ − ( α + δ ) S α S − δ V − δ R ,

and we have

d d t S V R = M ( X , 0 ) = Λ − ( α + δ ) S α S − δ V − δ R .

By solving the above system of equations analytically, we obtain

S ( t ) = Λ α + δ + S ( 0 ) α + δ e − ( α + δ ) t , V ( t ) = α Λ δ ( α + δ ) + V ( 0 ) δ ( α + δ ) e − ( α + δ ) t ,

and R ( t ) = R ( 0 ) e − δ t .

Then, as t → ∞ , S ( t ) = Λ α + δ , V ( t ) = α Λ δ ( α + δ ) , and R ( t ) = 0

Hence, X * = Λ ( α + δ ) , α Λ δ ( α + δ ) , 0 is globally asymptotically stable for d X d t = M ( X , 0 )

Hence, the first condition holds.

Now for system (1), the matrices U and W ¯ ( X , Y ) can be expressed as follows:

U = − ( ψ + δ ) β Λ ζ ( δ + η α ) δ ( α + δ ) β Λ ρ ( δ + η α ) δ ( α + δ ) 0 0 σ ψ − ( π a + γ a + δ a + δ ) 0 0 0 ( 1 − σ ) ψ 0 − ( π i + κ i + δ i + δ ) 0 0 0 π a π i − ( κ q + γ q + δ q + δ ) 0 0 0 κ i κ q − ( γ h + δ h + δ )

and

W ¯ ( X , Y ) = β ζ A 1 − S N + β ρ I 1 − S N 0 0 0 0 .

Since all nondiagonal elements of matrix U are nonnegative, U is a Metzler matrix.

Since S ( t ) ≤ N ( t ) , 1 − S N ≥ 0 .

Thus, W ¯ ( X , Y ) ≥ 0 in the region Δ as S ( t ) ≤ N ( t ) .

Therefore, the second condition is also satisfied.

Since the two conditions are satisfied, the IFE E 0 of system (1) is globally asymptotically stable for R 0 < 1 .□

3.4 Existence of endemic equilibrium

The EE Σ * = ( S * , V * , E * , A * , I * , Q * , H * , R * ) of proposed model will be determined by equating all equations of system (1) to zero and then find S * , V * , E * , A * , I * , Q * , H * and R * .

Let χ * = β ( ζ A + ρ I ) 1 N . Then,

S * = Λ ( 1 − σ ) ψ ( π a + γ a + δ a + δ ) ( β ζ σ ψ ( π i + κ i + δ i + δ ) + β ρ ( 1 − σ ) ψ ( π a + γ a + δ a + δ ) ) I * + ( 1 − σ ) ψ ( α + δ ) , V * = ( π i + κ i + δ i + δ ) ( ψ + δ ) η ( β ζ ρ ψ ( π i + κ i + δ i + δ ) + β ρ ( 1 − σ ) ψ ( π a + γ a + δ a + δ ) ) I * + Λ ( 1 − σ ) ψ ( π a + γ a + δ a + δ ) η ( β ζ ρ ψ ( π i + κ i + δ i + δ ) + β ρ ( 1 − ψ ) ψ ( π a + γ a + δ a + δ ) ) I * + ( 1 − σ ) ψ ( α + δ ) , E * = ( π i + κ i + δ i + δ ) I * ( 1 − σ ) ψ , A * = σ ψ ( π i + κ i + δ i + δ ) ( 1 − σ ) ψ ( π a + γ a + δ a + δ ) , Q * = ( σ ψ ( π i + κ i + δ i + δ ) + ( 1 − σ ) ψ ( π a + γ a + δ a + δ ) ) I * ( 1 − σ ) ψ ( π a + γ a + δ a + δ ) ( κ q + γ q + δ q + δ ) , H * = ( ( 1 − σ ) ψ ( π a + γ a + δ a + δ ) ( κ q + γ q + δ q + δ ) κ i ) I * ( 1 − σ ) ψ ( π a + γ a + δ a + δ ) ( κ q + γ q + δ q + δ ) ( γ h + δ h + δ ) , R * = γ a ( σ ψ ( π i + κ i + δ i I * + δ ) ) δ ( 1 − σ ) ψ ( π a + γ a + δ a + δ ) + γ q ( π a σ ψ ( π i + κ i + δ i + δ ) ) + I * δ ( 1 − σ ) ψ ( π a + γ a + δ a + δ ) ( κ q + γ q + δ q + δ ) + γ h ( 1 − σ ) ψ ( π a + γ a + δ a + δ ) ( κ q + γ q + δ q + δ ) + γ h σ ψ ( π i + κ i + δ i ) I * δ ( γ h + δ h + δ ) ( 1 − σ ) ψ ( π a + γ a + δ a + δ ) ( κ q + γ q + δ q + δ ) .

For our convenience, we consider

J 1 = ( π a + γ a + δ a + δ ) , J 2 = ( π i + κ i + δ i + δ ) , J 3 = σ ψ , J 4 = ( 1 − σ ) ψ , J 5 = ( α + δ ) , J 6 = ( ψ + δ ) .

Then , Y 2 ( I * ) 2 + Y 1 I * + Y 0 = 0 ,

where Y 2 = η J 2 J 6 ( β ζ J 2 J 3 + ρ J 1 J 4 ) , Y 1 = ( J 2 J 6 ) ( β ζ J 2 J 3 + ρ J 1 J 4 ) [ J 1 ( J 4 δ + J 2 ) − J 3 η ] , and Y 0 = δ J 1 J 2 J 5 J 6 ( δ + η α ) ( 1 − R 0 ) .

The system (1) has

at least one EE that exists if Y 0 < 0 iff R 0 > 1 ,
single EE occurs at Y 1 < 0 and Y 0 = 0 iff R 0 = 1 or Y 1 2 − 4 Y 0 Y 2 = 0 ,
two EEs at Y 1 < 0 and Y 0 > 0 iff R 0 < 1 or Y 1 2 − 4 Y 0 Y 2 > 0 ,
no EE that exists if Y 1 > 0 and Y 0 > 0 iff R 0 < 1 (the polynomial has only positive coefficients).

From above conditions, backward bifurcation existence will be possible by coexistence of stable IFE and EE for R 0 < 1 . We must set Y 1 2 − 4 Y 0 Y 2 = 0 in order to analyze the backward bifurcation for the immunization model and discover a crucial value R 0 * of R 0 for unique EE as follows: R 0 * = 1 − Y 1 2 4 Y 2 J 1 J 2 J 5 J 6

As a result, backward bifurcation might happen under certain condition R 0 * < R 0 < 1 .

4 Numerical simulation

4.1 Model calibration

We fit the model to daily confirmed and cumulative confirmed COVID-19 cases, which are collected from India [42] between 30 January 2020 and 24 November 2020. Ten parameters β , α , η , π a , π i , κ i , κ q , γ a , γ q , and γ h are estimated by using nonlinear least squares technique (lsqnonline function) in MATLAB. The initial population size that we used for this simulation are S ( 0 ) = 1,24,71,85,021 , V ( 0 ) = 0 , E ( 0 ) = 6,000 , A ( 0 ) = 1 , I ( 0 ) = 1 , Q ( 0 ) = 1 , H ( 0 ) = 1 , and R ( 0 ) = 0 . The estimated and fixed parameters values were presented in Table 2. Figure 2 displays the model fit with the daily confirmed and cumulative confirmed COVID-19 cases. Here, the blue dots (blue curve) represent the reported data, while the black denotes the model simulation. Using the model parameters that are listed in Table 2, the basic reproduction number R 0 is computed as 1.971. Figure 3 represents the variation in basic reproduction number with respect to parameters σ and ψ . This Figure illustrates that the reproduction number value rises as the transition from exposed people to infected individuals increases.

Table 2

Model parameter values and their sensitive indices

Parameter	Value	Source	Sensitive indices
Λ	Varies	—	1.0000
σ	0.2	[33]	0.3221
ψ	0.5	[20]	0.2376
α	0.0033	Estimated	‒ 0.7598
η	0.0042	Estimated	0.4942
ζ	0.4	[12]	0.6995
ρ	0.4	[23]	0.5404
β	0.7448	Estimated	1.0000
π a	0.0621	Estimated	‒ 0.1028
π i	0.0634	Estimated	‒ 0.4227
κ q	0.2719	Estimated	—
κ i	0.0561	Estimated	‒ 0.3740
γ a	0.0578	Estimated	‒ 0.0961
γ q	0.2719	Estimated	—
γ h	0.0658	Estimated	—
δ a	0.0001945	[15,28]	‒ 0.0217
δ i	0.0001945	[15,28]	‒ 0.0013
δ q	0.0001945	[15,28]	—
δ h	0.0001945	[15,28]	—
δ	0.000391	[41]	‒ 0.1364

Figure 2

Model fit with daily confirmed and cumulative confirmed COVID-19 cases.

$Figure 3 Basic reproduction number R 0 {R}_{0} when σ \sigma (fraction of exposed population to obtain infected) and ψ \psi (rate of transition from exposed to asymptomatic population) differ.$

Figure 3

Basic reproduction number R 0 when σ (fraction of exposed population to obtain infected) and ψ (rate of transition from exposed to asymptomatic population) differ.

4.2 Sensitivity analysis

Sensitivity analysis plays an important role in defining the importance of biological parameters in the spread of infectious diseases. It aids to identify the parameters that have a high and low impact on the threshold R 0 , respectively. Whenever these parameters are determined, multiple approaches will be applied to achieve the optimal result. Marino et al. [22] and Rodrgues et al. [32] made a thorough research on the sensitivity for the dengue epidemic. According to Chitnis et al. [9] the normalized forward sensitivity index of R 0 with respect to factor q is the ratio of the relative change in R 0 to the relative change in q and represented as Γ q R 0 = ∂ R 0 ∂ β × q R 0 .

Figure 4 demonstrates the sensitivity indices of R 0 with respect to model parameters Λ , ζ , ρ , η , σ , α , ψ , β , π a , π i , κ i , γ a , δ a , δ i , and δ . The sensitivity indices for these parameters are listed in Table 2. From this Figure 4, we identify that the parameters Λ , ζ , ρ , η , σ , ψ , and β share positive correlation with R 0 , whereas remaining parameters π a , κ i , γ a , γ q , δ a , δ i , and δ have negative correlation with R 0 . It is found that the three parameters α , π a , and κ i have strong impact in diminishing the value R 0 than the remaining parameters. That means if the vaccination rate of susceptible individuals, as well as the quarantine rate and hospitalization rate of symptomatic infected individuals, increases, then the COVID-19 disease spread reduces.

$Figure 4 The sensitivity indices of various parameters with respect to R 0 {R}_{0} .$

Figure 4

The sensitivity indices of various parameters with respect to R 0 .

Figure 5(a) displays that R 0 will raise as disease transmission rate β rises and vaccine parameter α falls. Therefore, to lower the fundamental reproduction number R 0 , it is necessary to control the spread of COVID-19 pandemic. To do this, we must take steps to increase the α value and decrease the rate β at which the disease is transmitted. In Figure 5(b), we note that R 0 can only shift from R 0 < 1 to R 0 > 1 when the value of η increases. So to control the disease transmission, we have to increase the vaccination rate rapidly.

$Figure 5 Contour graphs of R 0 {R}_{0} with respect to parameters (a) ( β \beta , α \alpha ) and (b) ( α \alpha , η \eta ).$

Figure 5

Contour graphs of R 0 with respect to parameters (a) ( β , α ) and (b) ( α , η ).

4.3 Bifurcation diagram

Through an examination of the existence of EE, we have determined the subthreshold range of bistable equilibrium in the S V E A I Q H R model system (1). As a result, within the domain R 0 * < R 0 < 1 ( R 0 * = 0.62 ), backward bifurcation can be observed with η ∈ [ 0.01 , 0.0033 ] . Figure 6(a) illustrates the coexistence of two equilibriums. The stability of equilibriums is denoted by blue line while the instability of equilibriums is denoted by red line. In contrast, the IFE exists for R 0 < R 0 * . Undoubtedly, this restriction proves to be adequate in achieving the complete elimination of infection from the population. Similarly, if vaccination population will not be infected, i.e., η = 0 then transcritical bifurcation occurs for the system (1) and is shown in Figure 6(b). Here, for R 0 > 1 , the stability of EE denoted by green line and the instability of IFE indicated by red line. Similarly, for R 0 < 1 , the stability of IFE is represented by blue line.

Figure 6

(a) Backward bifurcation and (b) Transcritical bifurcation of system (1).

4.4 Effect of COVID-19 vaccination

In this section, we verify the trajectories of infected populations for various vaccination rates. Figure 7 demonstrates that when the vaccination rate α increased, the number of individuals in asymptomatic infected, symptomatic infected, quarantined, and hospitalized classes reduces, whereas the number of individuals in vaccinated and recovered classes rises.

$Figure 7 Effect of vaccination rate α \alpha on infected populations.$

Figure 7

Effect of vaccination rate α on infected populations.

5 Optimal control studies

5.1 Optimal control problem

Optimal control strategies are crucial in preventing the spread of the COVID-19 disease. It is necessary to develop a control strategy that reduces the proportion of infected population and related costs. In this section, we develop an optimal control problem by adding four time-depdendent control variables u 1 ( t ) , u 2 ( t ) , u 3 ( t ) , and u 4 ( t ) to system (1). In order to protect the susceptible individuals from the COVID-19 infection, the first control u 1 comprises the effort of awareness campaigns, security initiatives, and social isolation measures. According to Watson et al. [39], vaccinations worldwide averted 19.8 million deaths during the first year of COVID-19 vaccination. Thus, the second control u 2 depicts the execution of continuous vaccination. Chhetri et al. [6] reveals that the better treatment policy to minimize the spread of COVID-19 infection is the combination of the immunotherapy treatment and antiviral drugs. Therefore, we designate that the treatments for symptomatic infected population and hospitalized population are in the form b 1 u 3 and b 2 u 4 , where b 1 and b 2 are treatment rates and u 3 and u 4 are control variables. The control variables exhibit a numerical range from 0 to 1, where a value of 0 indicates no efforts are placed in these controls. Furthermore, maximal effort is associated to controls whose values are equal to 1.

The following optimal control model was developed by considering all the above presumptions:

(3) d S d t = Λ − β ( 1 − u 1 ) ( ζ A + ρ I ) S N − ( α u 2 + δ ) S d V d t = α u 2 S − β η ( ζ A + ρ I ) V N − δ V d E d t = β ( 1 − u 1 ) ( ζ A + ρ I ) S N + β η ( ζ A + ρ I ) V N − ( ψ + δ ) E d A d t = σ ψ E − ( π a + γ a + δ a + δ ) A d I d t = ( 1 − σ ) ψ E − ( π i + κ i + b 1 u 3 + δ i + δ ) I d Q d t = π a A + π i I − ( κ q + γ q + δ q + δ ) Q d H d t = κ i A + κ q Q − ( γ h + b 2 u 4 + δ h + δ ) H d R d t = γ a A + γ q Q + γ h H + b 1 u 3 I + b 2 u 4 H − δ R

subject to minimize the objective functional

J ( u 1 ( t ) , u 2 ( t ) , u 3 ( t ) , u 4 ( t ) ) = ∫ 0 T ( C 1 I + C 2 H + 1 2 ( C 3 u 1 2 + C 4 u 2 2 + C 5 u 3 2 + C 6 u 4 2 ) ) d t

with the primary conditions (2). The constants C i for i = 1 , 2 , 3 , … , 6 correspond to the weight constants, which maintain the integrands’ units in balance.

We must now determine the optimal control values u 1 * , u 2 * , u 3 * , and u 4 * such that

J ( u 1 * , u 2 * , u 3 * , u 4 * ) = min u 1 , u 2 , u 3 , u 4 ∈ U J ( u 1 , u 2 , u 3 , u 4 ) ,

where

U = { u 1 ( t ) , u 2 ( t ) , u 3 ( t ) , u 4 ( t ) : measurable and 0 ≤ u 1 ( t ) , u 2 ( t ) , u 3 ( t ) , u 4 ( t ) ≤ 1 , t ∈ [ 0 , T ] } .

Here, we find out the necessary conditions for our optimal control problem using the Pontryagin’s maximum principle [18,29].

The Lagrangian function is given as follows:

ℒ ( I , H , u 1 , u 2 , u 3 , u 4 ) = C 1 I + C 2 H + 1 2 ( C 3 u 1 2 + C 4 u 2 2 + C 5 u 3 2 + C 6 u 4 2 ) .

The Hamiltonian function ℋ for system (3) is defined as follows:

ℋ = ℒ ( I , H , u 1 , u 2 , u 3 , u 4 ) + λ 1 d S d t + λ 2 d V d t + λ 3 d E d t + λ 4 d A d t + λ 5 d I d t + λ 6 d Q d t + λ 7 d H d t + λ 8 d R d t

where λ i for i = 1 , 2 , 3 , … , 8 are the adjoint variables related to the state variables in system (3).

The following system of adjoint equations can be derived by using the partial derivatives of ℋ with respect to state variables S , V , E , A , I , Q , H , and R :

(4) d λ 1 d t = − ∂ ℋ ∂ S = ( λ 1 − λ 3 ) β ( 1 − u 1 ) ( ζ A + ρ I ) 1 N + ( λ 1 − λ 2 ) ( α u 2 ) + λ 1 δ , d λ 2 d t = − ∂ ℋ ∂ V = ( λ 2 − λ 3 ) β η ( ζ A + ρ I ) 1 N + λ 2 δ , d λ 3 d t = − ∂ ℋ ∂ E = ( λ 3 − λ 5 ) ψ + ( λ 3 − λ 5 ) σ ψ + λ 3 δ , d λ 4 d t = − ∂ ℋ ∂ A = + ( λ 1 − λ 3 ) β ( 1 − u 1 ) ρ S N + ( λ 2 − λ 3 ) β η ρ V N + ( λ 4 − λ 6 ) π a + ( λ 4 − λ 8 ) γ a + λ 4 ( δ a + δ ) , d λ 5 d t = − ∂ ℋ ∂ I = − C 1 + ( λ 1 − λ 3 ) β ( 1 − u 1 ) ρ S N + ( λ 2 − λ 3 ) β η ρ V N + ( λ 5 − λ 6 ) π i + ( λ 5 − λ 7 ) κ i + ( λ 5 − λ 8 ) b 1 u 3 + λ 5 ( δ a + δ ) , d λ 6 d t = − ∂ ℋ ∂ Q = ( λ 6 − λ 7 ) κ q + ( λ 6 − λ 8 ) γ q + λ 6 ( δ q + δ ) , d λ 7 d t = − ∂ ℋ ∂ H = − C 2 + ( λ 7 − λ 8 ) ( γ h + b 2 u 4 ) + λ 7 ( δ h + δ ) , d λ 8 d t = − ∂ ℋ ∂ R = λ 8 δ .

Now we minimize the Hamilton function relating to the optimal control variables u 1 * ( t ) , u 2 * ( t ) , u 3 * ( t ) , and u 4 * ( t ) .

Theorem 5.1

The optimal controls u 1 ∗ ( t ) , u 2 ∗ ( t ) , u 3 ∗ ( t ) , and u 4 ∗ ( t ) which minimize J over the region U are given as follows:

u 1 ∗ ( t ) = min { 1 , max { 0 , u 1 ¯ } } , u 2 ∗ ( t ) = min { 1 , max { 0 , u 2 ¯ } } , u 3 ∗ ( t ) = min { 1 , max { 0 , u 3 ¯ } } a n d u 4 ∗ ( t ) = min { 1 , max { 0 , u 4 ¯ } } ,

where u 1 ¯ = ( λ 1 − λ 3 ) ( ζ A + ρ I ) S N C 3 , u 2 ¯ = ( λ 1 − λ 2 ) α S C 4 , u 3 ¯ = ( λ 5 − λ 8 ) b 1 I C 5 , and u 4 ¯ = ( λ 7 − λ 8 ) b 2 H C 6 .

Proof

Through the optimal conditions ∂ ℋ ∂ u 1 = 0 , ∂ ℋ ∂ u 2 = 0 , ∂ ℋ ∂ u 3 = 0 , and ∂ ℋ ∂ u 4 = 0 , we obtain

∂ ℋ ∂ u 1 = C 3 u 1 − λ 1 ( ζ A + ρ I ) S N + λ 3 ( ζ A + ρ I ) S N ⇒ ∂ ℋ ∂ u 1 = C 3 u 1 − ( λ 1 − λ 3 ) ( ζ A + ρ I ) S N

∴ u 1 = ( λ 1 − λ 3 ) ( ζ A + ρ I ) S N C 3 = u 1 ¯

∂ ℋ ∂ u 2 = C 4 u 2 − λ 1 α S + λ 2 α S = 0

⇒ ∂ ℋ ∂ u 2 = C 4 u 2 − ( λ 1 − λ 2 ) α S = 0

∴ u 2 = ( λ 1 − λ 2 ) α S C 4 = u 2 ¯

∂ ℋ ∂ u 3 = C 5 u 3 − λ 5 b 1 I + λ 8 b 1 I = 0

⇒ ∂ ℋ ∂ u 3 = C 5 u 3 − ( λ 5 − λ 8 ) b 1 I = 0

∴ u 3 = ( λ 5 − λ 8 ) b 1 I C 5 = u 3 ¯

∂ ℋ ∂ u 4 = C 6 u 4 − λ 7 b 2 H + λ 8 b 2 H = 0

⇒ ∂ ℋ ∂ u 4 = C 6 u 4 − ( λ 7 − λ 8 ) b 2 H = 0

∴ u 4 = ( λ 7 − λ 8 ) b 2 H C 6 = u 4 ¯ .

Since u 1 , u 2 , u 3 , u 4 ∈ [ 0 , 1 ] , we have

u 1 = u 2 = u 3 = u 4 = 0 for u 1 ¯ < 0 , u 2 ¯ < 0 , u 3 ¯ < 0 , u 4 ¯ < 0 , u 1 = u 2 = u 3 = u 4 = 1 for u 1 ¯ > 1 , u 2 ¯ > 1 , u 3 ¯ > 1 , u 4 ¯ > 1 .

Otherwise, u 1 = u 1 ¯ , u 2 = u 2 ¯ , u 3 = u 3 ¯ , and u 4 = u 4 ¯

Hence, we got optimum values of J for these optimal controls.□

5.2 Numerical simulation of optimal control problem

The numerical simulations of developed optimal control problem are executed in MATLAB over a period of 400 days. Using the forward and backward fourth-order R-K technique, the system of equation (3) and adjoint equation (4) are computed. For this simulation, we used the parameter values mentioned in Table 2. Assume that C 1 = 1 , C 2 = 1 , C 3 = 1 , C 4 = 30 , C 5 = 40 , and C 6 = 40 . Suppose that b 1 and b 2 are both equal to 0.2. The numerical simulations employed four different control strategies:

Single control interventions.
Two control interventions.
Three control interventions.
All the four control interventions.

5.2.1 Single control interventions

The dynamics of the infected and disinfected classes are shown in Figure 8 with respect to single control interventions in the time range [ 0 , 400 ] , and their average values are shown in Table 3. Figure 8 and Table 3, both demonstrate that in comparison to other control measures, vaccination control u 2 ≠ 0 significantly reduces the number of COVID-19 infections.

Figure 8

Dynamics of infected classes and disinfected classes with respect to single control.

Table 3

The average values of vaccination, asymptomatic infected, symptomatic infected, quarantine, hospitalization, and recovered individuals with respect to single control interventions

Combination of controls	Avg. V	Avg. A	Avg. I	Avg. Q	Avg. H	Avg. R
u 1 ≠ 0	3205.158	475.471	575.954	652.580	679.108	3972.540
u 2 ≠ 0	3317.720	438.658	527.628	629.848	651.894	4381.719
u 3 ≠ 0	2784.821	507.182	599.407	738.492	741.587	3727.925
u 4 ≠ 0	2644.079	546.534	637.839	761.056	772.604	3591.438

5.2.2 Two control interventions

The dynamics of the infected and disinfected classes are shown in Figure 9 with respect to two control interventions in the time range [ 0 , 400 ] , and their average values are shown in Table 4. Both Figure 9 and Table 4 show that, compared to other combinations of control measures, the combination of controls u 1 ≠ 0 and u 2 ≠ 0 effectively lowers the number of COVID-19 infections.

Figure 9

Dynamics of infected classes and disinfected classes with respect to two controls.

Table 4

The average values of vaccination, asymptomatic infected, symptomatic infected, quarantined, hospitalized, and recovered individuals with respect to two control interventions

Combination of controls	Avg. V	Avg. A	Avg. I	Avg. Q	Avg. H	Avg. R
u 1 ≠ 0 , u 2 ≠ 0	4695.580	431.549	420.672	432.018	561.640	5701.664
u 2 ≠ 0 , u 3 ≠ 0	4401.601	458.480	461.482	410.647	614.079	5485.190
u 3 ≠ 0 , u 4 ≠ 0	3109.483	514.018	492.618	557.284	645.437	4361.615
u 1 ≠ 0 , u 4 ≠ 0	3271.927	521.374	501.906	524.049	658.728	4547.077
u 1 ≠ 0 , u 3 ≠ 0	3540.167	507.682	525.005	499.525	670.501	4681.542
u 2 ≠ 0 , u 4 ≠ 0	4147.281	482.167	449.097	469.438	592.649	5064.679

5.2.3 Three control interventions

The dynamics of the infected and disinfected classes are shown in Figure 10 with respect to three control interventions in the time range [ 0 , 400 ] , and their average values are shown in Table 5. The combination of controls u 1 ≠ 0 , u 2 ≠ 0 and u 3 ≠ 0 drastically decreases the number of COVID-19 infections individuals as compared to other combination control measures, as shown in Figure 10 and Table 5, respectively.

Figure 10

Dynamics of infected classes and disinfected classes with respect to three controls.

Table 5

The average values of vaccination, asymptomatic infected, symptomatic infected, quarantine, hospitalization, and recovered individuals with respect to three control interventions

Combination of controls	Avg. V	Avg. A	Avg. I	Avg. Q	Avg. H	Avg. R
u 1 ≠ 0 , u 2 ≠ 0 & u 3 ≠ 0	6425.225	258.742	245.404	297.057	374.409	7053.281
u 2 ≠ 0 , u 3 ≠ 0 & u 4 ≠ 0	5829.570	314.749	308.199	330.545	425.516	6473.164
u 1 ≠ 0 , u 3 ≠ 0 & u 4 ≠ 0	5240.372	362.030	357.472	389.672	492.044	6105.738
u 1 ≠ 0 , u 2 ≠ 0 & u 4 ≠ 0	6071.185	280.669	262.908	354.431	410.725	6842.019

5.2.4 Four control interventions

The dynamics of the infected and disinfected classes are shown in Figure 9 with respect to four control interventions in the time range [ 0 , 400 ] , and their average values are shown in Table 4. Figure 11 and Table 6 illustrate that the combination of four controls u 1 ≠ 0 , u 2 ≠ 0 , u 3 ≠ 0 , and u 4 ≠ 0 significantly reduces the number of COVID-19 infections compared to the no control intervention case.

Figure 11

Dynamics of infected classes and disinfected classes with respect to four controls.

Table 6

The average values of vaccination, asymptomatic infected, symptomatic infected, quarantine, hospitalization, and recovered individuals with respect to two control interventions

Combination of controls	Avg. V	Avg. A	Avg. I	Avg. Q	Avg. H	Avg. R
u 1 ≠ 0 , u 2 ≠ 0 ,	6954.349	127.627	151.408	137.715	139.373	7901.659
u 3 ≠ 0 & u 4 ≠ 0
u 1 = 0 , u 2 = 0 ,	2921.125	1363.193	656.145	926.809	951.694	2759.896
u 3 = 0 & u 4 = 0

Compared to the use of a single control strategy, a double control strategy, and a triple control strategy, the four control strategies significantly reduced the COVID-19 infections and enhanced recoveries.

6 Conclusion

In this work, by adding vaccination, quarantine, and hospitalization classes, we created a deterministic mathematical model for the dynamics of COVID-19 transmission. We examined the stability of the IFE, and the results imply that disease can be eliminated if the basic reproduction number is less than unity. Due to less surveillance and protection, such as mask use, social distance, and good sanitation, vaccination receivers are frequently reinfected. So that if the transition rate η of vaccinated population become exposed is zero, then the EE was became asymptotically unstable. This means no vaccinated individual will be reinfected and remain healthy with high surveillance. By fitting the proposed model to daily confirmed and cumulative confirmed COVID-19 cases in India, ten biological significant parameters are estimated and the basic reproduction number value computed as 1.971. The sensitive analysis results evaluated that the significance of vaccination for susceptible individuals as well as the quarantine and hospitalization of symptomatic infected individuals in reducing the COVID-19 disease transmission. The impact of vaccination on the infected individuals reveals that if vaccine effectiveness is higher, the disease reproduction will be reduced and possibly even eliminated. The optimal control system was developed by using both pharmaceutical and non-pharmaceutical control strategies such as preventive control u 1 ( t ) related to public awareness programs, using a face mask, social distance etc., vaccine control u 2 ( t ) related to continuous vaccination for susceptible individuals and treatment controls u 3 ( t ) and u 4 ( t ) related to the most effective therapy, such as the use of immunomodulators like intraneural facilitation (INF) or zinc to boost the immunological response, or any prescribed antiviral drugs like remdesivir, arbidol, etc. for symptomatic infected and hospitalized individuals. Numerical simulation of optimal control problem demonstrated that in comparison to combinations of three controls strategies , two controls strategies as well as single control strategies , combination of four controls strategy is more effective in mitigating the spread of COVID-19 disease.

Acknowledgements

We would like to thank reviewers for taking the time and effort necessary to review the manuscript. We sincerely appreciate all valuable comments and suggestions, which help us to improve the quality of the manuscript.

Funding information: This research received no specific grant from any funding agency, commercial or nonprofit sectors.
Author contributions: Mallela Ankamma Rao: conceptualization, writing – original draft, software, investigation, validation and writing – review and editing. Murugadoss Prakash Raj and Karuppusamy Arun Kumar: software, methodology, and validation and writing. Ambalarajan Venkatesh: conceptualization, writing – original draft, validation and review and editing. D. K. K. Vamsi: conceptualization of optimal control studies and structuring the draft and editing.
Conflict of interest: The authors have no conflicts of interest to disclose.
Ethical approval: This research did not required ethical approval.
Data availability statement: This manuscript has associated data in a data repository. [Authors comment: The data used to carry out this study are openly available at https://www.covid19india.org/ [42]].

References

[1] Acuña-Zegarra, M. A., Díaz-Infante, S., Baca-Carrasco, D., & Olmos-Liceaga, D. (2021). COVID-19 optimal vaccination policies: A modeling study on efficacy, natural and vaccine-induced immunity responses. Mathematical Biosciences, 337, 108614. 10.1016/j.mbs.2021.108614Search in Google Scholar PubMed PubMed Central

[2] Algarni, A. D., Ben Hamed, A., Hamdi, M., Elmannai, H., & Meshoul, S. (2022). Mathematical COVID-19 model with vaccination: a case study in Saudi Arabia. PeerJ Computer Science, 8, e959. 10.7717/peerj-cs.959Search in Google Scholar PubMed PubMed Central

[3] Ankamma Rao, M., & Venkatesh, A. (2023). SEAIQHRDP mathematical model Analysis for the transmission dynamics of COVID-19 in India. Journal of Computational Analysis and Applications, 31(1), 96–116. Search in Google Scholar

[4] Carvalho, S. A., da Silva, S. O., & Charret, I. D. C. (2019). Mathematical modeling of dengue epidemic: control methods and vaccination strategies. Theory in Biosciences, 138(2), 223–239. 10.1007/s12064-019-00273-7Search in Google Scholar PubMed

[5] Castillo-Chavez, C., & Song, B. (2004). Dynamical models of tuberculosis and their applications. Mathematical Biosciences and Engineering, 1(2), 361–404. 10.3934/mbe.2004.1.361Search in Google Scholar PubMed

[6] Chhetri, B., Bhagat, V. M., Vamsi, D., Ananth, V., Prakash D. B., Mandale R., …, Sanjeevi C. B. Within-host mathematical modeling on crucial inflammatory mediators and drug interventions in COVID-19 identifies combination therapy to be most effective and optimal. (2021). Alexandria Engineering Journal, 60(2), 2491–2512. 10.1016/j.aej.2020.12.011Search in Google Scholar

[7] Chhetri, B., Vamsi, D. K. K., Prakash, D. B., Balasubramanian, S., & Sanjeevi, C. B. (2022). Age structured mathematical modeling studies on COVID-19 with respect to combined vaccination and medical treatment strategies. Computational and Mathematical Biophysics, 10(1), 281–303. 10.1515/cmb-2022-0143Search in Google Scholar

[8] Das, P., Nadim, S.S., & Das, S. (2021). Dynamics of COVID-19 transmission with comorbidity: a data driven modelling based approach. Nonlinear Dynamics, 106, 1197–1211. 10.1007/s11071-021-06324-3Search in Google Scholar PubMed PubMed Central

[9] Chitnis, N., Hyman, J. M., & Cushing, J. M. (2008). Determining important parameters in the spread of Malaria through the sensitivity analysis of a mathematical model. Bulletin of Mathematical Biology, 70(5), 1272–1296. 10.1007/s11538-008-9299-0Search in Google Scholar PubMed

[10] Diagne, M. L., Rwezaura, H., Tchoumi, S. Y., & Tchuenche, J. M. (2021). A mathematical model of COVID-19 with vaccination and treatment. Computational and Mathematical Methods in Medicine, 2021, 1–16. 10.1155/2021/1250129Search in Google Scholar PubMed PubMed Central

[11] Gantmacher, F.R., & Mirsky, L. (1961). The theory of matrices. The Mathematical Gazette, 45(353), 267–267. 10.2307/3612823Search in Google Scholar

[12] Gumel, A. B., Ruan, S., Day, T., Watmough, J., Brauer, F., van den Driessche, P., …, Sahai, B. M. (2004). Modelling strategies for controlling SARS outbreaks. Proceedings of the Royal Society of London. Series B: Biological Sciences, 271(1554), 2223–2232. 10.1098/rspb.2004.2800Search in Google Scholar PubMed PubMed Central

[13] Imai, N., Cori, A., Dorigatti, I., Baguelin, M., Christl, A., Riley, S., & Ferguson, M. (2020). Report 3: Transmissibility of 2019-nCoV. https://www.imperial.ac.uk/mrc-global-infectious-disease-analysis/news-wuhan-coronavirus/. Search in Google Scholar

[14] Klaus, D., & Heesterbeek, J. (2002). Daniel Bernoulli’s epidemiological model revisited. Mathematical Biosciences, 180(1–2), 1–21. 10.1016/S0025-5564(02)00122-0Search in Google Scholar

[15] Khajanchi, S., Sarkar, K., Mondal, J., Nisar, K. S., & Abdelwahab, S. F. (2021). Mathematical modeling of the COVID-19 pandemic with intervention strategies. Results in Physics, 25, 104285. 10.1016/j.rinp.2021.104285Search in Google Scholar PubMed PubMed Central

[16] Khajanchi, S., Sarkar, K., & Banerjee, S. (2022). Modeling the dynamics of COVID-19 pandemic with implementation of intervention strategies. The European Physical Journal Plus, 137(1), 129. 10.1140/epjp/s13360-022-02347-wSearch in Google Scholar PubMed PubMed Central

[17] Kurmi, S., & Chouhan, U. (2022). A multicompartment mathematical model to study the dynamic behaviour of COVID-19 using vaccination as control parameter. Nonlinear Dynamics, 109(3), 2185–2201. 10.1007/s11071-022-07591-4Search in Google Scholar PubMed PubMed Central

[18] Lenhart, S., & Workman, J. T. (2007). Optimal control applied to biological models. (1st ed.). Chapman and Hall/CRC. https://doi.org/10.1201/9781420011418. 10.1201/9781420011418Search in Google Scholar

[19] Lienart, A., Chipart, H.(1914)Sur le signe de la partie reelle des racines daune equation algebrique. J. Math. Pures Appl, 10, 291–346. Search in Google Scholar

[20] Li, Q., Guan, X., Wu, P., Wang, X., Zhou, L., Tong, Y., …, Feng, Z. (2020). Early transmission dynamics in Wuhan, China, of novel coronavirus-infected pneumonia. New England Journal of Medicine, 382(13), 1199–1207. 10.1056/NEJMoa2001316Search in Google Scholar PubMed PubMed Central

[21] Mandal, M., Jana, S., Nandi, S. K., Khatua, A., Adak, S., & Kar, T. (2020). A model based study on the dynamics of COVID-19: Prediction and control. Chaos, Solitons and Fractals, 136, 109889. 10.1016/j.chaos.2020.109889Search in Google Scholar PubMed PubMed Central

[22] Marino, S., Hogue, I. B., Ray, C. J., & Kirschner, D. E. (2008). A methodology for performing global uncertainty and sensitivity analysis in systems biology. Journal of Theoretical Biology, 254(1), 178–196. 10.1016/j.jtbi.2008.04.011Search in Google Scholar PubMed PubMed Central

[23] Nadim, S. S., Ghosh, I., & Chattopadhyay, J. (2021). Short-term predictions and prevention strategies for COVID-19: A model-based study. Applied Mathematics and Computation, 404, 126251. 10.1016/j.amc.2021.126251Search in Google Scholar PubMed PubMed Central

[24] Ndaïrou, F., Area, I., Nieto, J. J., & Torres, D. F. (2020June). Mathematical modeling of COVID-19 transmission dynamics with a case study of Wuhan. Chaos Solitons Fractals, 135, 109846. 10.1016/j.chaos.2020.109846Search in Google Scholar PubMed PubMed Central

[25] Negi, S. S., Rana, P. S., Sharma, N., & Khatri, M. S. (2022). A novel SEIAHR compartment model for accessing the impact of vaccination, intervention policies, and quarantine on the COVID-19 pandemic: a case study of most affected countries Brazil, India, Italy, and USA. Computational and Applied Mathematics, 41(7) , 305. 10.1007/s40314-022-01993-1Search in Google Scholar

[26] Ngonghala, C. N., Iboi, E., Eikenberry, S., Scotch, M., MacIntyre, C. R., Bonds, M. H., & Gumel, A. B. (2020). Mathematical assessment of the impact of non-pharmaceutical interventions on curtailing the 2019 novel Coronavirus. Mathematical Biosciences, 325, 108364. 10.1016/j.mbs.2020.108364Search in Google Scholar PubMed PubMed Central

[27] Okuonghae, D., & Omame, A. (2020). Analysis of a mathematical model for COVID-19 population dynamics in Lagos, Nigeria. Chaos Solitons Fractals, 139, 110032. 10.1016/j.chaos.2020.110032Search in Google Scholar PubMed PubMed Central

[28] Open government data (OGD) platform India or data.gov.in. 2020, https://data.gov.in/resources/crude-death-rate-india-2011. Search in Google Scholar

[29] Pontryagin, L. (1987). Mathematical theory of optimal processes. Classics of Soviet Mathematics. CRC Press. Search in Google Scholar

[30] Prathumwan, D., Trachoo, K., & Chaiya, I. (2020). Mathematical modeling for prediction dynamics of the coronavirus disease 2019 (COVID-19) pandemic, quarantine control measures. Symmetry, 12(9), 1404. 10.3390/sym12091404Search in Google Scholar

[31] Rai, R. K., Khajanchi, S., Tiwari, P. K., Venturino, E., & Misra, A. K. (2021). Impact of social media advertisements on the transmission dynamics of COVID-19 pandemic in India. Journal of Applied Mathematics and Computing, 68(1), 19–44. 10.1007/s12190-021-01507-ySearch in Google Scholar PubMed PubMed Central

[32] Rodrigues, H. S., Monteiro, M. T. T., & Torres, D. F. M. (2013). Sensitivity analysis in a dengue epidemiological model. Conference Papers in Mathematics, 2013, 1–7. 10.1155/2013/721406Search in Google Scholar

[33] Sardar, T., Nadim, S. S., Rana, S., & Chattopadhyay, J. (2020). Assessment of lockdown effect in some states and overall India: A predictive mathematical study on COVID-19 outbreak. Chaos, Solitons and Fractals, 139, 110078. 10.1016/j.chaos.2020.110078Search in Google Scholar PubMed PubMed Central

[34] Shyamsunder, Bhatter, S., Jangid, K., Abidemi, A., Owolabi, K., & Purohit, S. (2023). A new fractional mathematical model to study the impact of vaccination on COVID-19 outbreaks. Decision Analytics Journal, 6, 100156. 10.1016/j.dajour.2022.100156Search in Google Scholar

[35] Srivastav, A. K., Tiwari, P. K., Srivastava, P. K., Ghosh, M., & Kang, Y. (2021). A mathematical model for the impacts of face mask, hospitalization and quarantine on the dynamics of COVID-19 in India: deterministic vs. stochastic. Mathematical Biosciences and Engineering, 18(1), 182–213. 10.3934/mbe.2021010Search in Google Scholar PubMed

[36] Van den Driessche, P., & Watmough, J. (2002). Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Mathematical Biosciences, 180(1–2), 29–48. 10.1016/S0025-5564(02)00108-6Search in Google Scholar PubMed

[37] Venkatesh, A., Ankamma Rao, M., & Vamsi, D. K. K.(2023). A comprehensive study of optimal control model simulation for COVID-19 infection with respect to multiple variants. Communications in Mathematical Biology and Neuroscience, Article ID 75, https://doi.org/10.28919/cmbn/8031. Search in Google Scholar

[38] Venkatesh, A., & Ankamma Rao, M. (2023). Mathematical Model for COVID-19 pandemic with implementation of intervention strategies and cost-effectiveness analysis, Results in Control and Optimization, 14, 100345. 10.1016/j.rico.2023.100345Search in Google Scholar

[39] Watson, O. J., Barnsley, G., Toor, J., Hogan, A. B., Winskill, P., & Ghani, A. C. (2022. September). Global impact of the first year of COVID-19 vaccination: a mathematical modelling study. The Lancet Infectious Diseases, 22(9), 1293–1302. 10.1016/S1473-3099(22)00320-6Search in Google Scholar PubMed PubMed Central

[40] WHO Coronavirus (COVID-19) dashboard. (2020). WHO Coronavirus (COVID-19) Dashboard With Vaccination Data. https://covid19.who.int/. Search in Google Scholar

[41] Worldometer, (COVID-19) CORONAVIRUS PANDEMIC. 2020. https://worldometers.info/coronavirus/countries/. Search in Google Scholar

[42] Coronavirus in India: Latest Map and Case Count. (2020). https://www.covid19india.org. Search in Google Scholar

Received: 2023-01-05

Revised: 2023-12-04

Accepted: 2023-12-04

Published Online: 2024-02-29

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/cmb-2023-0112

Keywords for this article

COVID-19; epidemic model; vaccination; bifurcation; sensitivity analysis; optimal control problem

Creative Commons

BY 4.0