Optimal control and bifurcation analysis of SEIHR model for COVID-19 with vaccination strategies and mask efficiency

Poosan Moopanar Muthu; Anagandula Praveen Kumar

doi:10.1515/cmb-2023-0113

Article Open Access

Optimal control and bifurcation analysis of SEIHR model for COVID-19 with vaccination strategies and mask efficiency

Poosan Moopanar Muthu and Anagandula Praveen Kumar

Published/Copyright: March 9, 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 article, we present a susceptible, exposed, infected, hospitalized and recovered compartmental model for COVID-19 with vaccination strategies and mask efficiency. Initially, we established the positivity and boundedness of the solutions to ensure realistic predictions. To assess the epidemiological relevance of the system, an examination is conducted to ascertain the local stability of the endemic equilibrium and the global stability across two equilibrium points are carried out. The global stability of the system is demonstrated using Lyapunov’s direct method. The disease-free equilibrium is globally asymptotically stable when the basic reproduction number (BRN) is less than one, whereas the endemic equilibrium is globally asymptotically stable when BRN is greater than one. A sensitivity analysis is performed to identify the influential factors in the BRN. The impact of various time-dependent strategies for managing and regulating the dynamic transmission of COVID-19 is investigated. In this study, Pontryagin’s maximum principle for optimal control analysis is used to identify the most effective strategy for controlling the disease, including single, coupled, and threefold interventions. Single-control interventions reveal physical distancing as the most effective strategy, coupled measures reduce exposed populations, and implementing all controls reduces susceptibility and infections.

Keywords: global stability; backward bifurcation; optimal control; Lyapunov’s direct method

MSC 2010: 92B05; 92D30; 93C15; 93D30; 34D23

1 Introduction

The world is continuously grappled by COVID-19 with different variants of concern. Scientific knowledge has been developed to gain a profound understanding of the intricacies encompassing spread patterns, the degree of severity, clinical characteristics, and the factors predisposing individuals to the infection. This formidable ailment, renowned for its highly contagious nature, stems from the pathogenic SARS-CoV-2 virus, primarily having an impact on the respiratory system [38]. The rapidity of its transmission led to the eventual classification of the outbreak as a pandemic [38].

The mathematical model is a relevant tool for assessing disease evolution. These models can aid in understanding how a disease is spread, forecasting the course of an outbreak, and evaluating the effectiveness of various therapies. Out of the several mathematical models, one common and important type is the compartmental model, which categorizes the population as susceptible, infected, and recovered or deceased depending on their disease condition. These models are useful for researchers to predict the transmission rate or the effectiveness of interventions and how treatments will affect the spread of the disease [20].

Zeb et al. [35] discussed the behavior of COVID-19 by proposing an isolation class and have shown that person-to-person contact is the possible root of the outbreak. They suggested that the seclusion of infected individuals can curtail the spread of infection to others. Arcede et al. [1] understood the scenario of COVID-19 and divided the infected population into two distinct classes, namely symptomatic and asymptomatic. Those who are symptomatic will move to the treatment compartment. Since there are no serological tests available to assess asymptomatic patients, no treatment is administered to them. With the data on confirmed cases and fatalities from the Philippines, Italy, the United Kingdom, China, and the United States, the model is tested. They found that partial or total restriction of population may reduce the transmission of disease. Cooper et al. [6] constructed a simple susceptible, infected, and recovered (SIR) model to discuss the spread of COVID-19 within the community. The model provides valuable predictions such as the number of susceptible, infected, and recovered persons in major countries all over the world.

Youssef et al. [34] developed an susceptible, exposed, infected, and recovered (SEIR) model to study the number of infected persons in any country. They used the data from Saudi Arabia to verify the correctness of the model and urged the Saudi Arabians to spend as much time at home as they could. Kumari et al. [16] designed a new susceptible, exposed, infected, asymptomatic, quarantined, recovered, deaths, and insusceptible model to describe overall disease dynamics. Novel models have been introduced to improve disease spread predictions, surpassing the accuracy of previous ones like SIR, SEIR, and susceptible, exposed, infected, quarantined, and recovered models. The diligent efforts of researchers have led to refined models that have been used to estimate the pandemic’s impact on India and its worst-hit states, enhancing our understanding of disease dynamics and guiding effective containment measures.

Srivastav et al. [32] proposed a model to study COVID-19, considering face masks, hospitalizations, and quarantine of asymptomatic individuals. They have taken the data of the most affected states in India to forecast the infection rate of susceptible and symptomatic populations and the rate of recovery in quarantined individuals [37]. Peter et al. [24] composed a model to implement social distancing, employing the use of face masks and hand sanitizer. One of the important aspects is that recovered people are going to be susceptible once their immunity is over. To verify the formulated model, they fitted it with COVID-19 data from Pakistan.

There are researchers who have focused on the seasonal spread of COVID-19. Li et al. [18] developed a time-periodic compartmental model to study the seasonal fluctuations in many countries. Seasonal pattern mainly focuses on public gatherings, international travel, and festivals. They showed that the positive periodic solution is globally asymptotically stable. This means that the solution consistently approached a stable state regardless of the initial conditions. Kamrujjaman et al. [10] formulated a simple SEIR model to speak mainly on subtle features like panic, tension, and anxiety of COVID-19. Although human behavior is not permanent to discuss the dynamics, the statistics show that 5% of people in Italy reported anxiety [26]. They compared the results with the data of COVID-19 in Italy and showed a good fit to predict the number of infected individuals. Three distinct studies, led by Umdekar et al. [33] and Sharma and Sharma [28,29], delve into critical aspects of epidemic modeling. Umdekar et al. [33] emphasize the Holling Type II treatment function’s effectiveness in reducing infections, highlighting its superiority over traditional methods and the importance of a higher cure rate for disease eradication.

Optimal control analysis is an indispensable tool in the field of disease transmission dynamics, offering powerful insights into the design and implementation of effective intervention measures. It provides a framework for optimizing the allocation of resources, timing of interventions, and intensity of control measures to minimize disease transmission and mitigate its impact on public health. By leveraging mathematical modeling, optimization techniques, and advanced analytical methods, optimal control strategies offer a comprehensive approach to understand the intricacies of pathogen propagation dynamics and devising evidence-based strategies for disease control [13].

Madubueze et al. [19] studied the optimal control analysis of the SEQIJR model for COVID-19 by including public health education, quarantining the population, and isolation of infectious but not hospitalized population. The quarantine compartment is designed to accommodate individuals arriving from high-risk regions, which provides a crucial measure to prevent the spread of disease. Previous studies have demonstrated that in the absence of quarantine for exposed individuals, an infected individual can disseminate the disease to two or more individuals, leading to a rapid transmission rate. However, the introduction of quarantine measures has shown promising results in halting the transmission of the disease.

Bandekar et al. [2] developed a susceptible, exposed, infected, hospitalized and recovered (SEIHR) model to describe all the concepts related to COVID-19 like face mask efficiency, mask compliance, and vaccination efficiency. They included the reinfection of recovered individuals who will move to the infected class after waning of immunity. They fitted the COVID-19 data of India with that proposed model, and they predicted the number of infected persons until August 2022.

In this manuscript, we extend the analysis of the COVID-19 model [2] to include global stability, bifurcation, and optimal control analysis. In global stability, we will discuss the long-term behavior of the disease dynamics by using Lyapunov’s function. Bifurcation analysis provides a deeper understanding of the complex dynamics exhibited by disease, allowing the exploration of a wide range of scenarios and parameters. Sensitivity analysis is also performed to choose the biological indicators that are the most sensitive so as to examine the importance of constant controls. Optimal control strategies are helping to scale down the severity of the disease by employing Pontryagin’s maximum principle (PMP).

The article is organized as follows: Section 2 elaborates on the mathematical formulation of the model; Section 3 examines the local and global stability at equilibrium points; Section 4 undertakes the analysis of the backward bifurcation of the model; Section 5 delves into the sensitivity analysis of the model’s parameters; Section 6 conducts an optimal control analysis of the model, and lastly, Section 7 provides a comprehensive discussion and conclusions summarizing the article’s content.

2 Mathematical formulation

The formulated SEIHR compartmental model divides the total population (N) into five distinct compartments: These compartments represent the five stages that individuals can go through during their experience with COVID-19.

S: Individuals susceptible to contracting the disease.
E: Individuals who have been infected but are not yet capable of transmitting the disease.
I: Infectious individuals who can transmit the disease to others.
H: Infective in-care individuals who are getting care either at home or in isolation due to the severity of their infection.
R: Recovered individuals who have successfully overcome the disease and acquired immunity.

To develop a mathematical model, there are certain assumptions need to be noted.

The persons are being recruited to the susceptible class S at a continuous pace Λ . According to the presumption that disease spreads through simple mass action, the contact frequency between infected class I and susceptible class S is directly proportional to the entire population.

In this model, β 1 S I represents the rate at which susceptible individuals acquire infections. It takes into account the transmission rate ( β ), the effectiveness of mask efficiency ( ε m ), mask compliance ( c m ), and the effectiveness of vaccination ( ε v ). The term S I represents the number of potential interactions between susceptible and infected individuals that could lead to new infections; ( 1 − ε m c m ) accounts for the reduction in the transmission rate due to individuals wearing masks, and the population is compliance with mask usage. The reduction in infection rate due to vaccination of people is represented by ( 1 − ε v ) . Next, the β 1 S H represents the rate at which exposed individuals progress to the infective in care, taking into account the reduced chance of infection from individuals in the H class by modification factor ϕ < 1 . This reduction is due to the assumption that individuals receiving care are relatively protected and have a lower chance of transmitting the infection compared to those without care.

The population is moving to infected class I from exposed class E at a rate ν . If the infection is severe, then they will move to infective in-care class H with a screening rate of θ . Those who belong to the infected class I are going to class H at a rate of γ . Individuals from H class who are critically ill may die from their illnesses at a rate of ψ , and those who survive will be transferred to class R at a rate δ . Those who are recovered may again get in touch with I or H classes, become reinfected, and join infected class I . Due to strong immunity, the reinfection rate is very low ( η < 1 ).

Mask efficiency ( ε m ) is the effectiveness of masks in preventing the transmission of infection, with values in the range of [0, 1], and it represents the fraction of transmission that masks can effectively block. A value of 0 indicates that masks have no effect, while a value of 1 indicates that masks are fully effective at preventing transmission.

Mask compliance ( c m ) means average community behavior for using masks with values in the range of [0, 1], and it represents the proportion of the population that is compliant with wearing masks. A value of 0 indicates no compliance (nobody wears masks), while a value of 1 indicates full compliance (everyone wears masks).

Vaccination efficiency ( ε v ) represents, the effectiveness of vaccination in preventing the transmission of infection, with values in the range of [0, 1], and it represents the fraction of transmission that vaccination can effectively prevent. A value of 0 indicates that vaccination has no effect, while a value of 1 indicates that vaccination is fully effective at preventing transmission.

The information provided at the above points has been visually represented in Figure 1, and their corresponding mathematical equations are given in equation (1). Figure 1 presents the dynamics of these population transitions and interactions among the different classes, providing a comprehensive understanding of the disease spread and care given in model (1).

Figure 1

Transmission figure of SEIHR model.

The SEIHR model is delineated by the subsequent system of equations:

(1) d S d t = Λ − β 1 S I − β 2 S H − μ S d E d t = β 1 S I + β 2 S H − β 4 E d I d t = ν E − β 6 I + β 3 R ( I + H ) d H d t = θ E + γ I − β 5 H d R d t = δ H − μ R − β 3 R ( I + H ) ,

where

β 1 = β ( 1 − ε m c m ) ( 1 − ε v ) β 2 = β ϕ ( 1 − ε m c m ) β 3 = β η ( 1 − ε m c m ) β 4 = θ + μ + ν β 5 = ψ + μ + δ β 6 = μ + γ .

3 Analysis of the model

3.1 Positivity analysis

It is imperative to establish the essentiality of showcasing the existence of positive and boundedness of solutions for the system of equation (1), thereby validating the mathematical model at hand.

Bandekar et al. [2] conducted an analysis on positivity and boundedness in relation to the model presented in this study.

We now demonstrate that the positive outcome persists throughout all future instances if the initial conditions for system (1) are positive. Using equation (1), we obtain

(2) d S d t S = 0 = Λ ≥ 0 , d E d t E = 0 = β 1 S I + β 2 S H ≥ 0 , d I d t I = 0 = v E + β 3 R H ≥ 0 d H d t H = 0 = θ E + γ I ≥ 0 , d R d t R = 0 = δ H ≥ 0 .

According to the above equations, the rates exhibit nonnegativity along the bounding planes of the nonnegative domain of R 5 defined by S = 0 , E = 0 , I = 0 , H = 0 , and R = 0 . Therefore, if a solution starts inside of this region, it will stay there for the entire duration of time t > 0 . This occurs because, as shown by the earlier inequalities, the direction of the vector field at the boundaries is always inward. Consequently, we draw conclusion that, as long as the initial conditions are positive, all of the system’s solutions are positive at any time t > 0 .

3.2 Boundedness of solutions

Theorem 1

The closed region

(3) Ω = ( S , E , I , H , R ) ∣ S , E , I , H , R ≥ 0 a n d 0 < N ( t ) ≤ Λ μ

is bounded set for system (1).

We have, N ( t ) = S ( t ) + E ( t ) + I ( t ) + H ( t ) + R ( t ) ,

(4) d N d t = d S d t + d E d t + d I d t + d H d t + d R d t d N d t = Λ − μ N − ψ H .

Integrating equation (4), we obtain

(5) N ( t ) = Λ μ − ψ e − μ t ∫ H ( y ) e μ y d y .

We obtain

(6) limsup t → ∞ N ( t ) ≤ Λ μ .

Hence, all solutions S ( t ) , E ( t ) , I ( t ) , H ( t ) , and R ( t ) of equation (1) are constrained within the bound of Λ μ . Consequently, this implies that N ( t ) is also bounded. Hence proved.

3.3 Existence of equilibrium points

In order to find the equilibrium points of equation (1), we are solving the following equations:

(7) d S d t = 0 , d E d t = 0 , d I d t = 0 , d H d t = 0 , and d R d t = 0 .

3.3.1 Disease-free equilibrium point ( E 0 )

The disease-free equilibrium point for model (1) is as follows:

(8) E 0 = ( S 0 , E 0 , I 0 , H 0 , R 0 ) = Λ μ , 0 , 0 , 0 , 0 .

3.3.2 Endemic equilibrium point ( E 1 )

From model (1), we obtain E 1 = ( S * , E * , I * , H * , R * ) . The components of the E 1 are

(9) S * = Λ β 1 I * + β 2 H * + μ H * = ( θ E * + γ I * ) β 5 , R * = δ H * μ + β 3 ( I * + H * ) ,

where E * and I * are derived from the following equations:

(10) B 1 I * 2 + B 2 E * 2 + B 3 E * I * + B 4 E * − B 5 I * = 0 ,

(11) K 1 E * 2 + K 2 E * I * + K 3 E * − K 4 I * = 0 ,

where

B 1 = β 3 ( β 5 + γ ) ( δ γ − β 5 β 6 ) , B 2 = β 3 [ δ θ 2 + β 5 θ ν ] B 3 = β 3 [ δ θ ( β 5 + 2 γ ) + ν ( β 5 + γ ) β 5 − β 5 β 6 θ ] , B 4 = μ β 5 2 ν B 5 = μ β 6 β 5 2 , K 1 = θ β 2 β 4 K 2 = β 1 β 4 [ ϕ γ + ( 1 − ε v ) β 5 ] , K 3 = μ β 4 β 5 − Λ β 2 θ K 4 = Λ β 2 γ + Λ β 1 β 5 .

The existence of an endemic equilibrium point is proved by Bandekar et al. [2].

3.4 Basic reproduction number (BRN, ℛ 0 )

The BRN is a fundamental parameter used to analyze the spread of infectious diseases. It establishes whether the infection will disappear over time or remain prevalent in the community. In a community where everyone is susceptible, it is described as the number of secondary transmissions caused by one infected individual and is indicated by ℛ 0 .

We are finding the next generation matrix (NGM) by using disease free equilibrium (DFE), and the largest eigenvalue of the matrix is said to be ℛ 0 [9].

From the definition of NGM, we will calculate ℱ and V as follows:

ℱ = β 1 S I + β 2 S H β 3 R ( I + H ) 0 and V = β 4 E − ν E + β 6 I − θ E − γ I + β 5 H .

Now,

F = Jacobian of ℱ at E 0 = 0 β 1 Λ μ β 2 Λ μ 0 0 0 0 0 0

V = Jacobian of V at E 0 = β 4 0 0 − ν β 6 0 − θ γ β 5

and

F V − 1 = b 11 b 12 b 13 0 0 0 0 0 0 ,

where

b 11 = Λ [ ν β 1 β 5 + β 2 γ ν + θ β 2 β 6 ] μ β 4 β 5 β 6 b 12 = Λ β 1 ( β 5 + β 2 γ ) μ β 4 β 5 β 6 b 13 = β 2 Λ μ β 4 β 5 β 6 .

Hence, ℛ 0 is the spectral radius of F V − 1 and is given as follows:

(12) ℛ 0 = Λ [ ν β 1 β 5 + β 2 γ ν + θ β 2 β 6 ] μ β 4 β 5 β 6 .

3.5 Local stability analysis

In our study, we investigate the local stability of model (1) at the endemic state, considering the ongoing presence of COVID-19 in the population.

The Jacobian of equation (1) is computed at E 1 , yielding the following expression:

(13) D E 1 = c 11 c 12 c 13 c 14 c 15 c 21 c 22 c 23 c 24 c 25 c 31 c 32 c 33 c 34 c 35 c 41 c 42 c 43 c 44 c 45 c 51 c 52 c 53 c 54 c 55 ,

where

c 11 = − β 1 I * − β 2 H * − μ , c 12 = 0 , c 13 = − β 1 S * , c 14 = − β 2 S * , c 15 = 0 , c 21 = β 1 I * + β 2 H * , c 22 = − β 4 , c 23 = β 1 S * , c 24 = β 2 S * , c 25 = 0 , c 31 = 0 , c 32 = ν , c 33 = − β 6 + β 3 R * , c 34 = β 3 R * , c 35 = β 3 R * ( I * + H * ) , c 41 = 0 , c 42 = θ , c 43 = γ , c 44 = − β 5 , c 45 = 0 , c 51 = 0 , c 52 = 0 , c 53 = − β 3 R * , c 54 = δ − β 3 R * , c 55 = − β 3 ( I * + H * ) .

The eigenvalues of Jacobian matrix ( D E 1 ) are calculated as follows:

(14) q ( λ ) ≔ λ 5 + D 1 λ 4 + D 2 λ 3 + D 3 λ 2 + D 4 λ + D 5 = 0 .

From Lienard-Chipart condition [8] for the polynomial q ( λ ) ,

if D 2 > 0 , D 4 > 0 , and

Δ 2 = D 1 D 3 1 D 2 = D 1 D 2 − D 3 > 0 Δ 4 = D 1 D 3 D 5 0 1 D 2 D 4 0 0 D 1 D 3 D 5 0 1 D 2 D 4 = D 1 D 2 D 3 D 4 − D 3 2 D 4 − D 1 2 D 4 2 − D 1 D 2 2 D 5 + D 2 D 3 D 5 + 2 D 1 D 4 D 5 − D 5 2 > 0 ,

then the E 1 is l.a.s.

The values of D 1 , D 2 , D 3 , D 4 , and D 5 are given in Appendix .

3.6 Global stability analysis

In epidemiological models, global stability is a desired quality since it signifies the long-term control or elimination of the illness. It means that the disease can be eliminated and will not reemerge if effective efforts are taken to stop the initial introduction or transmission of the illness.

The global stability of the E 0 of system (1) is given by the following theorem.

Theorem 2

If ℛ 0 ≤ 1 , then the disease-free state E 0 of system (1) is globally asymptotically stable (g.a.s) on Ω .

Proof

We are proceeding to construct the ensuing Lyapunov function,

(15) L 1 ( S ) = S − S * − S * log S S * ,

where S * = Λ μ corresponding to the E 0 .

We obtain L 1 = 0 if S = S * .

By differentiating equation (15) w.r.t t , we obtain

(16) d L 1 d t = S ′ − S * S S ′ = 1 − S * S S ′ = 1 − S * S [ Λ − β 1 S I − β 2 S H − μ S ] = μ S − Λ μ S [ Λ − β 1 S I − β 2 S H − μ S ] = − ( Λ − μ S ) 2 μ S − 1 − Λ μ S β 1 S I − 1 − Λ μ S β 2 S H ≤ 0 .

By using the Lyapunov theorem, E 0 is g.a.s.□

Next, the global stability of the E 1 of system (1) is provided.

Theorem 3

Let ℛ 0 > 1 , whenever the sign ( N − N * ) = sign ( H − H * ) , the endemic state E 1 of system (1) is g.a.s on Ω .

Proof

In order to demonstrate the global stability of E 1 , we shall endeavor to devise the ensuing Lyapunov function.

(17) L 2 [ S , E , I , H , R ] = ( S + E + I + H + R ) − ( S * + E * + I * + H * + R * ) − ( S * + E * + I * + H * + R * ) log S + E + I + H + R S * + E * + I * + H * + R * = N − N * − N * log N N * .

By computing the time derivative of L 2 along the solutions of system (1), we derive the following expression:

(18) d L 2 d t = N ′ 1 − N * N = ( Λ − μ N − ψ H ) 1 − N * N

We have Λ = μ N * + ψ H *

(19) d L 2 d t = 1 − N * N [ μ N * + ψ H * − μ N − ψ H ] = − μ ( N − N * ) 2 N − ψ ( N − N * ) ( H − H * ) N ≤ 0 .

Therefore, it can be observed that d L 2 d t is negative and d L 2 d t = 0 iff S = S * , E = E * , I = I * , H = H * , and R = R * within the domain Ω . As a result, the singleton set { E 1 } is the greatest positively invariant set contained within { ( S , E , I , H , R ) ∈ Ω : L 2 = 0 } . By virtue of the Lyapunov-Lasalle theorem, we can conclude that the E 1 is g.a.s.□

4 Backward bifurcation analysis

To determine the conditions for a backward bifurcation, we utilize the central manifold theory and normal forms. Our approach involves employing a well-known result from published resources. Let x 1 = S , x 2 = E , x 3 = I , x 4 = H , and x 5 = R . Consider the following system

(20) f 1 = Λ − β 1 x 1 x 3 − β 2 x 1 x 4 − μ x 1 f 2 = β 1 x 1 x 3 + β 2 x 1 x 4 − β 4 x 2 f 3 = ν x 2 − β 6 x 3 + β 3 x 5 ( x 3 + x 4 ) f 4 = θ x 2 + γ x 3 − β 5 x 4 f 5 = δ x 4 − μ x 5 − β 3 x 5 ( x 3 + x 4 ) .

From system (1), we observe that β is most influential factor on BRN ( ℛ 0 ). Henceforth, we consider β as the bifurcation parameter.

Hence, when ℛ 0 = 1 and β = β * , we have

(21) β * = μ β 4 β 5 Λ ( ( 1 − ε m c m ) ( 1 − ε v ) ν β 5 + ( 1 − ε m c m ) ϕ γ ν + ( 1 − ε m c m ) ϕ θ β 6 ) .

The Jacobian matrix assessed at the disease-free state ( E 0 ) and β = β * for system (1) is given as follows:

(22) D E 0 = − μ 0 − β * ( 1 − ε m c m ) ( 1 − ε v ) Λ μ − β * ( 1 − ε m c m ) Λ μ 0 0 β 4 β * ( 1 − ε m c m ) ( 1 − ε v ) Λ μ β * ( 1 − ε m c m ) Λ μ 0 0 ν − β 6 0 0 0 θ γ − β 5 0 0 0 0 δ − μ .

At β = β * , it is clear that 0 is a simple eigenvalue of D E 0 . A corresponding right eigenvector associated with 0 is denoted as w = ( w 1 , w 2 , w 3 , w 4 , w 5 ) .

(23) w = − β 4 β 5 β 6 [ ( 1 − ε v ) β 5 + ϕ ( ν 2 + θ β 6 ) ] δ [ ( 1 − ε v ) ν β 5 + ϕ γ ν + ϕ θ β 6 ] [ ν 2 + θ β 6 ] μ β 5 β 6 δ ( ν 2 + θ β 6 ) μ ν β 5 δ [ ν 2 + θ β 6 ] μ δ 1

and a left eigenvector v = ( v 1 , v 2 , v 3 , v 4 , v 5 ) satisfying v . w = 1 is

(24) v = 0 [ γ β * ( 1 − ε m c m ) + β * ( 1 − ε m c m ) ( 1 − ε v ) β 5 + θ β 6 ] v 4 β 4 β 6 [ γ β * ( 1 − ε m c m ) + β * ( 1 − ε m c m ) ( 1 − ε v ) β 5 ] v 4 β 6 δ [ ν 2 + θ β 6 ] β 4 β 6 β 6 ( β 5 β 6 + 1 ) [ γ + ( 1 − ε v ) β 5 ] [ θ β 5 β 6 + δ β 4 ( ν 2 + θ β 6 ) ] 0 .

Now, we calculate the bifurcation parameters

(25) a = ∑ k , i , j = 1 5 v k w i w j ∂ 2 f k ( E 0 , β * ) ∂ x i ∂ x j and b = ∑ k , i = 1 5 v k w i ∂ 2 f k ( E 0 , β * ) ∂ x i ∂ β .

Algebraic calculations show that

(26) ∂ 2 f 2 ∂ x 1 ∂ x 3 = β * ( 1 − ε m c m ) ( 1 − ε v ) , ∂ 2 f 2 ∂ x 1 ∂ x 4 = β * ϕ ( 1 − ε m c m ) , ∂ 2 f 2 ∂ x 3 ∂ x 1 = β * ( 1 − ε m ) ( 1 − ε v ) , ∂ 2 f 2 ∂ x 4 ∂ x 1 = β * ϕ ( 1 − ε m c m ) , ∂ 2 f 3 ∂ x 3 ∂ x 5 = β * η ( 1 − ε m c m ) , ∂ 2 f 2 ∂ x 4 ∂ x 5 = β * η ( 1 − ε m c m ) , ∂ 2 f 2 ∂ x 5 ∂ x 3 = β * η ( 1 − ε m c m ) , ∂ 2 f 2 ∂ x 5 ∂ x 4 = β * η ( 1 − ε m c m ) , ∂ 2 f 2 ∂ x 3 ∂ β = ( 1 − ε m c m ) ( 1 − ε v ) Λ μ , ∂ 2 f 2 ∂ x 4 ∂ β = ( 1 − ε m c m ) Λ μ .

By utilizing both the left and right eigenvectors, along with the mixed partial derivatives that are non-zero, we obtain

(27) a = v 2 w 1 β * ( 1 − ε m c m ) [ w 3 ( 1 − ε v ) + 2 ϕ w 4 ] + v 2 w 3 w 4 β * ( 1 − ε m c m ) + v 3 β * η ( 1 − ε m c m ) [ 2 w 3 + 2 w 4 ]

(28) b = v 2 ( 1 − ε m c m ) ( 1 − ε v ) Λ μ + v 2 ( 1 − ε m c m ) ϕ Λ μ .

From above, it can be readily observed that the value of b is consistently positive, and similarly, the value of a is positive if

(29) β * μ ν ( 1 − ε m c m ) ( 1 − ε v ) β 5 δ [ ν 2 + θ β 6 ] + 2 β * ϕ ( 1 − ε m c m ) μ δ < 0 .

Now, for model (1) to maintain backward bifurcation, both a and b values must be positive at the same time. The backward bifurcation occurs when a stable DFE coexists with a stable endemic equilibrium, leading to sustained transmission of the disease. Understanding the conditions for a backward bifurcation is essential in assessing the stability and long-term behavior of disease dynamics.

5 Sensitivity analysis

In this segment, we will conduct an analysis of the effectiveness of different model parameters on infection spread and prevalence. To conduct our analysis, we numerically simulated the model using the parameters listed in Table 1. Our initial population size was set to S ( 0 ) = 2,00,00,000, E ( 0 ) = 60, I ( 0 ) = 10, H ( 0 ) = 10, and R ( 0 ) = 10 . All numerical simulations were performed using MATLAB.

Table 1

The description of the model parameters

Parameter	Interpretation	Value	Units	Reference
Λ	Inflow rate	1,500	persons day ‒ 1	[2]
β	Transmission rate	0.0000013	persons. day ‒ 1	-do-
μ	Natural death rate	0.0028	day ‒ 1	-do-
c m	Mask compliance	0.1	—	-do-
ε m	Face mask effectiveness	0.5	—	-do-
ε v	Vaccine effectiveness	0.5	—	-do-
ϕ	Modification factor	0.2	day ‒ 1	-do-
δ	Recovery rate	1 14	day ‒ 1	-do-
γ	Rate of infective in care ( I to H class)	0.1	day ‒ 1	-do-
ψ	Mortality rate in H class	0.0245	day ‒ 1	-do-
η	Modification parameter	0.01	day ‒ 1	-do-
ν	Transition rate ( E to I class)	0.0014	day ‒ 1	-do-
θ	Screening rate	0.000167	day ‒ 1	-do-

This study involves conducting a local sensitivity analysis through numerical simulations. The analysis involves observing how changing a single input parameter value affects the output values while keeping all other parameter values constant. Specifically, we investigate the sensitivity of five parameters β , ν , γ , and ε v . It is simple to see in Figure 2(a)–(e) as the transmission rate rises up, for instance, comparing β = 0.000001 and β = 0.000003 curves represents a decrease in susceptible population, and the population is moving to all other classes. When the infection rate ( ν ) increases, the population in E class is moving to class I and the number of infected persons is increasing at its peak in Figure 3(a)–(e). In Figure 4(a)–(e), we can observe that the number of people in class H is increasing due to the infection rate from class I. If we vaccinate and provide care for infected individuals, the graphs in Figure 5(a)–(e) clearly show a diminishing in the number of infections. In Figure 6, the bar plot shows that if we increment β by 1% while holding the remaining variables constant, it results in a 1% rise in ℛ 0 . This pattern is demonstrated in the other factors as well.

$Figure 2 Sensitivity of system (1) with varying transmission rate ( β \beta ) in all classes.$

Figure 2

Sensitivity of system (1) with varying transmission rate ( β ) in all classes.

$Figure 3 Sensitivity of system (1) with varying transition rate of infective class ( ν \nu ) in all classes.$

Figure 3

Sensitivity of system (1) with varying transition rate of infective class ( ν ) in all classes.

$Figure 4 Sensitivity of system (1) with infection rate ( γ \gamma ).$

Figure 4

Sensitivity of system (1) with infection rate ( γ ).

$Figure 5 Sensitivity of system (1) with varying rate of vaccine efficacy ( ε v {\varepsilon }_{v} ) in classes.$

Figure 5

Sensitivity of system (1) with varying rate of vaccine efficacy ( ε v ) in classes.

$Figure 6 In this figure, ℛ 0 β {{\mathcal{ {\mathcal R} }}}_{0}^{\beta } = 1, ℛ 0 ν {{\mathcal{ {\mathcal R} }}}_{0}^{\nu } = 0.2786, ℛ 0 c m = ‒ 0.1052 {{\mathcal{ {\mathcal R} }}}_{0}^{{c}_{m}}=‒0.1052 , ℛ 0 ε m = ‒ 0.5263 {{\mathcal{ {\mathcal R} }}}_{0}^{{\varepsilon }_{m}}=‒0.5263 , and ℛ 0 ε v = ‒ 0.7104 {{\mathcal{ {\mathcal R} }}}_{0}^{{\varepsilon }_{v}}=‒0.7104 .$

Figure 6

In this figure, ℛ 0 β = 1, ℛ 0 ν = 0.2786, ℛ 0 c m = ‒ 0.1052 , ℛ 0 ε m = ‒ 0.5263 , and ℛ 0 ε v = ‒ 0.7104 .

Our analysis reveals the intricate interplay between mask efficiency and mask compliance. The contour plot of ℛ 0 in Figure 7 over a range of mask efficiency and mask compliance values provides valuable insights into the potential effectiveness of these measures in controlling disease spread. We observe that higher vaccination efficiency ( ε v ) and mask compliance ( c m ) lead to lower ℛ 0 values, indicating reduced transmission potential. The contour lines on the plot, labeled with ℛ 0 values, help visualize the regions of varying disease transmission risk.

$Figure 7 Contour plot of ℛ 0 {{\mathcal{ {\mathcal R} }}}_{0} vs mask efficiency and mask compliance.$

Figure 7

Contour plot of ℛ 0 vs mask efficiency and mask compliance.

6 Optimal control analysis

Effective management techniques are essential for controlling the transmission of infectious diseases and lessening their effects on the general populace. Although conventional strategies like immunization and treatment are crucial, other nonpharmaceutical treatments are frequently used to further limit the dynamics of transmission. This study attempts to optimize the implementation of three important techniques, including physical or social distance measures, contact tracing, and quarantine and isolation of symptomatically infected persons with appropriate treatment represented by the control variables.

The SEIHR model serves as the foundation for understanding disease dynamics and encompassing the impacts of the implemented control strategies. The control problem is formulated to determine the optimal time-varying functions u 1 ( t ) , u 2 ( t ) , and u 3 ( t ) that minimize an objective function capturing the overall disease burden or maximizing a desired outcome. The objective function includes factors such as the overall number of exposed individuals, diseased individuals, and economic costs pertained to the outbreak. By optimizing the control strategies, we aim to find the optimal balance between minimizing disease transmission and reducing the negative impacts of the interventions on the population and society.

The control variable u 1 ( t ) represents the practice of physical or social distancing measures, which aim to reduce close contact and interactions among individuals. This control variable affects the transmission rate in the SEIHR model, reflecting the effectiveness of the implemented measures in reducing the spread of the disease.
If u 1 = 0 , it indicates that the disease transmission dynamics are not affected by the control technique u 1 for implementing physical or social distance measures. This suggests that no concrete steps have been taken to stop the spread of the disease through physical or interpersonal isolation. Any u 1 -related control efforts do not impact the disease’s ability to spread freely among the susceptible population.
If u 1 = 1 , it signifies that the disease transmission prevention method of using physical or social distance measures is fully implemented and successful. This suggests that the measures of social or physical isolation are extremely effective in completely reducing the rate of disease transmission. The prevention or decrease of intimate contacts and interactions between people is achieved by the control method executed by u 1 , which significantly slows the spread of the illness.
The control variable u 2 ( t ) captures the impact of contact tracing and quarantine measures. Contact tracing entails the verification of individuals who have met infected individuals and the subsequent isolation of such contacts to hinder the onward transmission of the infection. This control variable influences the transition rates from E to I state, reflecting the effectiveness of contact tracing efforts in detecting and isolating potential cases.
If u 2 = 0 , the contact tracing and quarantine controlling technique are either not being used or have no impact on the dynamics of the disease. The equations stay the same in this situation, and no further precautions are taken to detect and separate those who have come into contact with sick people.
If u 2 = 1 , it means that the contact tracing and quarantine control technique u 2 are completely implemented and successful in identifying and isolating those who have come into contact with infected people. This can aid in stopping the illness from spreading further.
The control variable u 3 ( t ) signifies the isolation of individuals displaying symptoms of infection, along with the administration of suitable treatment. This control variable affects the infected population (I), reflecting the impact of isolating and treating symptomatic cases on the overall disease dynamics.
If u 3 = 0 , it indicates that the control method ( u 3 ) for isolating symptomatically infected persons and treating them appropriately is not used or has no impact on the dynamics of the illness. This suggests that no special precautions have been taken to segregate and treat those with symptomatic infections. Without any focused measures for their isolation and treatment, the infected people (I) continue to contribute to the dynamics of disease transmission.
When u 3 = 1 , it suggests that focused measures are in place to identify, isolate, and treat individuals with symptoms. This can have a positive impact on the dynamics of the illness by limiting the transmission and potentially slowing down the infection.

By solving the optimal control problem, we seek to determine the time-dependent profiles that minimize the objective function while considering the constraints and dynamics of the SEIHR model. The optimal control strategies obtained from this analysis can guide policymakers and public health officials in making informed decisions about the implementation and timing of interventions to effectively manage disease outbreaks.

The optimal control problem is modeled as follows:

(30) d S d t = Λ − β 1 S I ( 1 − u 1 ) − β 2 ( 1 − u 1 ) S H − μ S d E d t = β 1 S I ( 1 − u 1 ) + β 2 S H ( 1 − u 1 ) − ( ν + u 2 ) E − θ E − μ E d I d t = ( ν + u 2 ) E + β 3 R ( I + H ) − ( γ + u 3 ) I − μ I d H d t = θ E − β 5 H + ( γ + u 3 ) I d R d t = δ H − μ R − β 3 R ( I + H )

Objective: To minimize the cost functional:

(31) Z ( u 1 , u 2 , u 3 ) = ∫ t 0 t f a 1 E ( t ) + a 2 I ( t ) + 1 2 ( w 1 u 1 2 + w 2 u 2 2 + w 3 u 3 2 ) d t

subject to system (30).

where a 1 and a 2 are weighting factors that assess the comparative significance of the exposed and infected individuals in the cost function and w 1 , w 2 , and w 3 are costs that determine the relative importance of each control strategy in the cost function, the cost associated with the implemented controls is presumed to follow a non-linear and quadratic relationship.

All control effects, namely u 1 ( t ) , u 2 ( t ) , and u 3 ( t ) , are presumed to be bounded and Lebesgue measurable time-dependent functions over the interval [ 0 , t f ] , where t f represents the final time. Henceforth, the objective lies in attaining optimal control u * = ( u 1 * , u 2 * , u 3 * ) , wherein the corresponding state trajectories serve as the solution to the system represented by the equation (30) for t ∈ [ 0 , t f ] . The primary aim is to minimize the specified cost functional as delineated in equation (31). That is,

(32) Z ( u * ( t ) ) = min ( Z ( u ) : u ∈ U ad ) ,

where Z ( u ) represents the cost functional defined in equation (31) and u belongs to the admissible control set U ad . The admissible control set is as follows:

(33) U ad = { u = ( u 1 , u 2 , u 3 ) : 0 < u i < 1 , t ∈ [ 0 , t f ] }

To find the Hamiltonian for the optimal control problem, we need to introduce the costate variables associated with each state variable in the SEIHR model. Let’s denote the costate variables as P S , P E , P I , P H , and P R corresponding to the state variables S , E , I , H , and R , respectively.

The PMP helps to reduce problems (30), (31), and (33) to a problem of minimizing the Hamiltonian. The Hamiltonian, denoted by ℋ is defined as the sum of the instantaneous cost and the inner product of the costate variables with the time derivatives of the state variables. In this case, the Hamiltonian for the SEIHR model with optimal control can be written as follows:

(34) ℋ = a 1 E + a 2 I + 1 2 ( w 1 u 1 2 + w 2 u 2 2 + w 3 u 3 2 ) + P S ( Λ − β 1 S I ( 1 − u 1 ) − β 2 ( 1 − u 1 ) S H − μ S ) + P E ( β 1 S I ( 1 − u 1 ) + β 2 S H ( 1 − u 1 ) − ( ν + u 2 ) E − θ E − μ E ) + P I ( ( ν + u 2 ) E + β 3 R ( I + H ) − ( γ + u 3 ) I − μ I ) + P H ( θ E − β 5 H + ( γ + u 3 ) I ) + P R ( δ H − μ R − β 3 R ( I + H ) )

Theorem 4

Let u * be the optimal solution to equations (30), and (31), (33), and ( S * , E * , I * , H * , R * ) be the associated optimal state. Then, there exist adjoint functions P S , P E , P I , P H , and P R that satisfy the adjoint system.

(35) d P S d t = μ P S + ( P S − P E ) ( β 1 I ( 1 − u 1 ) + β 2 H ( 1 − u 1 ) ) d P E d t = − a 1 + μ P E + ( P E − P I ) ( ν + u 2 ) + ( P E − P H ) θ d P I d t = − a 2 + μ P I + ( P S − P E ) β 1 S ( 1 − u 1 ) + ( P R − P I ) β 3 R + ( P I − P H ) ( γ + u 3 ) d P H d t = μ P H + β 2 S ( 1 − u 1 ) ( P S − P E ) + ( P H − P R ) β 3 R ( I + H ) + ( P H − P R ) δ + ψ P H d P R d t = μ P R + ( P R − P I ) β 3 ( I + H ) ,

subject to final time conditions P S ( t f ) = 0 , P E ( t f ) = 0 , and P I ( t f ) = 0 , P H ( t f ) = 0 , P R ( t f ) = 0 .

Furthermore, the optimal control, denoted as u * = ( u 1 * , u 2 * , u 3 * ) is expressed as follows:

(36) u 1 * = min max 0 , ( β 1 S I + β 2 S H ) ( P E − P S ) w 1 , 1 u 2 * = min max 0 , ( P E − P I ) E w 2 , 1 u 3 * = min max 0 , ( P I − P H ) I w 3 , 1 .

Proof

This result can be established through PMP, wherein:

(37) d P S d t = − ∂ ℋ ∂ S , d P E d t = − ∂ ℋ ∂ E , d P I d t = − ∂ ℋ ∂ I , d P H d t = − ∂ ℋ ∂ H , d P R d t = − ∂ ℋ ∂ R .

with final time conditions P S ( t f ) = 0 , P E ( t f ) = 0 , P I ( t f ) = 0 , P H ( t f ) = 0 , and P R ( t f ) = 0 .

The characterization of controls is derived by solving u 1 * and u 2 * from the following:

(38) ∂ ℋ ∂ u 1 = 0 , ∂ ℋ ∂ u 2 = 0 , ∂ ℋ ∂ u 3 = 0 ,

which implies that

(39) u 1 * = ( β 1 S I + β 2 S H ) ( P E − P S ) w 1 u 2 * = ( P E − P I ) E w 2 u 3 * = ( P I − P H ) I w 3

By virtue of the boundedness of u i * ( t ) on the interval (0, 1) and the condition of minimality, it follows that

(40) u 1 * ( t ) = 0 , if ∂ H ∂ u 1 > 0 ( β 1 S I + β 2 S H ) ( P E − P S ) w 1 if ∂ H ∂ u 1 = 0 1 , if ∂ H ∂ u 1 < 0 u 2 * ( t ) = 0 , if ∂ H ∂ u 2 > 0 ( P E − P I ) E w 2 if ∂ H ∂ u 2 = 0 1 , if ∂ H ∂ u 2 < 0 u 3 * ( t ) = 0 , if ∂ H ∂ u 3 > 0 ( P I − P H ) I w 3 if ∂ H ∂ u 3 = 0 1 , if ∂ H ∂ u 3 < 0 .

Hence, the proof is complete.□

6.1 Simulation of the control problem

To evaluate the efficacy of the recommended preventive measures, numerical simulations are performed on the proposed COVID-19 model, both with and without control measures. The widely utilized Runge-Kutta fourth order (RK-4) iterative scheme is employed to calculate the numerical solutions for the model.

The simulation utilizes the parameter values as presented in Table 1. The weights and balancing constants have been chosen in the following manner: a 1 = 8 , a 2 = 1 , w 1 = 4,500, w 2 = 10,000, and w 3 = 1,000. These values are chosen to appropriately weigh the different factors and components within the model. The dynamics pertaining to discrete populations, both with and without control measures, are presented in Figures 8–15. The simulations with control measures, are represented by dashed curves, while the simulations without control measures are depicted with solid curves.

$Figure 8 Graphical results of model (1) with practicing social distancing protocols, u 1 ≠ 0 , u 2 = 0 {u}_{1}\ne 0,{u}_{2}=0 , and u 3 = 0 {u}_{3}=0 .$

Figure 8

Graphical results of model (1) with practicing social distancing protocols, u 1 ≠ 0 , u 2 = 0 , and u 3 = 0 .

To evaluate each control strategy’s efficacy and influence on disease incidence, the control model considers three different scenarios. These situations are designed based on multiple control intervention strategies, namely single, coupled, and threefold control variables.

In the initial case, a one-control strategy is employed at a time. This scenario investigates the impact of individual control measures when utilized in isolation. The second scenario, referred to as coupled controls, examines the combined effects of two control measures simultaneously. By considering two controls together, the study evaluates their joint impact on the dynamics of the model state variables.

Ultimately, the final scenario evaluates the ramifications of all implemented control measures, where all control strategies are simultaneously implemented. This scenario provides insights into the collective influence of multiple regulatory measures influencing the model’s dynamics. A detailed analysis is provided for the graphical impact of each scenario, highlighting the observed changes in the state variables.

6.1.1 Scenario 1: Strategies with one control variable

In this situation, three control strategies u 1 (practicing physical or social distancing protocols), u 2 (effective contact tracing and quarantine), and u 3 (isolation of symptomatically infected individuals with treatment) are examined in separate discussions to delve into the distinct influence of each approach on the disease dynamics, we conduct an investigation. In the initial scenario, we examine the control set where u 1 ≠ 0 , u 2 = 0 , and u 3 = 0 . In simpler terms, only the time-dependent isolation control variable is applied, while the other two variables are not considered.

The above control strategy aligns with the implementation of physical distancing protocols. The comprehensive details can be observed in Figure 8(a)–(e). Notably, with this control strategy, there is a reduction in the count of exposed, infected, and recovered individuals. This suggests that reinforcing physical distancing measures aids in reducing interpersonal contact within the community, subsequently minimizing the susceptibility to infection, such as by avoiding public gatherings. On the other hand, Figure 9(a)–(e) portrays the impact of a single-control strategy, assuming u 2 ≠ 0 . However, this single strategy does not contribute significantly to reducing the number of infected and infective individuals receiving care.

$Figure 9 Graphical results of model (1) with effective contact tracing and quarantine, u 1 = 0 , u 2 ≠ 0 {u}_{1}=0,{u}_{2}\ne 0 , and u 3 = 0 {u}_{3}=0 .$

Figure 9

Graphical results of model (1) with effective contact tracing and quarantine, u 1 = 0 , u 2 ≠ 0 , and u 3 = 0 .

Comparatively, Figure 10(a)–(e) showcases the influence of symptomatically infected individuals receiving treatment, demonstrating a significant decrease in population in all classes. Among all the strategies examined, the implementation of u 1 ≠ 0 , u 2 = 0 , and u 3 = 0 prove to be the most effective in reducing the number of infections within the population.

$Figure 10 Graphical results of model (1) with isolation of infected with treatment, u 1 = 0 , u 2 = 0 {u}_{1}=0,{u}_{2}=0 , and u 3 ≠ 0 . {u}_{3}\ne 0.$

Figure 10

Graphical results of model (1) with isolation of infected with treatment, u 1 = 0 , u 2 = 0 , and u 3 ≠ 0 .

The impact of physical distancing ( u 1 ) reveals a reduction in exposed, infected, and recovered individuals, indicating its effectiveness. In contrast, contact tracing and quarantine ( u 2 ) shows limited efficacy in reducing infected individuals. Isolating and treating symptomatic cases ( u 3 ) leads to a significant decrease in the overall population across all classes. The most effective strategy observed is a combination of u 1 ≠ 0 , u 2 = 0 , and u 3 = 0 , emphasizing the importance of implementing physical distancing measures for comprehensive disease control.

6.1.2 Scenario 2: Strategies with coupled-control variables

The second scenario involves conducting an analysis of the utilization of coupled-control strategies to examine their impact on the dynamics of COVID-19. The efficient contact tracing, quarantine of individuals, and medical treatment to the symptomatically infected individuals are used as control strategies in the beginning, i.e., u 1 = 0 , u 2 ≠ 0 , and u 3 ≠ 0 , while social distancing is not taken into account.

The biological impact is shown in Figure 11(a)–(e), the strict implementation of quarantine and treatment for symptomatic will decrease the exposed population and a large amount of the population will be recovered. Next, we are implementing social distancing protocols and treatment for symptomatically infected people leads to a reduction in the number of infections in Figure 12(a)–(e). The physical or social distancing measure and effective contact tracing and quarantine in Figure 13(a)–(e) will have the same effect as the previous one.

$Figure 11 Graphical results of model (1) with effective contact tracing and quarantine and isolation of infected with treatment, u 1 = 0 , u 2 ≠ 0 {u}_{1}=0,{u}_{2}\ne 0 , and u 3 ≠ 0 . {u}_{3}\ne 0.$

Figure 11

Graphical results of model (1) with effective contact tracing and quarantine and isolation of infected with treatment, u 1 = 0 , u 2 ≠ 0 , and u 3 ≠ 0 .

$Figure 12 Graphical results of model (1) with practicing social distances and isolation of infected with treatment, u 1 ≠ 0 , u 2 = 0 {u}_{1}\ne 0,{u}_{2}=0 , and u 3 ≠ 0 . {u}_{3}\ne 0.$

Figure 12

Graphical results of model (1) with practicing social distances and isolation of infected with treatment, u 1 ≠ 0 , u 2 = 0 , and u 3 ≠ 0 .

$Figure 13 Graphical results of the model (1) with practicing social distances and effective contact tracing and quarantine, u 1 ≠ 0 , u 2 ≠ 0 {u}_{1}\ne 0,{u}_{2}\ne 0 , and u 3 = 0 . {u}_{3}=0.$

Figure 13

Graphical results of the model (1) with practicing social distances and effective contact tracing and quarantine, u 1 ≠ 0 , u 2 ≠ 0 , and u 3 = 0 .

In this scenario, the coupled controls of contact tracing, quarantine ( u 2 ), and treatment ( u 3 ) exhibit a biologically meaningful reduction in the exposed population. This suggests that targeted interventions, such as isolating and treating infected individuals, can effectively break the chain of transmission, leading to a decline in the number of individuals at risk of infection and promoting a robust recovery phase.

6.1.3 Scenario 3: Strategies with three control variables

In the third scenario, we conducted an analysis to explore the combined impact of social distancing protocols, effective contact tracing, quarantine of individuals, and treatment for symptomatically infected individuals on the dynamics of COVID-19 infections. Our objective was to examine the consequences of the simultaneous implementation of all control measures on disease incidence. With this objective in mind, we considered the implementations u 1 ≠ 0 , u 2 ≠ 0 , and u 3 ≠ 0 and simulated the model with all possible cases. The dynamics observed within the respective population compartments, both in the presence and absence of control measures, are illustrated in Figure 14(a)–(e). From Figure 15(a)–(c), control profile graphs give specific interventions or control strategies to influence the spread and impact of the disease over time. It is evident that with the implementation of all controls, the susceptible population decreases rapidly, while the vaccinated population experiences significant growth compared to the scenario without any control strategy or with constant control. Furthermore, the exposed and infected populations also exhibit a significant decline. At last, the simulation in this scenario indicates that the simultaneous implementation of the proposed control measures is highly effective and essential in minimizing the spread of infection within a community, thereby safeguarding against future disease incidence.

$Figure 14 Graphical results of model (1) with practicing social distances, effective contact tracing, and quarantine and isolating infected with treatment, u 1 ≠ 0 , u 2 ≠ 0 {u}_{1}\ne 0,{u}_{2}\ne 0 , and u 3 ≠ 0 . {u}_{3}\ne 0.$

Figure 14

Graphical results of model (1) with practicing social distances, effective contact tracing, and quarantine and isolating infected with treatment, u 1 ≠ 0 , u 2 ≠ 0 , and u 3 ≠ 0 .

Figure 15

(a) Impact of Physical Distancing (u ₁); (b) Impact of effective contact tracing and quarantine (u ₂); and (c) Impact of Isolation of symptomatically infected individuals with treatment (u ₃).

In this scenario 3, the simultaneous implementation of all controls manifests a biologically favorable outcome. The rapid decrease in susceptibility indicates the combined measures’ effectiveness in preventing new infections, while the significant increase in vaccinations reflects an enhanced immune response within the population. The substantial decline in exposed and infected populations further emphasizes the comprehensive control measures ability to mitigate the overall impact of the virus on public health.

7 Conclusions

In conclusion, this article comprehensively examines an SEIHR compartmental model to dissect various facets of COVID-19 dynamics.

For the long-term behavior of the disease, we analyze the global stability analysis to know the overall trends of the pandemic.
We investigate the impact of different model parameters on the spread and prevalence of infection.
We show the existence of backward bifurcation to challenges the assumption even if ℛ 0 < 1 , disease can persist in the population.
In optimal control analysis, a control profile graph represents the behavior or values of the control variables over a given time period.
It provides a visual representation of how the control variables change over time to achieve the desired objectives.

Madubueze et al.’s [19] work on the optimal control analysis of the SEQIJR model for COVID-19, incorporating public health education, population quarantine, and isolation of infectious individuals is commendable. While their study encompasses crucial measures, it lacks the explicit consideration of physical or social distancing and detailed control variables like u 1 ( t ) , u 2 ( t ) , and u 3 ( t ) . In contrast, our research offers a granular approach, delineating specific control strategies for physical distancing, contact tracing, and symptomatic isolation with treatment. This finer analysis in our work provides a comprehensive understanding of various intervention measures, enhancing the applicability and insights derived from the model.

In conclusion, the analysis of different control strategies for COVID-19 dynamics reveals that implementing physical or social distancing protocols, effective contact tracing and quarantine, and treatment for symptomatically infected individuals all contribute to reducing infections. Among the individual control strategies, physical or social distancing measures show the most significant impact. However, the combined implementation of all three control measures simultaneously yields the most favorable results, leading to a rapid decrease in susceptible individuals, significant recovery, and enhanced protection against future disease incidence.

Acknowledgements

The authors express gratitude to the editor and anonymous reviewers for their valuable suggestions, which have significantly enhanced the manuscript’s presentation.

Funding information: This research received no specific grant from any funding agency, commercial, or nonprofit sectors.
Conflict of interest: The authors have no conflicts of interest to disclose.
Ethical approval: This research did not require ethical approval.
Data availability statement: This article does not involve data sharing, as no datasets were generated or analyzed during the current study.

Appendix

The coefficients of the characteristic polynomial in equation (12) are given as

D 1 = − ( c 11 + c 22 + c 33 + c 44 + c 55 ) , D 2 = c 11 c 22 − c 23 c 32 + c 11 c 33 + c 22 c 33 − c 24 c 42 − c 34 c 43 + c 11 c 44 + c 22 c 44 + c 33 c 44 − c 35 c 53 + c 11 c 55 + c 22 c 55 + c 33 c 55 + c 44 c 55 , D 3 = − c 13 c 21 c 32 + c 11 c 23 c 32 − c 11 c 22 c 33 − c 14 c 21 c 42 + c 11 c 24 c 42 + c 24 c 33 c 42 − c 23 c 34 c 42 − c 24 c 32 c 43 + c 11 c 34 c 43 + c 22 c 34 c 43 − c 11 c 22 c 44 + c 23 c 32 c 44 − c 11 c 33 c 44 − c 22 c 33 c 44 + c 11 c 35 c 53 + c 22 c 35 c 53 + c 35 c 44 c 53 − c 35 c 43 c 54 − c 11 c 22 c 55 + c 23 c 32 c 55 − c 11 c 33 c 55 − c 22 c 33 c 55 + c 24 c 42 c 55 + c 34 c 43 c 55 − c 11 c 44 c 55 − c 22 c 44 c 55 − c 33 c 44 c 55 , D 4 = c 14 c 21 c 33 c 42 − c 11 c 24 c 33 c 42 − c 13 c 21 c 34 c 42 + c 11 c 23 c 34 c 42 − c 14 c 21 c 32 c 43 + c 11 c 24 c 32 c 43 − c 11 c 22 c 34 c 43 + c 13 c 21 c 32 c 44 − c 11 c 23 c 32 c 44 + c 11 c 22 c 33 c 44 − c 11 c 22 c 35 c 53 + c 24 c 35 c 42 c 53 − c 11 c 35 c 44 c 53 − c 22 c 35 c 44 c 53 − c 23 c 35 c 42 c 54 + c 11 c 35 c 43 c 54 + c 22 c 35 c 43 c 54 + c 13 c 21 c 32 c 55 − c 11 c 23 c 32 c 55 + c 11 c 22 c 33 c 55 + c 14 c 21 c 42 c 55 − c 11 c 24 c 42 c 55 − c 24 c 33 c 42 c 55 + c 23 c 34 c 42 c 55 + c 24 c 32 c 43 c 55 − c 11 c 34 c 43 c 55 − c 22 c 34 c 43 c 55 + c 11 c 22 c 44 c 55 − c 23 c 32 c 44 c 55 + c 11 c 33 c 44 c 55 + c 22 c 33 c 44 c 55 , D 5 = − c 14 c 21 c 35 c 42 c 53 − c 11 c 24 c 35 c 42 c 53 + c 11 c 22 c 35 c 44 c 53 − c 13 c 21 c 35 c 42 c 54 + c 11 c 23 c 35 c 42 c 54 − c 11 c 22 c 35 c 43 c 54 − c 14 c 21 c 33 c 42 c 55 + c 11 c 24 c 33 c 42 c 55 + c 13 c 21 c 34 c 42 c 55 − c 11 c 23 c 34 c 42 c 55 + c 14 c 21 c 32 c 43 c 55 − c 11 c 24 c 32 c 43 c 55 + c 11 c 22 c 34 c 43 c 55 − c 13 c 21 c 32 c 44 c 55 + c 11 c 23 c 32 c 44 c 55 − c 11 c 22 c 33 c 44 c 55 .

References

[1] Arcede, J. P., Caga-Anan, R. L., Mentuda, C. Q., & Mammeri, Y. (2020). Accounting for symptomatic and asymptomatic in a SEIR-type model of COVID-19. Mathematical Modelling of Natural Phenomena, 15, 34. 10.1051/mmnp/2020021Search in Google Scholar

[2] Bandekar, S. R., Das, T., Srivastav, A. K., Yadav, A., Kumar, A., Srivastava, P. K., & Ghosh, M. (2022). Modeling and prediction of the third wave of COVID-19 spread in India. Computational and Mathematical Biophysics, 10(1), 231–248. 10.1515/cmb-2022-0138Search in Google Scholar

[3] Bandekar, S. R., Ghosh, M., & Bi, K. (2023). Impact of high-risk and low-risk population on COVID-19 dynamics considering antimicrobial resistance and control strategies. The European Physical Journal Plus, 138(8), 697. 10.1140/epjp/s13360-023-04328-zSearch in Google Scholar

[4] Brauer, F., Castillo-Chavez, C., & Castillo-Chavez, C. (2012). Mathematical models in population biology and epidemiology (Vol. 2, No. 40). New York: Springer. 10.1007/978-1-4614-1686-9Search in Google Scholar

[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] Cooper, I., Mondal, A., & Antonopoulos, C. G. (2020). A SIR model assumption for the spread of COVID-19 in different communities. Chaos, Solitons & Fractals, 139, 110057. 10.1016/j.chaos.2020.110057Search in Google Scholar PubMed PubMed Central

[7] Das, T., Bandekar, S. R., Srivastav, A. K., Srivastava, P. K., & Ghosh, M. (2023). Role of immigration and emigration on the spread of COVID-19 in a multipatch environment: a case study of India. Scientific Reports, 13(1), 10546. 10.1038/s41598-023-37192-zSearch in Google Scholar PubMed PubMed Central

[8] Daud, A. A. M. (2021). A note on Lienard-Chipart criteria and its application to epidemic models. Mathematics and Statistics, 9(1), 41–45. 10.13189/ms.2021.090107Search in Google Scholar

[9] Diekmann, O., Heesterbeek, J. A. P., & Metz, J. A. (1990). On the definition and the computation of the basic reproduction ratio R 0 in models for infectious diseases in heterogeneous populations. Journal of Mathematical Biology, 28, 365–382. 10.1007/BF00178324Search in Google Scholar PubMed

[10] Kamrujjaman, M., Saha, P., Islam, M. S., & Ghosh, U. (2022). Dynamics of SEIR model: A case study of COVID-19 in Italy. Results in Control and Optimization, 7, 100119. 10.1016/j.rico.2022.100119Search in Google Scholar

[11] Kapur, J. N. (2008). Mathematical models in biology and medicine. East-West Press Private Limited. Search in Google Scholar

[12] Kaur, S., Bherwani, H., Gulia, S., Vijay, R., & Kumar, R. (2021). Understanding COVID-19 transmission, health impacts and mitigation: timely social distancing is the key. Environment, Development and Sustainability, 23, 6681–6697. 10.1007/s10668-020-00884-xSearch in Google Scholar PubMed PubMed Central

[13] Khan, A. A., Ullah, S., & Amin, R. (2022). Optimal control analysis of COVID-19 vaccine epidemic model: a case study. The European Physical Journal Plus, 137(1), 1–25. 10.1140/epjp/s13360-022-02365-8Search in Google Scholar PubMed PubMed Central

[14] Koutou, O., & Sangaré, B. (2023). Mathematical analysis of the impact of the media coverage in mitigating the outbreak of COVID-19. Mathematics and Computers in Simulation, 205, 600–618. 10.1016/j.matcom.2022.10.017Search in Google Scholar PubMed PubMed Central

[15] Kumar, A., & Srivastava, P. K. (2023). Role of optimal screening and treatment on infectious diseases dynamics in presence of self-protection of susceptible. Differential Equations and Dynamical Systems, 31(1), 135–163. 10.1007/s12591-019-00467-xSearch in Google Scholar

[16] Kumari, P., Singh, H. P., & Singh, S. (2021). SEIAQRDT model for the spread of novel coronavirus (COVID-19): A case study in India. Applied Intelligence, 51, 2818–2837. 10.1007/s10489-020-01929-4Search in Google Scholar PubMed PubMed Central

[17] Lamba, S. & Srivastava, P. (2023). Cost-effective optimal control analysis of a COVID-19 transmission model incorporating community awareness and waning immunity. Computational and Mathematical Biophysics, 11(1), 20230154. 10.1515/cmb-2023-0154Search in Google Scholar

[18] Li, Z., & Zhang, T. (2022). Analysis of a COVID-19 epidemic model with seasonality. Bulletin of Mathematical Biology, 84(12), 146. 10.1007/s11538-022-01105-4Search in Google Scholar PubMed PubMed Central

[19] Madubueze, C. E., Dachollom, S., & Onwubuya, I. O. (2020). Controlling the spread of COVID-19: optimal control analysis. Computational and Mathematical methods in Medicine, 2020, Article ID 6862516. 10.1101/2020.06.08.20125393Search in Google Scholar

[20] Martcheva, M. (2015). An introduction to mathematical epidemiology (Vol. 61, pp. 9–31). New York: Springer. 10.1007/978-1-4899-7612-3_2Search in Google Scholar

[21] Nayak, D., Chhetri, B., Vamsi Dasu, K., Muthusamy, S. & Bhagat, V. (2023). A comprehensive and detailed within-host modeling study involving crucial biomarkers and optimal drug regimen for type I Lepra reaction: A deterministic approach. Computational and Mathematical Biophysics, 11(1), 20220148. 10.1515/cmb-2022-0148Search in Google Scholar

[22] Paul, S., Mahata, A., Ghosh, U., & Roy, B. (2021). Study of SEIR epidemic model and scenario analysis of COVID-19 pandemic. Ecological Genetics and Genomics, 19, 100087. 10.1016/j.egg.2021.100087Search in Google Scholar PubMed PubMed Central

[23] Perko, L. (2013). Differential equations and dynamical systems (Vol. 7). Springer Science & Business Media. Search in Google Scholar

[24] Peter, O. J., Qureshi, S., Yusuf, A., Al-Shomrani, M., & Idowu, A. A. (2021). A new mathematical model of COVID-19 using real data from Pakistan. Results in Physics, 24, 104098. 10.1016/j.rinp.2021.104098Search in Google Scholar PubMed PubMed Central

[25] Prakash, D. B., Chhetri, B., Vamsi, D. K. K., Balasubramanian, S., & Sanjeevi, C. B. (2021). Low temperatures or high isolation delay increases the average COVID-19 infections in India: A Mathematical modeling approach. Computational and Mathematical Biophysics, 9(1), 146–174. 10.1515/cmb-2020-0122Search in Google Scholar

[26] Prete, G., Fontanesi, L., Porcelli, P., & Tommasi, L. (2020). The psychological impact of COVID-19 in Italy: worry leads to protective behavior, but at the cost of anxiety. Frontiers in Psychology, 11, 566659. 10.3389/fpsyg.2020.566659Search in Google Scholar PubMed PubMed Central

[27] Rana, K. & Kumari, N. (2023). Application of dynamic mode decomposition and compatible window-wise dynamic mode decomposition in deciphering COVID-19 dynamics of India. Computational and Mathematical Biophysics, 11(1), 20220152. 10.1515/cmb-2022-0152Search in Google Scholar

[28] Sharma, S., & Sharma, P. K. (2021). A study of SIQR model with Holling type-II incidence rate. Malaya Journal of Matematik, 9(1), 305–311. 10.26637/MJM0901/0052Search in Google Scholar

[29] Sharma, S., & Sharma, P.K. (2023). Stability analysis of an SIR model with alert class modified saturated incidence rate and Holling functional type-II treatment. Computational and Mathematical Biophysics, 11. 10.1515/cmb-2022-0145Search in Google Scholar

[30] Sharma, S., & Singh, F. (2021). Bifurcation and stability analysis of a cholera model with vaccination and saturated treatment. Chaos, Solitons & Fractals, 146, 110912. 10.1016/j.chaos.2021.110912Search in Google Scholar

[31] Srivastav, A. K., & Ghosh, M. (2023). Optimal control analysis of the malaria model with saturated treatment. TWMS Journal of Applied & Engineering Mathematics, 13(1), 265–275. Search in Google Scholar

[32] 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

[33] Umdekar, S., Sharma, P. & Sharma, S. (2023). An SEIR model with modified saturated incidence rate and Holling type II treatment function. Computational and Mathematical Biophysics, 11(1), 20220146. 10.1515/cmb-2022-0146Search in Google Scholar

[34] Youssef, H. M., Alghamdi, N. A., Ezzat, M. A., El-Bary, A. A., & Shawky, A. M. (2020). A new dynamical modeling SEIR with global analysis applied to the real data of spreading COVID-19 in Saudi Arabia. Mathematical Biosciences and Engineering, 17(6), 7018–7044. 10.3934/mbe.2020362Search in Google Scholar PubMed

[35] Zeb, A., Alzahrani, E., Erturk, V. S., & Zaman, G. (2020). Mathematical model for coronavirus disease 2019 (COVID-19) containing isolation class. BioMed Research International, 2020, Article ID 3452402. 10.1155/2020/3452402Search in Google Scholar PubMed PubMed Central

[36] URL: https://www.cdc.gov/Search in Google Scholar

[37] URL: https://www.mohfw.gov.in/. Search in Google Scholar

[38] URL: https://covid19.who.int/. Search in Google Scholar

[39] World Health Organization. (2022). A global analysis of COVID-19 intra-action reviews: Reflecting on, adjusting and improving country emergency preparedness and response during a pandemic. Search in Google Scholar

Received: 2023-09-26

Revised: 2023-11-30

Accepted: 2023-12-07

Published Online: 2024-03-09

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-0113

Keywords for this article

global stability; backward bifurcation; optimal control; Lyapunov’s direct method

Creative Commons

BY 4.0