Transmission of COVID-19 in the presence of single-dose and double-dose vaccines with hesitancy: mathematical modeling and optimal control analysis

Introduction: The unexpected emergence of novel coronavirus identified as SAR-CoV-2 virus (severe acute respiratory syndrome corona virus 2) disrupted the world order to an extent that the human activities that are core to survival came almost to a halt. The COVID-19 pandemic created an insurmountable global health crisis that led to a united front among all nations to research on effective pharmaceutical measures that could stop COVID-19 proliferation. Consequently, different types of vaccines were discovered (single-dose and double-dose vaccines). However, the speed at which these vaccines were developed and approved to be administered created other challenges (vaccine skepticism and hesitancy).

Method: This paper therefore tracks the transmission dynamics of COVID-19 using a non-linear deterministic system that accounts for the unwillingness of both susceptible and partially vaccinated individuals to receive either single-dose or double-dose vaccines (vaccine hesitancy). Further the model is extended to incorporate three time-dependent non-pharmaceutical and pharmaceutical intervention controls, namely preventive control, control associated with screening-management of both truly asymptomatic and symptomatic infectious individuals and control associated with vaccination of susceptible individuals with a single dose vaccine. The Pontryagin's Maximum Principle is applied to establish the optimality conditions associated with the optimal controls.

Results: If COVID-19 vaccines administered are imperfect and transient then there exist a parameter space where backward bifurcation occurs. Time profile projections depict that in a setting where vaccine hesitancy is present, administering single dose vaccines leads to a significant reduction of COVID-19 prevalence than when double dose vaccines are administered. Comparison of the impact of vaccine hesitancy against either single dose or double dose on COVID-19 prevalence reveals that vaccine hesitancy against single dose is more detrimental than vaccine hesitancy against a double dose vaccine. Optimal analysis results reveal that non-pharmaceutical time-dependent control significantly flattens the COVID-19 epidemic curve when compared with pharmaceutical controls. Cost-effectiveness assessment suggest that non-pharmaceutical control is the most cost-effective COVID-19 mitigation strategy that should be implemented in a setting where resources are limited.

Discussion: Policy makers and medical practitioners should assess the level of COVID-19 vaccine hesitancy inorder to decide on the type of vaccine (single-dose or double-dose) to administer to the population.

1 Introduction

The advent of COVID-19 infection whose aetiological agent is SARS-CoV-2 virus (severe acute respiratory syndrome corona virus 2) caught the world off guard. The index case of COVID-19 was first confirmed in Wuhan city, Hubei province in China as early as December 2019. Within a period of about 60 days, COVID-19 spread from the epicenter (Wuhan city) to other countries in the world, leading to WHO declaring COVID-19 a global pandemic on March 2020 [1, 2]. The main routes of COVID-19 transmission are direct spread, contact spread, and aerosol spread [3]. Direct spread of COVID-19 occurs as a result of exhaling respiratory droplets from sneezing or coughing. Individuals who are in close contact with other individuals who are suspected or confirmed being infected with COVID-19 can as well be infected [4, 5]. As for the indirect transmission, it may occur when individuals touch eyes, nose, or mouth after contacting contaminated surfaces or objects [4, 5]. Infected individuals are broadly categorized as either asymptomatic (which include both pre-symptomatic and true asymptomatic) or symptomatic. Symptomatic individuals manifest various symptoms which include viral pneumonia, dry cough, fever, myalgia, sore throat, fatigue, nasal congestion, and diarrhea [6, 7]. There is concrete evidence in the sequel that shows that pre-symptomatic transmission significantly contributes to the spread of COVID-19 [8]. Although the field of corona virology was well documented in the sequel [9], the interaction of corona viruses strains with human population on a global scale was not witnessed before. Even when there was SARS outbreak between 2002 and 2003 in North America, South America, Europe, and Asia, it caused about 8,098 cases and 774 SARS-induced fatalities [9] which was significantly low in comparison with 524,467,084 COVID-19-confirmed cases and over 6,915,286 (and increasing) COVID-19-related deaths reported globally [10].

The rapid speed at which COVID-19 established itself within communities triggered myriad challenges such as disruption of normal life, loss of source of livelihood, shrinking of economies (in particular in developing countries), and insurmountable straining of the existing health infrastructure [11]. Exponential increase of critically ill patients in need of urgent medical attention particularly requiring ventilators created a burden within the health system that has never been witnessed in the 21st century. Governments across the world implemented stringent countermeasures to halt the spread of the COVID-19 pandemic. Initially, these suppressive intervention measures were non-pharmaceutical in nature, and they included lock down, contact tracing and quarantine, a ban on air traffic, social distancing, wearing face masks and washing hands regularly, suspending traveling in and out of countries, banning all social gatherings, closing schools, airing awareness programs, and PCR testing for case detection [5, 12]. However, these non-pharmaceutical interventions (NPIs) were short-term measures against proliferation of COVID-19 [13]. A sudden rebound of COVID-19 spread occurred once these intervention measures were relaxed [14]. Consequently, there was a global consensus that pharmaceutical measures such as development of COVID-19 vaccines and cure could mitigate the spread of COVID-19 pandemic.

The proliferation of COVID-19 pandemic compelled multilateral global organizations and private funders to collectively raise funds for researching on an effective vaccine against the novel SARS-CoV-2 [15]. On 11 January 2020, the genetic sequence of SARS-CoV-2 (the causative agent of COVID-19) was made public [16]. Following this breakthrough, the Coalition for Epidemic Preparedness Innovation (CEPI) tasked the global pharmaceutical companies and virology research institutes to conduct a cutting edge research on COVID-19 vaccines. Within a period of about 4 months after emergence of COVID-19, there were about 115 COVID-19 vaccine candidates that were either at exploratory or pre-clinical stages [16]. Shortly afterward, some of these vaccines in particular mRNA-1273 from Moderna, Ad5-nCoV from CanSino Biologicals, INO-4800 from Inovio, and LV-SMENP-DC and pathogen-specific aAPC from Shenzhen Geno-Immune Medical Institute moved into clinical development, and some developed countries started administering them to the general populace [16]. By early 2021, WHO and Food and Drug Administration (FDA) approved more COVID-19 vaccines such as AstraZeneca, Johnson and Johnson, and Pfizer, for the management of COVID-19 pandemic.

It is imperative to note that most of these vaccine innovation breakthroughs and manufacturing were initially based in developed countries. Given most of these high-income countries (HICs) were experiencing an exponential increase of COVID-19 cases on a daily basis, they implemented vaccine nationalism. According to Riaz et al. [15], vaccine nationalism is referred to as an economic strategy that involve hoarding vaccines from manufacturers so that the respective countries can have enough vaccines to vaccinate a large proportion of their populace without being concerned with the limited vaccine production and distribution for the rest of the world. Vaccine nationalism is well known to cause devastating consequences on affected countries. These include leaving many disadvantaged countries with low economic status unable to pay for vaccines as prices of vaccines hike over time; it decreases the likelihood of developing countries to access vaccines as they are likely to be out of stock even when they accumulate enough resources to buy them. Moreover, vaccine nationalism prolongs the duration of the pandemic as it will continue spreading within developing countries unabated [15]. One of the key factors that underpinned vaccine nationalism was that many pharmaceutical companies in various countries were competing in vaccine production. However, the competition was not for the sole purpose of supplying the vaccine to the “ailing humanity” at a subsidized rate but with the aim of selling the vaccines to countries which could offer “strong financial stock upholds” [17]. This policy made it harder for the low- and middle-income countries (LMICs) to afford vaccines, hence alienating them from the global target of vaccination against the lethal SARS-CoV-2. Nevertheless, these pharmaceuticals failed to recognize that the circulation of COVID-19 in developing countries still posed a potential threat to HICs as the cross-border traffic exacerbated COVID-19 transmission. Moreover, many COVID-19 experts in HICs argued that vaccinating some countries while neglecting others burdened with COVID-19 created a huge gap in attaining global immunity [18]. Furthermore, a surge in the population of infected individuals could render an advantageous environment for the mutations of the novel corona virus [19]. Hence, vaccine nationalism was regarded as a self-defeating policy [15].

WHO criticized vaccine nationalism and called for a global funding to increase vaccine production at unprecedented speed. In addition, WHO stressed on the importance of equitable distribution of vaccines among the nations (i.e., both HICs and LMICs) at a subsidized rate. This was in agreement with WHO's slogan that “no one is safe until everyone is safe” [20–23]. WHO vaccine equity sensitization suggested that equitable roll-out of COVID-19 vaccine will not only speed up the end of the pandemic but will considerably increase population immunity worldwide, relieve burdening of health systems, and facilitate rapid recovery of economies and reduction of the likelihood of new variants emerging [24]. By September 2021, COVID-19 vaccines were widely accessible by both HICs and LMICs. However, even with vaccines being available in both HICs and LMICs, other pertinent issues adversely impacted vaccine roll-out within the countries populace. Some of the fundamental issues being inequality in vaccine distribution, vaccine prioritization, and skepticism among different groups. Limited production of COVID-19 vaccines led to many countries providing vaccines to the most vulnerable cohorts. For instance, WHO, United Kingdom, and the United States prioritized offering vaccines to healthcare personnel, individuals with high risk of developing severe symptoms that can lead to death, the elderly, and individuals with co-morbidities [13, 25, 26]. Vaccine distribution inequality was also observed even within a particular cohort (e.g., among health care workers). Studies by Chen and Krieger [27] Ogedegbe et al. [28] and Vahidy et al. [29] signaled that there is considerable inequalities in COVID-19 burden in terms of race, ethnicity, and socio-economic status which is impacted by distribution and prioritization of vaccinations. They suggested that barriers that hinder vaccine distribution should be identified and removed [30]. Although there is gradual reduction of vaccine distribution inequality, the number of individuals vaccinated in developed countries is 69 times higher than in developing countries [22].

In addition to COVID-19 vaccine being attributed to COVID-19 reduction, individual willingness to receive COVID-19 jab significantly influenced vaccine distribution [15]. Although the desire to receive COVID-19 vaccine remains relatively high globally, obstacles such as vaccine skepticism and hesitancy hamper efforts to curb COVID-19 pandemic [31]. Vaccine skepticism is based on several factors. Some individuals refuse to be vaccinated as a result of vaccine safety concerns, especially the conviction that the vaccine is dangerous due to the speed at which it was produced [15]. Other individuals argued that the COVID-19 was harmless, and therefore, they did not see the need of receiving COVID-19 vaccine. Factors such as doubt of COVID-19 vaccine effectiveness, belief in pre-existing immunity, and the origin of the vaccine also contributed in vaccine skepticism [15, 31]. Some individuals reasoned that the long-term impact of COVID-19 vaccines was unknown and therefore refused to be vaccinated. Moreover, political reasons and conspiracy theories found in many social media fuelled vaccine hesitancy [32]. Factors such as detrimental impact of COVID-19 on the lives of milliards of people around the globe and devastation of the world economy were thought to be determinants that will fuel ubiquitous acceptance of COVID-19 vaccine [33]. However, this was not the case as suggested by many findings. The study by Neumann-Bohme and co-workers [34] assessed the level of anti-COVID-19 vaccination among citizens on seven European countries and found that hesitancy and resistance against the vaccines were substantial across all cohorts, age groups, and sexes. For instance, 38% of those interviewed in France were hesitant toward COVID-19 vaccine where 10% were opposed to COVID-19 vaccination while the remaining 28% were undecided about getting COVID-19 vaccination [33]. In the United States (US), a report by the American Medical Association Journal hinted that skepticism toward vaccines that were being administered was on the rise [35]. This was further corroborated by the Kaiser Family Foundation COVID-19 vaccine monitoring which indicated that approximately 29% of health practitioners were reluctant to accept the vaccine [36].

During the early onset of COVID-19 (i.e., before accumulation of enough data to predict and project on COVID-19 prevalence), COVID-19 mathematical models proved to be an invaluable tool (and they still are) on providing insights on the appropriate mitigation strategy to adopt to prevent COVID-19 spiraling. Many models that were formulated focused on elucidating how NPIs will impact COVID-19 transmission dynamics. For instance, COVID-19 models by Ngonghala et al. [37], Srivastav et al. [88], and Wangari et al. [38] showed that correct use of face masks by public is always beneficial. However, the studies emphasized that the beneficial impact of face masks coverage at population-level is linearly dependent on face mask efficacy. Furthermore, their results showed that if the public combined both face masks and physical/social distancing, then a significant decline in COVID-19 burden occurred. The article by Tiwari et al. [39] analyzed whether disseminating COVID-19 awareness at individual, community, and global level can influence general populace viewpoint and behavior toward COVID-19. Their findings showed that higher dissemination of COVID-19 awareness among susceptible population (in particular mouth-to mouth awareness dissemination) played a key role in preventing COVID-19 spread. Moreover, an increase in awareness rate led to a significant decline of asymptomatic individuals as well as number of symptomatic individuals requiring hospitalization. A closely related study by Rai et al. [40] assessed the effectiveness of social media advertisements in suppressing rates of interaction between susceptible individuals and COVID-19 infectious individuals in India. Their results showed that airing COVID-19 information through social media platform could play a vital role in suppressing COVID-19 transmission and therefore advocated such approach to be included as a mitigation strategy. Given the study conducted by Tiwari et al. [39] was done during early development of COVID-19 when there was insufficient information on COVID-19 epidemiological dynamics, their study did not capture the role pre-symptomatic infectious individuals could have played as far as COVID-19 prevalence was concerned. Moreover, the study by Tiwari et al. [39] did not include vaccination and vaccine hesitancy, perhaps due to the reason that vaccines were in their early clinical trial stage. Similarly, the article Rai et al. [40] did not consider a scenario where pre-symptomatic infectious individuals could contribute in COVID-19 transmission and also neglected vaccination as an intervention strategy, and therefore, they could not analyze impact of vaccine hesitancy in the general populace.

Shortly after WHO and Center for Disease Control (CDC) declared availability of COVID-19 vaccines, several epidemiological mathematical frameworks began to emerge with the aim of investigating whether vaccines will slow or end the COVID-19 pandemic as well as determining strategies that ensured optimal distribution of COVID-19 vaccines [41–44]. In the sequel, few epidemiological mathematical models have attempted to elucidate on the impact vaccine hesitancy played in the early roll-out of COVID-19 vaccines [33, 35]. Some pertinent questions that relate to vaccine hesitancy have been investigated as far as COVID-19 vaccine roll-out is concerned. For example, authors in Oduro et al. [35] proposed a model that stratified the general populace depending on their willingness and unwillingness to receive COVID-19 jab. Their study exclusively investigated how educating the unwilling population about the importance of COVID-19 vaccine will impact COVID-19 transmission dynamics. Their findings suggested that education of the unwilling individuals such that they accept to get COVID-19 vaccines led to a significant reduction of control reproduction number. Although their study mentioned about vaccine hesitancy as the motivation behind incorporating education parameter among unwilling cohorts, their study did not incorporate parameters that explicitly modeled vaccine hesitancy. Moreover, vaccine hesitancy may occur in several stages especially in a scenario where COVID-19 vaccine are offered at certain intervals. For example, some vaccines such as AstraZeneca required individuals to get two doses of vaccines before one is considered to be fully vaccinated. In some cases, some individuals refused to get second dose of vaccine with the justification that they experienced severe side effects. The study by Oduro et al. [35] never captured vaccine hesitancy among partially vaccinated individuals. Buonomo et al. [33] considered a behavioral epidemiology approach to model how vaccination campaign impacted COVID-19 dynamics. The study incorporated vaccination rate as a function of the present and past information that the general populace have on the spread of the COVID-19 pandemic. Their approach was similar to the one used by dOnofrio and Manfredi [45] where they showed that incorporation of information-related changes in contact patterns induces oscillations in the endemic prevalence of infectious diseases. Similarly, the study by Buonomo et al. [33] did not consider a possibility where vaccine hesitancy could occur among different cohorts, given vaccines are also administered in phases (first dose, second dose and booster shots).

The recent study by Peter et al. [46] proposed a mathematical model that incorporated a double-dose vaccination strategy. In their study Peter et al. [46], they assumed that vaccine is administered in two phases (first dose followed by the second dose). One major assumption in their study was that fully vaccinated individuals acquired a complete vaccine-induced immunity, and therefore, these individuals could neither experience breakthrough infections nor waning of the vaccine-induced immunity. Moreover, vaccine hesitancy was not accounted for in their model. Furthermore, despite mentioning that pre-symptomatic individuals greatly contribute in the spread of COVID-19, their force of infection did not explicitly capture transmission by pre-symptomatic individuals. A single-dose vaccination strategy was not considered in Peter et al. [46], despite existence of a concrete evidence that some vaccines require individuals to receive a single dose to be considered fully vaccinated against COVID-19 [47]. For example, Johnson & Johnson COVID-19 vaccine required individuals to receive a single dose. Scientific evidence demonstrated that Johnson & Johnson vaccine rendered 85% effectiveness in preventing individuals from developing severe COVID-19 symptoms that led to COVID-19-related hospitalizations and mortality [47]. Kumar et al. [48] have recently investigated how contamination of environment and surfaces with COVID-19 could influence COVID-19 transmission dynamics in the presence of vaccination. Although the study [48] incorporated vaccination (as a mitigation strategy) without distinguishing between a single-dose or double-dose vaccine, they did not consider the epidemiological implication of vaccine hesitancy as far as COVID-19 pandemic is concerned.

The novelty of this study stems from first recognizing that pre-symptomatic infectious individuals contribute significantly in the spread of COVID-19 [8], and therefore, the impact on how pre-symptomatic individuals affect COVID-19 transmission dynamics needs to be closely examined. Second, our study is the first to consider a scenario where vaccine hesitancy is assumed to occur in two phases depending on individual's choice of vaccine type (single dose or double dose), that is, first phase of vaccine hesitancy is assumed to occur among all susceptible individuals regardless of the choice of vaccine type. The second phase of vaccine hesitancy is only assumed to occur among individuals who choose double-dose vaccine. We will, therefore, investigate epidemiological implication likely to be witnessed if individuals who choose double-dose vaccine refuse second dose and therefore remain partially vaccinated. Noteworthy partially vaccinated individuals potentially pose a greater risk to COVID-19 mitigation measures than fully vaccinated individuals [49, 50]. Moreover, this study comprehensively compares how vaccine hesitancy against either single dose or double dose influences COVID-19 transmission dynamics. To the best of our knowledge, this comparison was not done in previous studies [33, 35, 41–44, 46, 48]. Findings in this study could form the basis for policymakers and medical practitioners in deciding the type of vaccine (single dose or double dose) to administer to a given community if the levels of vaccine hesitancy are known to be either relatively low/high.

Thus, our present study recognizes the gaps and relax the assumptions made in the studies [33, 35, 41–44, 46]. Consequently, the study considers a scenario where both single-dose vaccine and a double-dose vaccination strategies are concurrently adopted. Both single- and double-dose vaccination strategies are assumed to be hindered by vaccine hesitancy in the general populace. In this study we aim to:

1. Investigate whether pre-symptomatic individuals played a significant role as far as COVID-19 transmission dynamics is concerned?

2. Determine which type of vaccine (single dose or double dose) is most suitable in a setting where vaccine hesitancy is relatively low/high?

3. Determine the most cost-effective COVID-19 mitigation strategy, especially in a scenario where both NPIs and pharmaceutical control measures are concurrently being advocated by policymakers and medical practitioners (this will be achieved by extending the proposed model using optimal control theory)?

2 Model construction

In this section, a model framework that describes COVID-19 transmission dynamics in a population where vaccine is available is presented. We assume that the government has not implemented any policy that requires all individuals to be vaccinated against COVID-19, that is, vaccination is non-mandatory. Therefore, we will consider a scenario where some individuals who have received the first COVID-19 vaccine have a choice to refuse the second dose (remain partially vaccinated). Moreover, we have considered a possibility where some individuals receive a single-dose vaccine rather than a two phase vaccine (double dose vaccine). Hence, we shall assume that COVID-19 vaccines available in the population are administered at different phases. First, we introduce notations that are necessary for the classification of individuals in the population depending on their health status. At any time t, we let N denotes the total population. N is divided into eight mutually exclusive compartments:

1. Susceptible, S: individuals who are considered to be healthy and are liable to contract COVID-19,

2. Vaccinated, V_p: individuals who have received the first dose of COVID-19 vaccine and are assumed to be partially protected. Hence, they can be infected with COVID-19 although at a reduced rate in comparison with susceptible individuals,

3. Vaccinated, V: individuals who have received either double-dose COVID-19 vaccine or the single-dose vaccine and are therefore considered to be fully vaccinated against COVID-19, but they can experience breakthrough infections at a relatively lower rate compared to both partially vaccinated and susceptible individuals,

4. Exposed E: consists of individuals who have been exposed to SARS-CoV-2, but are not capable of passing the virus to other healthy individuals, reason being that the virus incubation period has not yet elapsed,

5. Pre-symptomatic, I_p: This cohort includes individuals who are at post-latency stage (i.e., individuals who are in the incubation phase following latency). They are asymptomatic and able to transmit SARS-CoV-2,

6. Asymptomatic infectious I_a: This cohort includes individuals who remain asymptomatic throughout the course of the infection, but are able to pass the SARS-CoV-2 to others. Often regarded as truly asymptomatic individuals,

7. Symptomatic infectious, I_s: Include individuals who manifest either mild or severe COVID-19 symptoms and are able to spread SARS-CoV-2,

8. Recovered, R: This cohort includes those who have recovered after infectious period.

The size of each class at any given time t describes the state variable of the model framework and N = S + V_p + V + E + I_p + I_a + I_s + R. For the proposed model framework, we shall assume both a double-dose vaccine and a single-dose vaccines are concurrently being administered, contrary to study [33] that considered a single-dose vaccine which rendered a lifelong protection against COVID-19. It is assumed that the vaccine-induced immunity is not perfect and can wane over time. Therefore, vaccine-induced immunity among individuals who are fully vaccinated against COVID-19 and individuals who have received one dose of COVID-19 vaccine can wane over time such that it is of the same level as susceptible individuals; hence, they are assumed to join S compartment. Regarding transmission of COVID-19 among susceptible, partially vaccinated and fully vaccinated individuals, we assume that infection occurs as a result of interaction with three infectious cohorts, namely, pre-symptomatic infectious individuals, truly asymptomatic infectious individuals, and symptomatic infectious individuals. Transmission of COVID-19 by pre-symptomatic individuals is based on the findings that individuals exposed to COVID-19 are most infectious during the post-latency stage when they are not manifesting COVID-19 symptoms (usually 1–2 days before symptoms set in) [2, 8]. Research on this finding is still ongoing. On the one hand, truly asymptomatic individuals transmit COVID-19 at a relatively low rate in comparison with pre-symptomatic individuals. This assumption corroborates well with the finding that truly asymptomatic patients cause little or no transmission of COVID-19 [51].

In the proposed COVID-19 model, the force of infection which represents the per capita rate at which susceptible, partially vaccinated individuals and fully vaccinated individuals contract COVID-19 is assumed to be a mass action incidence defined as λ = β(I_p + η₁I_a + η₂I_s). The justification behind choosing a mass action incidence is a 2-fold. One being that during the COVID-19 pandemic, the total population remained relatively constant [33]. Second, mass action incidence is chosen for mathematical simplification and tractability. Note that in the force of infection, parameters η₁ and η₂ are modification coefficients that account for reduced infectiousness of truly asymptomatic and symptomatic individuals when compared with infectious post-latency individuals (pre-symptomatic individuals). Parameter β represents the effective contact rate. All subgroups decrease due to natural mortality at per capita rate μ per day, while infectious symptomatic individuals experience an extra mortality due to COVID-19-induced death at per capita rate d per day.

Recruitment into susceptible subgroup S occurs at a rate ∧ humans per day. Since first dose vaccine-induced immunity is not lifelong (imperfect to protect against COVID-19), partially vaccinated individual immunity can wane at a rate ω₁ per day and these individuals revert to susceptible class. Similarly, fully vaccinated individual immunity can wane over time, and these individuals return to susceptible subgroup at a rate ω₂ per day. Due to myriad challenges that faced COVID-19 vaccine distribution during the early days, among them being vaccine hesitancy (just to mention a few) the rates at which susceptible individuals in the class S, partially vaccinated individuals in the subgroup V_p, and fully vaccinated individuals V are vaccinated are assumed to occur at rates (1 − h_i)θ_i, where θ_i,   i ∈ (1, 2, 3) represent rates at which vaccine is administered (respectively, among subgroups S, V_p, and V) and 0 ≤ h_i < 1,   i ∈ (1, 2, 3) represent fractions of susceptible and partially vaccinated individuals who are unwilling to get vaccinated (vaccine hesitancy) by the double dose vaccine (which includes first dose and second dose) or single dose vaccine. Thus, susceptible sub-population decreases as a result of vaccination at rates (1–h₁)θ₁ and (1–h₃)θ₃ per day and also due to infection at a rate λ per day. The subgroup V_p increases due to vaccination of susceptible individuals at a rate (1 − h₁)θ₁ per day. The subgroup V_p decreases as a result of receiving second dose of COVID-19 vaccine at a rate (1 − h₂)θ₂ per day, waning of vaccine-induced immunity at a rate ω₁ per day and also due to breakthrough infections at a rate (1 − σ)λ per day. The parameter σ, 0 ≤ σ < 1 represents vaccine efficacy among subgroup V_p who have received the first dose. The sub-population V is generated as a result of vaccination of both partially vaccinated individuals (double dose) and susceptible individuals (single dose) at rates (1 − h₂)θ₂ and (1 − h₃)θ₃ per day, respectively. They decrease due to vaccine waning at a rate ω₂ per day and breakthrough infections at a rate (1 − ψ)λ per day. The parameter ψ, 0 ≤ ψ < 1 represents vaccine efficacy among subgroup V (who have either received a complete double-dose vaccine or single-dose vaccine). The exposed sub-population arise due to infection of susceptible population at a rate λ and also as a result of breakthrough infections among partially vaccinated individuals at a rate (1 − σ)λ and fully vaccinated individuals at a rate (1 − ψ)λ per day. The exposed individuals progress to the infectious post-latency compartment I_p at per capita rate ϕ per day (note, 1/ϕ is the average incubation period/latent period). Pre-symptomatic infectious individuals progress to either infectious truly asymptomatic individuals at per capita rate (1 − v)α per day or to the infectious symptomatic individuals at per capita rate vα per day. Note that 0 < v < 1(0 < 1 − v < 1) represents the fraction of post-latency individuals (infectious pre-symptomatic persons) who either manifest or do not manifest COVID-19 symptoms at the end of the intrinsic incubation period and 1/α represents the average pre-symptomatic infectious period [52]. Truly infectious asymptomatic and infectious symptomatic individuals recover from COVID-19 and progress to subgroup R at per capita rates γ_a and γ_s, respectively (where 1/γ_a and 1/γ_s are the average durations of infectiousness for I_a and I_s subgroups, respectively).

For biological plausibility of our model, we will assume that vaccine-induced immunity strength among fully vaccinated individuals in the subgroup V is stronger in comparison with partially vaccinated individuals in the class V_p. Consequently, vaccine-induced immunity waning rate among fully vaccinated individuals will be considered to be slightly lower/equal to vaccine waning rate among partially vaccinated individuals (i.e., ω₁ ≥ ω₂). This assumption is made based on recent evidence regarding breakthrough infections among fully vaccinated individuals [53]. A schematic representation of the model is shown in Figure 1.

FIGURE 1

Figure 1. Schematic representation of COVID-19 model that incorporates a single-dose and a double-dose vaccination strategies that are coupled with vaccine hesitancy. The population N(t) is stratified into eight mutually exclusive classes: susceptible S(t), partially vaccinated with one dose V_p, fully vaccinated with either double-dose or single-dose vaccine V(t), exposed E(t), pre-symptomatic I_p(t), truly asymptomatic I_a(t), symptomatic I_s(t), and recovered R(t). Blue dashed arrows indicate breakthrough infections.

2.1 The model equations

Considering the aforementioned model description and assumptions, the system of non-linear differential equations (1) governs the transmission dynamics of COVID-19 in the presence of a single-dose and two-dose vaccines that are non-mandatory and takes into account vaccine hesitancy. A detailed description of the model parameters is captured in Table 1.

$\begin{array}{l}\frac{dS}{dt}=\wedge +{\omega }_1{V}_p+{\omega }_{2}V-(\mu +\lambda +(1-h_1){\theta }_1+(1-h_3){\theta }_3)S,\\ \frac{d{V}_p}{dt}=(1-h_1){\theta }_{1}S-(\mu +{\omega }_1+(1-h_2){\theta }_2+(1-\sigma )\lambda )V_p,\\ \frac{dV}{dt}=(1-h_2){\theta }_2{V}_p+(1-h_3){\theta }_{3}S-(\mu +{\omega }_2+(1-\psi )\lambda )V,\end{array}$

$\begin{array}{ll}\frac{dE}{dt}=\lambda S+(1-\sigma )\lambda V_p+(1-\psi )\lambda V-(\mu +\varphi )E,& (1)\end{array}$

$\begin{array}{l}\frac{d{I}_p}{dt}=\varphi E-(\mu +\alpha )I_p,\\ \frac{d{I}_a}{dt}=(1-v)\alpha I_p-(\mu +{\gamma }_a)I_a,\\ \frac{d{I}_s}{dt}=v\alpha I_p-(\mu +d+{\gamma }_s)I_s,\\ \frac{dR}{dt}={\gamma }_a{I}_a+{\gamma }_s{I}_s-\mu R,\end{array}$

where λ remains as previously defined. The model equation (1) is subject to the initial conditions

$\begin{array}{l}{S}_0=S(0)\ge 0,V_{p0}=V_p(0)\ge 0,V_0=V(0)\ge 0,E_0=E(0)\ge 0,\end{array}$

$\begin{array}{ll}{I}_{p0}=I_p(0)\ge 0,I_{a0}=I_a(0)\ge 0,I_{s0}=I_s(0)\ge 0,R_0=R(0)\ge 0.& (2)\end{array}$

Note that the dynamics of model equation (1) is not affected by the last equation that models the time evolution of the recovered class. Hence, the qualitative dynamics of the subgroup R will not be given much attention in the theoretical analysis.

TABLE 1

Table 1. Description of model baseline parameter values.

3 Analytical results of the model

In this section, the model is qualitatively analyzed. First, it is shown that the model state variables remain non-negative for all time t > 0 whenever non-negative initial conditions are supplied into the system (1). Boundedness of the model trajectories will be established. The intrinsic COVID-free equilibrium and the fundamental threshold that determine the severity of COVID-19 (control reproduction number) will be determined. It will be shown that the disease-free equilibrium of the model described in section 2 is locally (globally) asymptotically stable whenever the control reproduction number is less than unity. Finally, the endemic equilibrium point together with the relevant properties (local stability) will be established.

3.1 Non-negativity of the model solutions

Given the proposed model considers a human population, it is important to show the well possessedness of the model. The following theorem establishes that the model solutions remain non-negative over all the time.

Theorem 1. If the model equations (1) is supplied with positive initial conditions; S(0), V_p(0), V(0), E(0), I_p(0), I_a(0), I_s(0), and R(0), then the model state variables; S(t), V_p(t), V(t), E(t), I_p(t), I_a(t), I_s(t), and R(t) will always remain positive for all time, t > 0.

Proof: Suppose

$\begin{array}{l}T=\sup \{\tau \ge 0,\mid \forall t,0\le t\le \tau \text{ such that }S(t)\ge 0,V_p(t)\ge 0,\\ V(t)\ge 0,E(t)\ge 0,I_p(t)\ge 0,I_a(t)\ge 0,I_s(t)\ge 0,R(t)\ge 0\}.\end{array}$

First let 0 < T < +∞. Now considering continuity of solutions, we can have S(T) = 0 or V_p(T) = 0 or V(T) = 0 or E(T) = 0 or I_p(T) = 0 or I_a(T) = 0 or I_s(T) = 0 or R(T) = 0. Supposing S(T) = 0 before all other state variables V_p, V, E, I_p, I_a, I_s, R become zero. Then from the first equation of model system (1) we have

$\begin{array}{l}\frac{dS(T)}{dt}{|}_{\left\{S(T)=0,V_p>0,V>0,E>0,I_p>0,I_a>0,I_s>0,R>0\right\}}\\ =\wedge +{\omega }_1{V}_p+{\omega }_{2}V>0\text{ since }\wedge ,{\omega }_1,{\omega }_2>0.& (3)\end{array}$

If V_p(T) becomes zero before all other state variables S, V, E, I_p, I_a, I_s, R, then from the second equation of model system (1) we have

$\begin{array}{l}\frac{d{V}_p(T)}{dt}{|}_{\left\{S>0,V_p(T)=0,V>0,E>0,I_p>0,I_a>0,I_s>0,R>0\right\}}\\ =(1-h_1){\theta }_{1}S>0\text{since }{h}_1<1,{\theta }_1>0.& (4)\end{array}$

Following a similar procedure as described in equations (3) and (4), it is not difficult to deduce from third, fourth, fifth, sixth, seventh, and eighth equations of model system (1) the following results;

$\begin{array}{r}\frac{dV(T)}{dt}{|}_{\left\{S>0,V_p>0,V(T)=0,E>0,I_p>0,I_a>0,I_s>0,R>0\right\}}\\ =(1-h_2){\theta }_2{V}_p+(1-h_3){\theta }_{3}S>0\text{since }\\ h_2,h_3<1,{\theta }_2,{\theta }_3>0,\\ \frac{dE(T)}{dt}{|}_{\left\{S>0,V_p>0,V>0,E(T)=0,I_p>0,I_a>0,I_s>0,R>0\right\}}\\ =\lambda S+(1-\sigma )\lambda V_p+(1-\psi )\lambda V>0\\ \text{since }0<\sigma ,\psi <1,\\ \frac{d{I}_p(T)}{dt}{|}_{\left\{S>0,V_p>0,V>0,E>0,I_p(T)=0,I_a>0,I_s>0,R>0\right\}}\\ =\varphi E>0,\\ \frac{d{I}_a(T)}{dt}{|}_{\left\{S>0,V_p>0,V>0,E>0,I_p>0,I_a(T)=0,I_s>0,R>0\right\}}\\ =(1-v)\alpha I_p>0,\text{since }0<v<1,\alpha >0,\\ \frac{d{I}_s(T)}{dt}{|}_{\left\{S>0,V_p>0,V>0,E>0,I_p>0,I_a>0,I_s(T)=0,R>0\right\}}\\ =v\alpha I_p>0,\\ \frac{dR(T)}{dt}{|}_{\left\{S>0,V_p>0,V>0,E>0,I_p>0,I_a>0,I_s>0,R(T)=0\right\}}\\ ={\gamma }_a{I}_a+{\gamma }_s{I}_s>0,\end{array}$

which completes the proof by showing that all the state variables of the model system (1) remain positive for all time.

3.2 Boundedness of model solutions

Theorem 2. Define the following biologically feasible region:

$\begin{array}{r}\Delta =\{(S,V_p,V,E,I_p,I_a,I_s,R)\in {ℝ}_{+}^8\mid S+V_p+V+E+I_p\\ +I_a+I_s+R\le \frac{\wedge }{\mu }\}.\end{array}$

Then, the region Δ is positively invariant with respect to the model system (1).

Proof: Noting that the total population is given by the sum of the state variables of the proposed model, we have

$\begin{array}{ll}N(t)=S(t)+V_p(t)+V(t)+E(t)+I_p(t)+I_a(t)+I_s(t)+R(t).& (5)\end{array}$

The derivative of equation (5) is given by

$\begin{array}{l}\frac{dN(t)}{dt}=\frac{dS(t)}{dt}+\frac{d{V}_p(t)}{dt}+\frac{dV(t)}{dt}+\frac{dE(t)}{dt}+\frac{d{I}_p(t)}{dt}+\frac{d{I}_a(t)}{dt}\\ +\frac{d{I}_s(t)}{dt}+\frac{dR(t)}{dt}\\ =\wedge -\mu N(t)-d_s{I}_s.& (6)\end{array}$

From equation (6), we have $\frac{dN(t)}{dt}\le \wedge -\mu N(t)$ which can be written as $\frac{dN(t)}{dt}+\mu N(t)\le \wedge .$ Applying integration factor method and after algebraic manipulation leads to

$\begin{array}{ll}N(t)\le \frac{\wedge }{\mu }+\left(N(0)-\frac{\wedge }{\mu }\right)e^{-\mu t}.& (7)\end{array}$

From (7), we can deduce that $\underset{t\to \infty }{lim}supN(t)\le \frac{\wedge }{\mu }.$ Hence, it follows that $N(t)\le \frac{\wedge }{\mu }$ as t → ∞. Thus, all the solutions of model equation (1) remain in Δ for all time which implies that the model system (1) is positively invariant. Hence, all the analysis conducted on the model system (1) are limited in the region Δ, where the model equation (1) is considered as being mathematically and epidemiologically meaningful.

3.3 COVID-free equilibrium and control reproduction number

In the absence of any COVID-19 infection in the community (i.e., E= I_p = I_a = I_s = 0), the model system (1) has an intrinsic COVID-free equilibrium $\mathcal{CFE}$ which is obtained by setting the right-hand side of equation (1) to zero. This results to

$\begin{array}{ll}\mathcal{CFE}=\left(S^0,V_p^0,V^0,0,0,0,0,0\right),& (8)\end{array}$

where

$\begin{array}{l}{S}^0=\frac{D_1{D}_2\wedge }{K_1+K_2{D}_2},V_p^0=\frac{(1-h_1){\theta }_1{D}_1\wedge }{K_1+K_2{D}_2},\\ V^0=\frac{(1-h_1){\theta }_1(1-h_2){\theta }_2\wedge +(1-h_3){\theta }_3{D}_2\wedge }{K_1+K_2{D}_2},\text{where}\end{array}$

$\begin{array}{ll}{D}_1=\mu +{\omega }_2,D_2=(\mu +{\omega }_1+(1-h_2){\theta }_2),& (9)\end{array}$

$\begin{array}{l}{K}_1=\mu (1-h_1){\theta }_1(\mu +{\omega }_2+(1-h_2){\theta }_2),\\ K_2=\mu (\mu +{\omega }_2+(1-h_3){\theta }_3).\end{array}$

The basic reproduction number often denoted by ${\mathcal{R}}_0$ is one of the important thresholds in mathematical epidemiology due to its ability in predicting the transmission dynamics of an infectious disease. It is defined as the number of secondary infectives triggered by a single primary infection during infectious period when introduced into an entirely susceptible population. For the proposed COVID-19 model equation (1), instead of computing reproduction number ${\mathcal{R}}_0$ we compute the control reproduction number ${\mathcal{R}}_c$ due to the reason that intervention measures for mitigating the spread of the pandemic are incorporated (vaccination strategies: single-dose and double-dose vaccines). Using the next-generation approach originally described in Van den Driessche and Watmough [65], we compute the control reproduction number as the dominant eigenvalue (spectral radius) of the next-generation matrix ${\mathcal{FV}}^{-1}.$ The matrices $\mathcal{F}$ and $\mathcal{V},$ respectively, represent Jacobian matrices of the terms that involve new infections and the transition rates (from one healthy status to another) computed at the COVID-free equilibrium $\mathcal{CFE}$ [see Van den Driessche and Watmough [65]]. Thus, we have

$\mathcal{F}=\left[\begin{array}{cccc}0& \beta S^0+(1-\sigma )\beta V_p^0+(1-\psi )\beta V^0& \beta \eta S^0+(1-\sigma )\beta {\eta }_1{V}_p^0+(1-\psi )\beta {\eta }_1{V}^0& \beta {\eta }_2{S}^0+(1-\sigma )\beta {\eta }_2{V}_p^0+(1-\psi )\beta {\eta }_2{V}^0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right]$

and

$\mathcal{V}=\left[\begin{array}{cccc}(\mu +\varphi )& 0& 0& 0\\ -\varphi & (\mu +\alpha )& 0& 0\\ 0& -(1-v)\alpha & (\mu +{\gamma }_a)& 0\\ 0& -v\alpha & 0& (\mu +d+{\gamma }_s)\end{array}\right].$

The control reproduction number of the model system (1) can now be determined from $\varrho ({\mathcal{FV}}^{-1})$, where ϱ is the spectral radius or the dominant eigenvalue.

Let

$\begin{array}{r}\varpi =D_2(D_1+(1-\psi )(1-h_3){\theta }_3)+(1-h_1){\theta }_1((1-\sigma )D_1\\ +(1-\psi )(1-h_2){\theta }_2),\text{ then}\end{array}$

$\begin{array}{ll}{\mathcal{R}}_c={\mathcal{R}}_p+{\mathcal{R}}_a+{\mathcal{R}}_s,& (10)\end{array}$

where

$\begin{array}{l}{\mathcal{R}}_p=\frac{\varpi \beta \wedge \varphi }{(\mu +\varphi )(\mu +\alpha )(K_1+K_2{D}_2)},\\ {\mathcal{R}}_a=\frac{\varpi \beta \wedge \varphi (1-v)\alpha {\eta }_1}{(\mu +\varphi )(\mu +\alpha )(\mu +{\gamma }_a)(K_1+K_2{D}_2)},\\ {\mathcal{R}}_s=\frac{\varpi \beta \wedge \varphi v\alpha {\eta }_2}{(\mu +\varphi )(\mu +\alpha )(\mu +d+{\gamma }_s)(K_1+K_2{D}_2)}.\end{array}$

Following the theorem proved in Van den Driessche and Watmough [65], we state the following lemma.

Lemma 1. The COVID-free equilibrium $\mathcal{CFE}$ of the model system (1) is LAS (locally asymptotically stable) when ${\mathcal{R}}_c=\varrho ({\mathcal{FV}}^{-1})<1$ and unstable if ${\mathcal{R}}_c=\varrho ({\mathcal{FV}}^{-1})>1.$

It is easy to note that the control reproduction number (10) is a sum of three terms, namely, ${\mathcal{R}}_p,{\mathcal{R}}_a$, and ${\mathcal{R}}_s.$ Epidemiologically, ${\mathcal{R}}_p$ represents the reproduction number of infectious individuals generated due to infection with pre-symptomatic infectious individuals in the subgroup I_p, ${\mathcal{R}}_a$, accounts for the reproduction number of infectious cases produced by infection with truly asymptomatic infectious individuals in the compartment I_a, and ${\mathcal{R}}_s$ is the reproduction number associated with the number of infection cases generated as a result of infection with symptomatic infectious individuals in the subgroup I_s. Mathematically, ${\mathcal{R}}_a$ and ${\mathcal{R}}_s$ can be expressed in terms of ${\mathcal{R}}_p$, that is, ${\mathcal{R}}_a=\frac{(1-v)\alpha {\eta }_1}{(\mu +{\gamma }_a)}{\mathcal{R}}_p$ and ${\mathcal{R}}_s=\frac{v\alpha {\eta }_2}{(\mu +d+{\gamma }_s)}{\mathcal{R}}_p.$ Consequently,

$\begin{array}{l}{\mathcal{R}}_c={\mathcal{R}}_p+\frac{(1-v)\alpha {\eta }_1}{(\mu +{\gamma }_a)}{\mathcal{R}}_p+\frac{v\alpha {\eta }_2}{(\mu +d+{\gamma }_s)}{\mathcal{R}}_p,\\ ={\mathcal{R}}_p\left(1+\frac{(1-v)\alpha {\eta }_1}{(\mu +{\gamma }_a)}+\frac{v\alpha {\eta }_2}{(\mu +d+{\gamma }_s)}\right),& (11)\end{array}$

$\begin{array}{l}=\frac{Q_0}{(\mu +{\gamma }_a)(\mu +d+{\gamma }_s)}{\mathcal{R}}_p,\text{where}\end{array}$

$\begin{array}{l}{Q}_0=(\mu +{\gamma }_a)(\mu +d+{\gamma }_s)+(1-v)\alpha {\eta }_1(\mu +d+{\gamma }_s)\\ +\mu \alpha {\eta }_2(\mu +{\gamma }_a)>0.\end{array}$

For biologically plausible non-negative parameter values (see Table 1), both $\frac{(1-v)\alpha {\eta }_1}{(\mu +{\gamma }_a)}$ and $\frac{v\alpha {\eta }_2}{(\mu +d+{\gamma }_s)}$ are less than one. If both $\frac{(1-v)\alpha {\eta }_1}{(\mu +{\gamma }_a)}⋘1$ and $\frac{v\alpha {\eta }_2}{(\mu +d+{\gamma }_s)}⋘1,$ then ${\mathcal{R}}_c\approxeq {\mathcal{R}}_p.$ This indicates that the reproduction number ${\mathcal{R}}_p$ associated with pre-symptomatic individuals is larger than both ${\mathcal{R}}_a$ and ${\mathcal{R}}_s.$ Epidemiologically, this implies that pre-symptomatic individuals pose a greater risk as far as the spread of COVID-19 is concerned. This was not captured in several mathematical models due to the fact that it was not clear whether pre-symptomatic individuals played any significant role in the spread of COVID-19, and therefore, transmission of COVID-19 by pre-symptomatic individuals was overlooked. However, now there is a global consensus (and research still being conducted) that pre-symptomatic individuals (individuals at post-latency stage) are more infectious before symptoms start to manifest [see WHO [2] and Alleman et al. [8] and also the recent research that incorporated pre-symptomatic transmission in the force of infection [33]].

Lemma 1 implies that there is a possibility of COVID-19 pandemic being eradicated in case the initial sizes of the state variables of model equation (1) are in the basin of attraction of the COVID-free equlibrium $\mathcal{CFE}.$ However, there is a caveat regarding local stability of $\mathcal{CFE}$, given the proposed model system (1) assumes the COVID-19 vaccines (both single dose and double dose) being administered are not perfect. To be certain COVID-19 pandemic eradication is not dependent on the initial sizes of the state variables, we investigate whether the $\mathcal{CFE}$ is globally asymptotically stable (if not, then we cannot rule out occurrence of other complex bifurcation structures, particularly backward bifurcation which occurs when ${\mathcal{R}}_c<1$).

3.4 Global stability of $\mathcal{CFE}$

To analyze the global stability of model equation (1), we apply the procedure described by Castillo-Chavez et al. [66], where the two conditions that need to be fulfilled for the global asymptotic stability of the disease-free equilibrium are outlined. First, we re-write the model system (1) into the following form:

$\begin{array}{ll}{Z}^{{}^{\prime }}=G_1(Z,P),& (12)\end{array}$

$\begin{array}{ll}{P}^{{}^{\prime }}=G_2(Z,P),G_2(Z,0)=0,& (13)\end{array}$

where $Z=(S,V_p,V,R)\in {ℝ}^4,$ represents the number of uninfected state variables, and $P=(E,I_p,I_a,I_s)\in {ℝ}^4$ represents the number of infected state variables for model system (1). Let $W_0=(Z^{*},0)$ define the COVID-free equilibrium of the proposed model (1). Then, the conditions outlined below must be fulfilled for the model system (1) to be globally asymptotically stable.

1. For $Z^{{}^{\prime }}=G_1(Z,0),$     Z^* is g.a.s (globally asymptotically stable).

2. G₂(Z, P) = AP − Ĝ₂(Z, P),     Ĝ₂(Z, P) ≥ 0, where

$A=D_P{G}_2(Z^{*},0)$ represents a matrix whose off-diagonal entries are non-negative (M-matrix). Supposing conditions (i) and (ii) stated above hold in model system (1), then the following theorem follows:

Theorem 3. Assuming conditions (i) and (ii) above hold and ${\mathcal{R}}_c<1,$ then the fixed point $W_0=(Z^{*},0)$ is a g.a.s equilibrium point of model system (1).

Proof: It is straight forward to observe that

$\begin{array}{l}{Z}^{{}^{\prime }}=G_1(Z,P)\\ =\left[\begin{array}{c}\wedge +{\omega }_1{V}_p+{\omega }_{2}V-(\mu +\lambda +(1-h_1)\theta +(1-h_3){\theta }_3)S\\ (1-h_1){\theta }_{1}S-(\mu +{\omega }_1+(1-h_2){\theta }_2+(1-\sigma )\lambda )V_p\\ (1-h_2){\theta }_2{V}_p+(1-h_3){\theta }_{3}S-(\mu +{\omega }_2+(1-\psi )\lambda )V\\ {\gamma }_a{I}_a+{\gamma }_s{I}_s-\mu R\end{array}\right],\end{array}$

$G_1(Z,0)=\left[\begin{array}{c}\wedge +{\omega }_1{V}_p+{\omega }_{2}V-(\mu +(1-h_1){\theta }_1+(1-h_3){\theta }_3)S^0\\ (1-h_1){\theta }_1{S}^0-(\mu +{\omega }_1+(1-h_2){\theta }_2)V_p^0\\ (1-h_2){\theta }_2{V}_p^0+(1-h_3){\theta }_3{S}^0-(\mu +{\omega }_2)V^0\\ 0\end{array}\right],$

$P^{{}^{\prime }}=G_2(Z,P)=\left[\begin{array}{c}\begin{array}{c}\beta (I_p+{\eta }_1{I}_a+{\eta }_2{I}_s)S+(1-\sigma )\beta (I_p+{\eta }_1{I}_a+{\eta }_2{I}_s)V_p\\ +(1-\psi )\beta (I_p+{\eta }_1{I}_a+{\eta }_2{I}_s)V-(\mu +\varphi )E\end{array}\\ \varphi E-(\mu +\alpha )I_p\\ (1-v)\alpha I_p-(\mu +{\gamma }_a)I_a\\ v\alpha I_p-\mu +d+{\gamma }_s{I}_s\end{array}\right]$

$A=\left[\begin{array}{cccc}-(\mu +\varphi )& A_{12}& A_{13}& A_{14}\\ \varphi & -(\mu +\alpha )& 0& 0\\ 0& (1-v)\alpha & -(\mu +{\gamma }_a)& 0\\ 0& v\alpha & 0& -(\mu +d+{\gamma }_s)\end{array}\right].$

where $A_{12}=\beta S^0+(1-\sigma )\beta V_p^0+(1-\psi )\beta V^0,A_{13}=\beta {\eta }_1{S}^0+(1-\sigma )\beta {\eta }_1{V}_p^0+(1-\psi )\beta V^0,$ $A_{14}=\beta {\eta }_2{S}^0+(1-\sigma )\beta {\eta }_2{V}_p^0+(1-\psi )\beta {\eta }_2{V}^0.$ Now Ĝ₂(Z, P) = AP − G₂(Z, P) is given as:

${Ĝ}_2(Z,P)=\left[\begin{array}{c}\begin{array}{c}(S^0-S)\beta (I_p+{\eta }_1{I}_a+{\eta }_2{I}_s)+(V_p^0-V_p)(1-\sigma )\beta (I_p+{\eta }_1{I}_a+{\eta }_2{I}_s)\\ +(V^0-V)(1-\psi )\beta (I_p+{\eta }_1{I}_a+{\eta }_2{I}_s)\end{array}\\ 0\\ 0\\ 0\end{array}\right]$

As t → ∞, the total population N is bounded above by $\frac{\wedge }{\mu }$ which implies that the state variables S, V_p, and V satisfy the inequalities: $0<S\le N\le \frac{\wedge }{\mu },$ $0<V_p\le N\le \frac{\wedge }{\mu },$ and $0<V\le N\le \frac{\wedge }{\mu }.$ Consequently, after algebraic manipulation, we have

$\begin{array}{l}{S}^0-S\approx S^0-\frac{\wedge }{\mu }\\ =\frac{\wedge }{\mu }\left(\frac{1}{1+\frac{(1-h_1){\theta }_1(\mu +{\omega }_2+(1-h_2){\theta }_2)+(1-h_3){\theta }_3(\mu +{\omega }_1+(1-h_2){\theta }_2)}{(\mu +{\omega }_2)(\mu +{\omega }_1+(1-h_2){\theta }_2)}}-1\right),& (14)\end{array}$

$\begin{array}{l}{V}_p^0-V_p\approx V_p^0-\frac{\wedge }{\mu }\\ =\frac{\wedge }{\mu }\left(\frac{1}{1+\frac{(1-h_1){\theta }_1(1-h_2){\theta }_2+(\mu +{\omega }_2)(1-h_3){\theta }_3(\mu +{\omega }_1+(1-h_2){\theta }_2)}{(1-h_1){\theta }_1(\mu +{\omega }_2)}}-1\right),& (15)\end{array}$

$\begin{array}{l}{V}^0-V\approx V^0-\frac{\wedge }{\mu }\\ =\frac{\wedge }{\mu }\left(\frac{1}{1+\frac{(1-h_1){\theta }_1(\mu +{\omega }_2)+(\mu +{\omega }_2)(\mu +{\omega }_1+(1-h_2){\theta }_2)}{(1-h_1){\theta }_1(1-h_2){\theta }_2+(1-h_3){\theta }_3(\mu +{\omega }_1+(1-h_2){\theta }_2)}}-1\right).& (16)\end{array}$

It is clear that the first terms inside the brackets in equations (14), (15), and (16) are all less than one, which implies $S^0-S<0,V_p^0-V_p<0,$ and V⁰ − V < 0. The implication of this is that Ĝ₂(Z, P) ≱ 0. Hence, condition (ii) for global asymptotic stability is not satisfied. Thus, the COVID-free equilibrium $\mathcal{CFE}$ is not globally asymptotically stable. This signals a possibility of occurrence of bistability phenomenon where multiple equilibria coexist when the control reproduction number is less than unity. In what follows, we investigate the persistence of COVID-19 pandemic in the population through determining the endemic equilibrium point of model system (1).

3.5 COVID-19 persistence equilibria

In this subsection, we determine the endemic equilibrium points of model system (1) by considering that the COVID-19 is persistent in the population. Thus, we set the right-hand side of equation (1) to zero and evaluate for steady states $S^{*},V_p^{*},V^{*},E^{*},I_p^{*},I_a^{*},I_s^{*}$ (note the steady state R^* has been omitted due to the fact that the model dynamics of equation (1) do not depend on state variable R). For the purpose of mathematical simplification, the endemic equilibrium is expressed in terms of the force of infection, which at equilibrium is denoted by ${\lambda }^{*}=\beta (I_p^{*}+{\eta }_1{I}_a^{*}+{\eta }_2{I}_s^{*}).$ Let D₁ and D₂ remain as previously defined in equation (9), then we define

$\begin{array}{l}{D}_3=(1-h_1){\theta }_1,D_4=\mu +(1-h_2){\theta }_2,D_5={\omega }_2(1-h_1){\theta }_1,\\ D_6=(1-h_3){\theta }_3,\\ D_7=(1-h_2){\theta }_2,\\ D_8=[(D_1+(1-\psi ){\lambda }^{*})(\mu +{\lambda }^{*})\\ +D_6(\mu +(1-\psi ){\lambda }^{*})][D_2+(1-\sigma ){\lambda }^{*}]\\ +D_3[\mu +(1-\psi ){\lambda }^{*}][D_4+(1-\sigma ){\lambda }^{*}]+D_5[\mu +(1-\sigma ){\lambda }^{*}],\\ D_9=(\mu +\varphi )(\mu +{\gamma }_a)(\mu +d+{\gamma }_s)(\mu +\alpha ).\end{array}$

Now, the steady states expressed in terms of force of infection at equilibrium are given as

$\begin{array}{l}{S}^{*}=\frac{\wedge (D_2+(1-\sigma ){\lambda }^{*})(D_1+(1-\psi ){\lambda }^{*})}{D_8},\\ V_p^{*}=\frac{\wedge D_3(D_1+(1-\psi ){\lambda }^{*})}{D_8},\\ V^{*}=\frac{\wedge (D_3{D}_7+D_6(D_2+(1-\sigma ){\lambda }^{*}))}{D_8},\end{array}$

$\begin{array}{ll}{E}^{*}={\lambda }^{*}\wedge (\frac{(D_2+(1-\sigma ){\lambda }^{*})(D_1+(1-\psi ){\lambda }^{*})+(1-\sigma )D_3(D_1+(1-\psi ){\lambda }^{*})+(1-\psi )(D_3{D}_7+D_6(D_2+(1-\sigma ){\lambda }^{*}))}{D_8}),& (17)\end{array}$

$\begin{array}{l}{I}_p^{*}=\frac{\varphi }{(\mu +\alpha )(\mu +\varphi )D_8}{E}^{*},\\ I_a^{*}=\frac{(1-v)\alpha }{\mu +{\gamma }_a}{I}_p^{*},\\ I_s^{*}=\frac{v\alpha }{(\mu +d+{\gamma }_s)}{I}_p^{*}.\end{array}$

Substituting $I_p^{*},I_a^{*}$, and $I_s^{*}$ into the force of infection at equilibrium yields the following fourth degree polynomial expressed as a function of λ^*

$\begin{array}{ll}g({\lambda }^{*})={\lambda }^{*}(C_3{\lambda }^{*3}+C_2{\lambda }^{*2}+C_1{\lambda }^{*}+C_0)=0,& (18)\end{array}$

where

$\begin{array}{ll}\begin{array}{l}{C}_3=D_9(1-\sigma )(1-\psi ),\\ C_2=\mu D_9(1-\psi )(1-\sigma )+D_1{D}_9(1-\sigma )+D_2{D}_9(1-\psi )\\ +(1-\psi )(1-\sigma )D_6{D}_9+D_3{D}_9(1-\psi )(1-\sigma )\\ -(\mu +{\gamma }_a)(\mu +d+{\gamma }_s)\beta \varphi \wedge (1-\sigma )(1-\psi )\\ -(\mu +d+{\gamma }_s)\beta {\eta }_1(1-v)\alpha \varphi \wedge (1-\sigma )(1-\psi )\\ -(\mu +{\gamma }_a)\beta {\eta }_{2}v\alpha \varphi \wedge (1-\sigma )(1-\psi ),\\ C_1=\mu D_9{D}_1(1-\sigma )+\mu D_9{D}_2(1-\psi )+D_1{D}_2{D}_9\\ +\mu (1-\sigma )D_6{D}_9+D_2{D}_6{D}_9(1-\psi )+D_3{D}_9\mu (1-\sigma )\\ +D_3{D}_4{D}_9(1-\psi )+D_5{D}_9(1-\sigma )\\ -(\mu +{\gamma }_a)(\mu +d+{\gamma }_s)\beta \varphi \wedge D_2(1-\psi )\\ -(\mu +{\gamma }_a)(\mu +d+{\gamma }_s)D_1\beta \varphi \wedge (1-\sigma )\\ -(\mu +{\gamma }_a)(\mu +d+{\gamma }_s)\beta \varphi \wedge (1-\sigma )(1-\psi )D_3\\ -(\mu +{\gamma }_a)(\mu +d+{\gamma }_s)\beta \varphi \wedge (1-\psi )(1-\sigma )D_6\\ -(\mu +d+{\gamma }_s)\beta {\eta }_1(1-v)\alpha \varphi D_2(1-\psi )\\ -(\mu +d+{\gamma }_s)\beta {\eta }_1(1-v)\alpha \varphi \wedge (1-\sigma )D_1\\ -(\mu +d+{\gamma }_s)\beta {\eta }_1(1-v)\alpha \varphi \wedge (1-\sigma )(1-\psi )D_3\\ -(\mu +d+{\gamma }_s)\beta {\eta }_1(1-v)\alpha \varphi \wedge (1-\psi )(1-\sigma )D_6\\ -(\mu +{\gamma }_a)\beta {\eta }_{2}v\alpha \varphi \wedge D_2(1-\psi )\\ -(\mu +{\gamma }_a)\beta {\eta }_{2}v\alpha \varphi \wedge (1-\sigma )D_1\\ -(\mu +{\gamma }_a)\beta {\eta }_{2}v\alpha \wedge (1-\sigma )(1-\psi )D_3\\ -(\mu +{\gamma }_a)\beta {\eta }_{2}v\alpha \varphi \wedge (1-\psi )(1-\sigma )D_6,\\ C_0=\mu D_9{D}_1{D}_2+\mu D_2{D}_6{D}_9+\mu D_3{D}_4{D}_9+\mu D_5{D}_9\\ -(\mu +{\gamma }_a)(\mu +d+{\gamma }_s)\beta \varphi \wedge D_1{D}_2\\ -(\mu +{\gamma }_a)(\mu +d+{\gamma }_s)\beta \varphi \wedge D_1{D}_3(1-\sigma )\\ -(\mu +{\gamma }_a)(\mu +d+{\gamma }_s)\beta \varphi \wedge (1-\psi )D_3{D}_7\\ -(\mu +{\gamma }_a)(\mu +d+{\gamma }_s)\beta \varphi \wedge (1-\psi )D_2{D}_6\\ -(\mu +d+{\gamma }_s)\beta {\eta }_1(1-v)\alpha \varphi \wedge D_1{D}_2\\ -(\mu +d+{\gamma }_s)\beta {\eta }_1(1-v)\alpha \varphi \wedge D_1{D}_3(1-\sigma )\\ -(\mu +d+{\gamma }_s)\beta {\eta }_1(1-v)\alpha \varphi \wedge (1-\psi )D_3{D}_7\\ -(\mu +d+{\gamma }_s)\beta {\eta }_1(1-v)\alpha \varphi \wedge (1-\psi )D_2{D}_6\\ -(\mu +{\gamma }_a)\beta {\eta }_{2}v\alpha \varphi \wedge D_1{D}_2\\ -(\mu +{\gamma }_a)\beta {\eta }_{2}v\alpha \varphi \wedge (1-\sigma )D_1{D}_3\\ -(\mu +{\gamma }_a)\beta {\eta }_{2}v\alpha \varphi \wedge (1-\psi )D_3{D}_7\\ -(\mu +{\gamma }_a)\beta {\eta }_{2}v\alpha \varphi \wedge (1-\psi )D_2{D}_6\\ \text{ which can be simplified to}\\ =D_9(K_1+K_2{D}_2)-(\mu +{\gamma }_a)(\mu +d+{\gamma }_s)\beta \varphi \wedge \varpi \\ -(\mu +d+{\gamma }_s)\beta {\eta }_1(1-v)\alpha \varphi \wedge \varpi \\ -(\mu +{\gamma }_a)\beta {\eta }_{2}v\alpha \varphi \wedge \varpi ,\\ =D_9(K_1+K_2{D}_2)(1-{\mathcal{R}}_c).\end{array}\hfill & (19)\end{array}$

In equation (18), λ^* = 0 corresponds to COVID-free equilibrium $\mathcal{CFE}.$ Hence, the possible number of endemic equilibria for model system (1) is obtained by evaluating for λ^* in the following third-degree polynomial

$\begin{array}{ll}{g}_1({\lambda }^{*})=C_3{\lambda }^{*3}+C_2{\lambda }^{*2}+C_1{\lambda }^{*}+C_0=0& (20)\end{array}$

and substituting the corresponding roots in (17). It is straightforward to note that if the vaccines administered are imperfect (i.e., σ, ψ < 1), then C₃ is always positive while C₀ can be positive/negative depending on whether the control reproduction number is less/greater than unity. Note that for ${\mathcal{R}}_c=1,$ C₀ reduces to zero. The number of endemic equilibria for the cubic polynomial (20) has been comprehensively studied using the technique by Descartes, usually referred as Descartes rule of signs [67]. Supplementary Table S1 (in the Supplementary Material S1) gives the possible number of COVID-19 persistent equilibria and the corresponding type of bifurcation structure.

For a perfect vaccine (i.e., σ = ψ = 1, ω₁ = ω₂ = 0), C₃ and C₂ become zero leading to the following linear equation:

$\begin{array}{ll}{g}_2({\lambda }^{*})=C_{11}{\lambda }^{**}+C_{00}=0,& (21)\end{array}$

where

$\begin{array}{l}{C}_{11}=\mu D_{22}{D}_9,C_{00}=D_9(K_{11}+K_{22}{D}_{22})(1-{\mathcal{R}}_{ca}),\\ {\mathcal{R}}_{ca}={\mathcal{R}}_c{|}_{(\sigma =\psi =1,{\omega }_1={\omega }_2=0)}=\frac{{\varpi }_1\beta \wedge \varphi Q_0}{(\mu +\varphi )(\mu +\alpha )(K_{11}+K_{22}{D}_{22})},\\ \text{where }{\varpi }_1=\mu D_{22},D_{22}=(\mu +(1-h_2){\theta }_2),\\ K_{11}=\mu (1-h_1){\theta }_1{D}_{22},K_{22}=\mu (\mu +(1-h_3){\theta }_3).\end{array}$

Solving linear equation (21) yields

$\begin{array}{l}{\lambda }^{**}=\frac{-D_9(K_{11}+K_{22}{D}_{22})(1-{\mathcal{R}}_{ca})}{C_{11}}.\end{array}$

Hence, equation (21) has a unique positive root λ^** whenever ${\mathcal{R}}_{ca}>1$ and a negative root whenever ${\mathcal{R}}_{ca}<1.$ The epidemiological implication of non-existence of a positive endemic equilibrium point when ${\mathcal{R}}_{ca}<1$ is that it is sufficient to eradicate COVID-19 pandemic by implementing control strategies that decrease control reproduction number below one. That is, for vaccines (both single-dose and double-dose vaccines) with 100% effectiveness against COVID-19, the model system (1) has a unique endemic equilibrium point which we shall denote by $E_1^{*}=(S^{**},V_p^{**},V^{**},E^{**}{I}_p^{**},I_a^{**},I_s^{**})$ that is globally asymptotically stable whenever ${\mathcal{R}}_{ca}>1.$ Thus, the following theorem follows.

Theorem 4. Provided the vaccine is perfect (i.e., σ = ψ = 1, ω₁ = ω₂ = 0) and ${\mathcal{R}}_{ca}>1$, the unique endemic equilibrium point $E_1^{*}$ is globally asymptotically stable.

Proof: Using the well-known scalar function for investigating global stability [68] [see also the recent study [69] where an auxiliary function has been used], we define the following Lyapunov candidate function:

$\begin{array}{l}\mathcal{W}\left(S,V_p,V,E,I_p,I_a,I_s\right)=\left(S-S^{**}-S^{**}ln\frac{S}{S^{**}}\right)\end{array}$

$\begin{array}{ll}+\left(V_p-V_p^{**}-V_p^{**}ln\frac{V_p}{V_p^{**}}\right)+\left(V-V^{**}-V^{**}ln\frac{V}{V^{**}}\right)& (22)\end{array}$

$\begin{array}{ll}+\left(E-E^{**}-E^{**}ln\frac{E}{E^{**}}\right)+\left(I_p-I_p^{**}-V_p^{**}ln\frac{V_p}{V_p^{**}}\right)& (23)\end{array}$

$\begin{array}{ll}+\left(I_a-I_a^{**}-I_a^{**}ln\frac{I_a}{I_a^{**}}\right)+\left(I_s-I_s^{**}-I_s^{**}ln\frac{I_s}{I_s^{**}}\right).& (24)\end{array}$

The orbital derivative of (22) is given by

$\begin{array}{l}\frac{d\mathcal{W}}{dt}=\left(1-\frac{S^{**}}{S}\right)\frac{dS}{dt}+\left(1-\frac{V_p^{**}}{V_p}\right)\frac{d{V}_p}{dt}\\ +\left(1-\frac{V^{**}}{V}\right)\frac{dV}{dt}+\left(1-\frac{E^{**}}{E}\right)\frac{dE}{dt}\\ +\left(1-\frac{I_p^{**}}{I_p}\right)\frac{d{I}_p}{dt}+\left(1-\frac{I_a^{**}}{I_a}\right)\frac{d{I}_a}{dt}+\left(1-\frac{I_s^{**}}{I_s}\right)\frac{d{I}_s}{dt},\\ =\left(1-\frac{S^{**}}{S}\right)(\wedge -(\mu +\lambda +(1-h_1){\theta }_1+(1-h_3){\theta }_3)S)\\ +\left(1-\frac{V_p^{**}}{V_p}\right)((1-h_1){\theta }_{1}S-(\mu +(1-h_2){\theta }_2)V_p)\\ +\left(1-\frac{V^{**}}{V}\right)((1-h_2){\theta }_2{V}_p+(1-h_3){\theta }_{3}S-\mu V)\\ +\left(1-\frac{E^{**}}{E}\right)(\lambda S-(\mu +\varphi )E)\\ +\left(1-\frac{I_p^{**}}{I_p}\right)(\varphi E-(\mu +\alpha )I_p)\\ +\left(1-\frac{I_a^{**}}{I_a}\right)((1-v)\alpha I_p-(\mu +{\gamma }_a)I_a)\\ +\left(1-\frac{I_s^{**}}{I_s}\right)(v\alpha I_p-(\mu +d+{\gamma }_s)I_s).& (25)\end{array}$

At equilibrium, the following equalities hold:

$\begin{array}{l}\wedge =(\mu +{\lambda }^{**}+(1-h_1){\theta }_1+(1-h_3){\theta }_3)S^{**},\\ (1-h_1){\theta }_1\frac{S^{**}}{V_p^{**}}=(\mu +(1-h_2){\theta }_2),\\ \mu =(1-h_2){\theta }_2\frac{V_p^{**}}{V^{**}}+(1-h_3){\theta }_3\frac{S^{**}}{V^{**}},\\ {\lambda }^{**}=(\mu +\varphi )\frac{E^{**}}{S^{**}},(\mu +\alpha )=\frac{\varphi E^{**}}{I_p^{**}},\\ (1-v)\alpha =(\mu +{\gamma }_a)\frac{I_a^{**}}{I_p^{**}},v\alpha =(\mu +d+{\gamma }_s)\frac{I_s^{**}}{I_p^{**}}.& (26)\end{array}$

Substituting (26) in (25) and simplifying yields

$\begin{array}{l}\frac{d\mathcal{W}}{dt}=\mu S^{**}\left(2-\frac{S}{S^{**}}-\frac{S^{**}}{S}\right)\\ +{\lambda }^{**}{S}^{**}\left(3-\frac{S^{**}}{S}-\frac{E}{E^{**}}-\frac{\lambda S{E}^{**}}{{\lambda }^{**}{S}^{**}E}\right)\\ +(1-h_1){\theta }_1{S}^{**}\left(3-\frac{S^{**}}{S}-\frac{V_p}{V_p^{**}}-\frac{S{V}_p^{**}}{S^{**}{V}_p}\right)\\ +(1-h_3){\theta }_3{S}^{**}\left(3-\frac{S^{**}}{S}-\frac{V}{V^{**}}-\frac{V^{**}S}{V{S}^{**}}\right)\\ +(1-h_2){\theta }_2{V}_p^{**}\left(1+\frac{V_p}{V_p^{**}}-\frac{V}{V^{**}}-\frac{V_p{V}^{**}}{V_p^{**}V}\right)\\ +\varphi E^{**}\left(1+\frac{E}{E^{**}}-\frac{I_p}{I_p^{**}}-\frac{I_p^{**}E}{I_p{E}^{**}}\right)\\ +(\mu +d+{\gamma }_s)I_s^{**}\left(1+\frac{I_p}{I_p^{**}}-\frac{I_s}{I_s^{**}}-\frac{I_s^{**}{I}_p}{I_s^{**}{I}_p^{**}}\right)\\ +(\mu +{\gamma }_a)I_a^{**}\left(1+\frac{I_p}{I_p^{**}}-\frac{I_a}{I_a^{**}}-\frac{I_a^{**}{I}_p}{I_a{I}_p^{**}}\right).& (27)\end{array}$

The fifth, sixth, seventh, and eighth terms in equation (27) together with further explanation in Supplementary Material [S2] can be re-written as follows:

$\begin{array}{l}(1-h_2){\theta }_2{V}_p^{**}\left(1+\frac{V_p}{V_p^{**}}-\frac{V}{V^{**}}-\frac{V_p{V}^{**}}{V_p^{**}V}\right)\\ =(1-h_2){\theta }_2{V}_p^{**}(2-\frac{V_p{V}^{**}}{V_p^{**}V}-\frac{V_p^{**}V}{V_p{V}^{**}}\\ +\left(\frac{V_p^{**}V}{V_p{V}^{**}}-1\right)\left(1-\frac{V_p}{V_p^{**}}\right))& (28)\end{array}$

$\begin{array}{ll}\le (1-h_2){\theta }_2{V}_p^{**}\left(2-\frac{V_p{V}^{**}}{V_p^{**}V}-\frac{V_p^{**}V}{V_p{V}^{**}}\right),& (29)\end{array}$

$\begin{array}{l}\varphi E^{**}\left(1+\frac{E}{E^{**}}-\frac{I_p}{I_p^{**}}-\frac{I_p^{**}E}{I_p{E}^{**}}\right)\\ =\varphi E^{**}\left(2-\frac{I_p^{**}E}{I_p{E}^{**}}-\frac{I_p{E}^{**}}{I_p^{**}E}+\left(\frac{I_p{E}^{**}}{I_p^{**}E}-1\right)\left(1-\frac{E}{E^{**}}\right)\right)\\ \le \varphi E^{**}\left(2-\frac{I_p^{**}E}{I_p{E}^{**}}-\frac{I_p{E}^{**}}{I_p^{**}E}\right),& (30)\end{array}$

$\begin{array}{l}(\mu +d+{\gamma }_s)I_s^{**}\left(1+\frac{I_p}{I_p^{**}}-\frac{I_s}{I_s^{**}}-\frac{I_s^{**}{I}_p}{I_s{I}_p^{**}}\right)\\ =(\mu +d+{\gamma }_s)I_s^{**}(2-\frac{I_s^{**}{I}_p}{I_s{I}_p^{**}}-\frac{I_s{I}_p^{**}}{I_s^{**}{I}_p}\\ +\left(\frac{I_s{I}_p^{**}}{I_s^{**}{I}_p}-1\right)\left(1-\frac{I_p}{I_p^{**}}\right))& (31)\end{array}$

$\begin{array}{ll}\le (\mu +d+{\gamma }_s)I_s^{**}\left(2-\frac{I{s}^{**}{I}_p}{I_s{I}_p^{**}}-\frac{I_s{I}_p^{**}}{I_s^{**}{I}_p}\right),& (32)\end{array}$

$\begin{array}{l}(\mu +{\gamma }_a)I_a^{**}\left(1+\frac{I_p}{I_p^{**}}-\frac{I_a}{I_a^{**}}-\frac{I_a^{**}{I}_p}{I_a{I}_p^{**}}\right)\\ =(\mu +{\gamma }_a)I_a^{**}\left(2-\frac{I_a^{**}{I}_p}{I_a{I}_p^{**}}-\frac{I_a{I}_p^{**}}{I_a^{**}{I}_p}+\left(\frac{I_a{I}_p^{**}}{I_a^{**}{I}_p}-1\right)\left(1-\frac{I_p}{I_p^{**}}\right)\right)\\ \le (\mu +{\gamma }_a)I_a^{**}\left(2-\frac{I_a^{**}{I}_p}{I_a{I}_p^{**}}-\frac{I_a{I}_p^{**}}{I_a^{**}{I}_p}\right).& (33)\end{array}$

Substituting inequalities (29)-(33) in (27) leads to

$\begin{array}{l}\frac{d\mathcal{W}}{dt}\le \mu S^{**}\left(2-\frac{S}{S^{**}}-\frac{S^{**}}{S}\right)\\ +{\lambda }^{**}{S}^{**}\left(3-\frac{S^{**}}{S}-\frac{E}{E^{**}}-\frac{\lambda S{E}^{**}}{{\lambda }^{**}{S}^{**}E}\right)\\ +(1-h_1){\theta }_1{S}^{**}\left(3-\frac{S^{**}}{S}-\frac{V_p}{V_p^{**}}-\frac{S{V}_p^{**}}{S^{**}{V}_p}\right)\end{array}$

$\begin{array}{ll}+(1-h_3){\theta }_3{S}^{**}\left(3-\frac{S^{**}}{S}-\frac{V}{V^{**}}-\frac{V^{**}S}{V{S}^{**}}\right)& (34)\end{array}$

$\begin{array}{l}+(1-h_2){\theta }_2{V}_p^{**}\left(2-\frac{V_p{V}^{**}}{V_p^{**}V}-\frac{V_p^{**}V}{V_p{V}^{**}}\right)\\ +\varphi E^{**}\left(2-\frac{I_p^{**}E}{I_p{E}^{**}}-\frac{I_p{E}^{**}}{I_p^{**}E}\right)\\ +(\mu +d+{\gamma }_s)I_s^{**}\left(2-\frac{I_s^{**}{I}_p}{I_s{I}_p^{**}}-\frac{I_s{I}_p^{**}}{I_s^{**}{I}_p}\right)\\ +(\mu +{\gamma }_a)I_a^{**}\left(2-\frac{I_a^{**}{I}_p}{I_a{I}_p^{**}}-\frac{I_a{I}_p^{**}}{I_a^{**}{I}_p}\right).\end{array}$

Considering the fact that the geometric mean is smaller than the arithmetic mean, it follows that $\frac{d\mathcal{W}}{dt}\le 0,$ with equality holding (i.e., $\frac{d\mathcal{W}}{dt}=0$) if and only if $S=S^{**},V_p=V_p^{**},V=V^{**},E=E^{**},I_p=I_p^{**},I_a=I_a^{**},$ and $I_s=I_s^{**}.$ Hence, $\mathcal{W}$ is a suitable Lyapunov function in the region Δ, and the unique endemic equilibrium $E_1^{*}$ is globally asymptotically stable when ${\mathcal{R}}_{ca}>1.$ This implies that if the vaccines being administered are perfect, all trajectories of model system (1) whose initial conditions are in the region Δ eventually converge to the unique endemic equilibrium point $E_1^{*}$ as t → ∞ provided ${\mathcal{R}}_{ca}>1.$ Epidemiologically, it means that model system (1) cannot exhibit the phenomenon of backward bifurcation if the vaccine efficacy is 100% and everlasting.

3.6 Analysis of non-specific case (imperfect vaccine)

Theorem 4 implies that in case the vaccines administered are perfect and livelong, then the model system (1) cannot have the phenomenon of backward bifurcation. In this subsection, we relax the possibility of a perfect and everlasting vaccine but instead investigate model dynamics in presence of an imperfect vaccine.

3.6.1 Occurence of the phenomenon of backward bifurcation

To investigate whether the model system (1) exhibits bistability phenomenon, we apply the center manifold theory as descrbed by Castillo-Chavez et al. [70]. Given at the point ${\mathcal{R}}_{ca}=1$ is where there is change in bifurcation structure, we choose β = β^* as our bifurcation parameter where β^* is obtained by setting ${\mathcal{R}}_{ca}=1$ and evaluating for β. That is

$\begin{array}{ll}\beta ={\beta }^{*}=\frac{(\mu +\varphi )(\mu +\alpha )(\mu +{\gamma }_a)(\mu +d+{\gamma }_s)(K_1+K_2{D}_2)}{\varpi \wedge \varphi Q_0}.& (35)\end{array}$

For the purpose of simplification, we re-write model system (1) using new variables denoted by S = x₁, V_p = x₂, V = x₃, E = x₄, I_p = x₅, I_a = x₆, I_s = x₇) and the total population size $N=\sum _i^7{x}_i.$ Let $X={(x_1,x_2,x_3,x_4,x_5,x_6,x_7)}^T$ (T denotes transpose), which implies $\frac{dX}{dt}=F(X)$ where F = (f₁, f₂, f₃, f₄, f₅, f₆, f₇). With the new variables, the $\mathcal{CFE}=(x_1^{*},x_2^{*},x_3^{*},0,0,0,0).$ Now, we define the model system (1) as follows:

$\begin{array}{l}\frac{d{x}_1}{dt}=f_1=\wedge +{\omega }_1{x}_2+{\omega }_2{x}_3-(\mu +\beta (x_5+{\eta }_1{x}_6+{\eta }_2{x}_7)\\ +(1-h_1){\theta }_1+(1-h_3){\theta }_3)x_1,\\ \frac{d{x}_2}{dt}=f_2=(1-h_1){\theta }_1{x}_1-(\mu +{\omega }_1+(1-h_2){\theta }_2\\ +(1-\sigma )\beta (x_5+{\eta }_1{x}_6+{\eta }_2{x}_7))x_2,\\ \frac{d{x}_3}{dt}=f_3=(1-h_2){\theta }_2{x}_2+(1-h_3){\theta }_3{x}_1-(\mu +{\omega }_2\\ +(1-\psi )\beta (x_5+{\eta }_1{x}_6+{\eta }_2{x}_7))x_3,& (36)\end{array}$

$\begin{array}{l}\frac{d{x}_4}{dt}=f_4=\beta (x_5+{\eta }_1{x}_6+{\eta }_2{x}_7)x_4\\ +(1-\sigma )\beta (x_5+{\eta }_1{x}_6+{\eta }_2{x}_7)x_2\\ +(1-\psi )\beta (x_5+{\eta }_1{x}_6+{\eta }_2{x}_7)x_3-(\mu +\varphi )x_4,\\ \frac{d{x}_5}{dt}=f_5=\varphi x_4-(\mu +\alpha )x_5,\\ \frac{d{x}_6}{dt}=f_6=(1-v)\alpha x_5-(\mu +{\gamma }_a)x_6,\\ \frac{d{x}_7}{dt}=f_7=v\alpha x_5-(\mu +d+{\gamma }_s)x_7.\end{array}$

The Jacobian matrix of model (36) evaluated at $\mathcal{CFE}$ is given by:

$J(\mathcal{CFE}){|}_{\beta ={\beta }^{*}}=\left(\begin{array}{ccccccc}-(\mu +D_3+D_6)& {\omega }_1& {\omega }_2& 0& -{\beta }^{*}{x}_1^{*}& -{\beta }^{*}{\eta }_1{x}_1^{*}& -\beta {\eta }_2{x}_1^{*}\\ D_7& -D_2& 0& 0& -(1-\sigma ){\beta }^{*}{x}_2^{*}& -(1-\sigma ){\beta }^{*}{\eta }_1{x}_2^{*}& -(1-\sigma ){\beta }^{*}{\eta }_2{x}_2^{*}\\ D_6& D_7& -D_1& 0& -(1-\psi ){\beta }^{*}{x}_3^{*}& -(1-\psi ){\beta }^{*}{\eta }_1{x}_3& -(1-\psi ){\beta }^{*}{\eta }_2{x}_3\\ 0& 0& 0& -(\mu +\varphi )& J_{45}& J_{46}& J_{47}\\ 0& 0& 0& \varphi & -(\mu +\alpha )& 0& 0\\ 0& 0& 0& 0& -(1-v)\alpha & -(\mu +{\gamma }_a)& 0\\ 0& 0& 0& 0& v\alpha & 0& -(\mu +d+{\gamma }_s)\end{array}\right),$

where $J_{45}={\beta }^{*}{x}_1^{*}+(1-\sigma ){\beta }^{*}{x}_2^{*}+(1-\psi ){\beta }^{*}{x}_3^{*},J_{46}={\beta }^{*}{\eta }_1{x}_1^{*}+(1-\sigma ){\beta }^{*}{\eta }_1{x}_2^{*}+(1-\psi ){\beta }^{*}{\eta }_1{x}_3^{*},$ and $J_{47}={\beta }^{*}{\eta }_2{x}_1^{*}+(1-\sigma ){\beta }^{*}{\eta }_2{x}_2^{*}+(1-\psi ){\beta }^{*}{\eta }_2{x}_3^{*}$

Now, the right eigenvector corresponding to the zero eigenvalue of the Jacobian matrix $J(\mathcal{CFE}){|}_{\beta ={\beta }^{*}}$ is obtained as $Ṽ={(\widetilde{v_1},\widetilde{v_2},\widetilde{v_3},\widetilde{v_4},\widetilde{v_5},\widetilde{v_6},\widetilde{v_7})}^T,$ where

$\begin{array}{l}{ṽ}_1=\left(\frac{D_2(P_1{D}_1+{\omega }_2{P}_2)}{D_7(Q_1{D}_1+{\omega }_2{Q}_2)}+\frac{(1-\sigma )x_2^{*}{Q}_0}{D_7(\mu +{\gamma }_a)(\mu +d+{\gamma }_s)}\right){\beta }^{*}\widetilde{v_5},\\ {ṽ}_2=\left(\frac{P_1{D}_1+{\omega }_2{P}_2}{Q_1{D}_1+{\omega }_2{Q}_2}\right){\beta }^{*}\widetilde{v_5},\\ \widetilde{v_3}=\left(\frac{(Q_2{P}_1{D}_1+Q_2{\omega }_2{P}_2)-P_2(Q_1{D}_1+{\omega }_2{Q}_2)}{D_1{D}_7(Q_1{D}_1+{\omega }_2{Q}_2)}\right){\beta }^{*}\widetilde{v_5},\\ \widetilde{v_4}=\frac{(\mu +\alpha )}{\varphi }\widetilde{v_5},\\ \widetilde{v_5}=\widetilde{v_5}>0,\\ \widetilde{v_6}=\frac{(1-v)\alpha }{(\mu +{\gamma }_a)}\widetilde{v_5},\\ {ṽ}_7=\frac{v\alpha }{(\mu +d+{\gamma }_s)}{ṽ}_5,\text{and}\end{array}$

$\begin{array}{l}{P}_1=\frac{(\mu +D_3+D_6)(1-\sigma )x_2^{*}{Q}_0+x_1^{*}{Q}_0{D}_7}{(\mu +{\gamma }_a)(\mu +d+{\gamma }_s)},\\ P_2=\frac{(1-\psi )x_3^{*}{D}_7{Q}_0-(1-\sigma )x_2^{*}{D}_6{Q}_0}{(\mu +{\gamma }_a)(\mu +d+{\gamma }_s)},\\ Q_1=({\omega }_1{D}_7-(\mu +D_3+D_6)D_2),Q_2=D_2{D}_6+D_7^2.\end{array}$

Similarly, the left eigenvector corresponding to the zero eigenvalue of $J(\mathcal{CFE}){|}_{\beta ={\beta }^{*}}$ is given by $\tilde{W}={({\tilde{w}}_1,{\tilde{w}}_2,{\tilde{w}}_3,{\tilde{w}}_4,{\tilde{w}}_5,{\tilde{w}}_6,{\tilde{w}}_7)}^T$ where

$\begin{array}{l}{\tilde{w}}_1={\tilde{w}}_2={\tilde{w}}_3=0,{\tilde{w}}_4=\frac{\varphi }{(\mu +\varphi )}{w}_5,{\tilde{w}}_5={\tilde{w}}_5>0,\\ {\tilde{w}}_6=\frac{J_{46}\varphi }{(\mu +\varphi )(\mu +\gamma )}{\tilde{w}}_5,\\ {\tilde{w}}_7=\frac{J_{47}\varphi }{(\mu +\varphi )(\mu +d+{\gamma }_s)}{\tilde{w}}_5.\end{array}$

As described in Theorem 4.1 of Castillo-Chavez [70], the type of bifurcation structure exhibited by the model is determined by the bifurcation coefficients, a and b which are defined as

$\begin{array}{l}a={\displaystyle \sum _{k,i,j=1}^7}{\tilde{w}}_k{ṽ}_i{ṽ}_j\frac{{\partial }^2{f}_k(0,0)}{\partial x_i\partial x_j},b={\displaystyle \sum _{k,j=1}^7}{\tilde{w}}_k{ṽ}_i\frac{{\partial }^2{f}_k(0,0)}{\partial x_i\partial {\beta }^{*}}.\end{array}$

Computations involving parameter a: Now, we proceed and determine the non-vanishing partial derivatives evaluated at the COVID-free equilibrium $\mathcal{CFE}$ of model equation (36). These partial derivatives include:

$\begin{array}{l}\frac{{\partial }^2{f}_4(0,0)}{\partial x_5\partial x_1}={\beta }^{*},\frac{{\partial }^2{f}_4(0,0)}{\partial x_6\partial x_1}={\eta }_1{\beta }^{*},\frac{{\partial }^2{f}_4(0,0)}{\partial x_7\partial x_1}\\ ={\eta }_2{\beta }^{*},\frac{{\partial }^2{f}_4(0,0)}{\partial x_5\partial x_2}=(1-\sigma ){\beta }^{*},\frac{{\partial }^2{f}_4(0,0)}{\partial x_6\partial x_2}\\ =(1-\sigma ){\beta }^{*}{\eta }_1,\frac{{\partial }^2{f}_4(0,0)}{\partial x_7\partial x_3}=(1-\sigma ){\beta }^{*}{\eta }_2,\frac{{\partial }^2{f}_4(0,0)}{\partial x_5\partial x_3}\\ =(1-\psi ){\beta }^{*},\frac{{\partial }^2{f}_4(0,0)}{\partial x_6\partial x_3}=(1-\psi ){\beta }^{*}{\eta }_1,\frac{{\partial }^2{f}_4(0,0)}{\partial x_7\partial x_3}=(1-\psi ){\beta }^{*}{\eta }_2.\end{array}$

Hence,

$\begin{array}{l}a={\tilde{w}}_4{ṽ}_1{ṽ}_2\frac{{\partial }^2{f}_4(0,0)}{\partial x_1\partial x_5}+{\tilde{w}}_4{ṽ}_1{ṽ}_6\frac{{\partial }^2{f}_4(0,0)}{\partial x_1\partial x_6}+{\tilde{w}}_4{ṽ}_1{ṽ}_7\frac{{\partial }^2{f}_4(0,0)}{\partial x_1\partial x_7}\\ +{\tilde{w}}_4{ṽ}_2{ṽ}_5\frac{{\partial }^2{f}_4(0,0)}{\partial x_2\partial x_5}+{\tilde{w}}_4{ṽ}_2{ṽ}_6\frac{{\partial }^2{f}_4(0,0)}{\partial x_2\partial x_6}\\ +{\tilde{w}}_4{ṽ}_2{ṽ}_7\frac{{\partial }^2{f}_4(0,0)}{\partial x_2\partial x_7}+{\tilde{w}}_4{ṽ}_3{ṽ}_5\frac{{\partial }^2{f}_4(0,0)}{\partial x_3\partial x_5}+{\tilde{w}}_4{ṽ}_3{ṽ}_6\frac{{\partial }^2{f}_4(0,0)}{\partial x_3\partial x_6}\\ +{\tilde{w}}_4{ṽ}_3{ṽ}_7\frac{{\partial }^2{f}_4(0,0)}{\partial x_3\partial x_7}\\ \begin{array}{c}=\frac{\varphi {\tilde{w}}_5{ṽ}_5^2{\beta }^{*2}{Q}_0}{(\mu +\varphi )(\mu +{\gamma }_a)(\mu +d+{\gamma }_s)}({\Gamma }_1-{\Gamma }_2),\text{where}\\ {\Gamma }_1=\left(\frac{D_2(P_1{D}_1+{\omega }_2{P}_2)}{D_7(Q_1{D}_1+{\omega }_2{Q}_2)}+\frac{(1-\sigma )x_2^{*}{Q}_0}{D_7(\mu +{\gamma }_a)(\mu +d+{\gamma }_s)}\right)\hfill \\ +\frac{(P_1{D}_1+{\omega }_2{P}_2)(1-\sigma )}{(Q_1{D}_1+{\omega }_2{Q}_2)}\hfill \\ +\frac{(Q_2{P}_1{D}_1+Q_2{\omega }_2{P}_2)(1-\psi )}{D_1{D}_7(Q_1{D}_1+{\omega }_2{Q}_2)}\hfill \\ {\Gamma }_2=\frac{P_2}{D_1{D}_7}.\hfill \end{array}\end{array}$

Computations involving bifurcation parameter b: The bifurcation parameter b is associated with the following non-vanishing partial derivatives evaluated at $\mathcal{CFE}:$

$\begin{array}{l}\frac{{\partial }^2{f}_4(0,0)}{\partial x_5\partial {\beta }^{*}}=\frac{\begin{array}{l}\wedge (D_1{D}_2+(1-\sigma )(1-h_1){\theta }_1{D}_1+(1-\psi )\hfill \\ ((1-h_1){\theta }_1(1-h_2){\theta }_2+(1-h_3){\theta }_3{D}_2))\hfill \end{array}}{K_1+K_2{D}_2}\\ =\frac{\wedge \varpi }{K_1+K_2{D}_2},\\ \frac{{\partial }^2{f}_4(0,0)}{\partial x_6\partial {\beta }^{*}}=\frac{\begin{array}{l}{\eta }_1\wedge (D_1{D}_2+(1-\sigma )(1-h_1){\theta }_1{D}_1\hfill \\ +(1-\psi )((1-h_1){\theta }_1(1-h_2){\theta }_2+(1-h_3){\theta }_3{D}_2))\hfill \end{array}}{K_1+K_2{D}_2}\\ =\frac{{\eta }_1\wedge \varpi }{K_1+K_2{D}_2},\\ \frac{{\partial }^2{f}_4(0,0)}{\partial x_7\partial {\beta }^{*}}=\frac{\begin{array}{l}{\eta }_2\wedge (D_1{D}_2+(1-\sigma )(1-h_1){\theta }_1{D}_1\hfill \\ +(1-\psi )((1-h_1){\theta }_1(1-h_2){\theta }_2+(1-h_3){\theta }_3{D}_2))\hfill \end{array}}{K_1+K_2{D}_2}\\ =\frac{{\eta }_2\wedge \varpi }{K_1+K_2{D}_2}.\end{array}$

Hence,

$\begin{array}{l}b={\tilde{w}}_4{ṽ}_5\frac{{\partial }^2{f}_4(0,0)}{\partial x_5\partial {\beta }^{*}}+{\tilde{w}}_4{ṽ}_6\frac{{\partial }^2{f}_4(0,0)}{\partial x_6\partial {\beta }^{*}}+{\tilde{w}}_4{ṽ}_7\frac{{\partial }^2{f}_4(0,0)}{\partial x_7\partial {\beta }^{*}}\\ =\left(\frac{\varphi {\tilde{w}}_5{ṽ}_5}{(\mu +\varphi )}\right)\left(\frac{\wedge \varpi }{K_1+K_2{D}_2}\right)\left(\frac{Q_0}{(\mu +{\gamma }_a)(\mu +d+{\gamma }_s)}\right)>0.\end{array}$

Note that the eigenvectors ${\tilde{w}}_4$ and ṽ₅ are always positive and Q₀ > 0; hence, b is non-negative. For backward bifurcation to occur, we need both a and b to be positive. Notice that the bifurcation parameter a is positive if Γ₁ > Γ₂. Thus, backward bifurcation can occur if the inequality

$\begin{array}{ll}{\Gamma }_1>{\Gamma }_2& (37)\end{array}$

holds. In case vaccines are perfect and permanent, (i.e., σ = ψ = 1, ω₁ = ω₂ = 0) Γ₂, P₂ are equal to zero while Γ₁ reduces to

$\begin{array}{l}{\Gamma }_1=-\frac{D_{22}{Q}_0{x}_{11}^{*}}{\begin{array}{l}(\mu +(1-h_2){\theta }_2)((1-h_3){\theta }_3+(\mu \hfill \\ +(1-h_1){\theta }_1))(\mu +{\gamma }_a)(\mu +d+{\gamma }_s)\hfill \end{array}}<0,\end{array}$

where

$\begin{array}{l}{x}_{11}^{*}=S^0{|}_{(\psi =\sigma =1,{\omega }_1={\omega }_2=0)}=\frac{\mu D_{11}}{(K_{11}+K_{22}{D}_{22})}>0.\end{array}$

The observations that a < 0 when vaccines administered are assumed to be perfect and everlasting (i.e., ψ = σ = 1, ω₁ = ω₂ = 0) agree with the global asymptotic stability results proved in Theorem 4. However, if the vaccines administered among susceptible and partially vaccinated cohorts are imperfect and non-permanent hence, allowing for breakthrough infections to occur, the well-known backward bifurcation phenomenon [71–73] may arise. The epidemiological implication of the occurrence of backward bifurcation is that it will be more difficult for implemented intervention measures to eliminate the COVID-19 pandemic. This is due to the fact that reducing ${\mathcal{R}}_c$ below unity will be necessary but not sufficient. Thus, following the above bifurcation analysis, we state the following result:

Theorem 5. Provided the inequality given by (37) holds, the COVID-19 model system (1) exhibits backward bifurcation phenomenon when ${\mathcal{R}}_c$ crosses unity.

4 Model incorporating optimal control theory

To gain insight on the most effective combination of mitigation strategies that can be adopted to significantly minimize the spread of COVID-19 pandemic as well as the cost, we introduce three time-dependent control measures that are broadly categorized as either non-pharmaceutical (NPIs) and pharmaceutical control measures. Thus, the proposed model (1) is transformed to an optimal control COVID-19 infection model. The non-pharmaceutical time-dependent COVID-19 control is denoted by u₁(t), which represents preventive control measures that involve social/physical distancing, wearing face mask, hand washing, and adherent to education on proper hygiene among susceptible, partially vaccinated (with a double-dose vaccine) and fully vaccinated individuals (with either double-dose or single-dose vaccine). The pharmaceutical time-dependent intervention measures include u₂(t) which represents a time-dependent screening-management control measure for infectious individuals in the subgroups I_a and I_s, and θ₃(t) which represents time-dependent vaccination control for susceptible individuals with a single-dose vaccine. Incorporating these time-dependent control measures, the model system (1) transforms to the following non-linear (also non-autonomous) system of ordinary differential equations:

$\begin{array}{l}\frac{dS}{dt}=\wedge +{\omega }_1{V}_p+{\omega }_{2}V-(\mu +(1-u_1(t))\lambda +(1-h_1){\theta }_1\\ +(1-h_3){\theta }_3(t))S,\frac{d{V}_p}{dt}=(1-h_1){\theta }_{1}S-(\mu +{\omega }_1\\ +(1-h_2){\theta }_2+(1-u_1(t))(1-\sigma )\lambda )V_p,\\ \frac{dV}{dt}=(1-h_2){\theta }_2{V}_p+(1-h_3){\theta }_3(t)S-(\mu +{\omega }_2\\ +(1-u_1(t))(1-\psi )\lambda )V,\end{array}$

$\begin{array}{ll}\frac{dE}{dt}=(1-u_1(t))\lambda (S+(1-\sigma )V_p+(1-\psi )V)-(\mu +\varphi )E,& (38)\end{array}$

$\begin{array}{l}\frac{d{I}_p}{dt}=\varphi E-(\mu +\alpha )I_p,\\ \frac{d{I}_a}{dt}=(1-v)\alpha I_p-(\mu +{\gamma }_a+c{u}_2(t))I_a,\\ \frac{d{I}_s}{dt}=v\alpha I_p-(\mu +d+{\gamma }_s+c{u}_2(t))I_s,\\ \frac{dR}{dt}=({\gamma }_a+c{u}_2(t))I_a+({\gamma }_s+c{u}_2(t))I_s-\mu R.\end{array}$

In model system (38), the parameter c accounts for control rate constant while the terms 1 − u₁(t), θ₃(t), and u₂(t) are the controls that aim to curb the spread of the COVID-19 pandemic over a long time scale. Note that u₁, θ₃, u₂ ∈ [0, 1]. Given the objective of time-dependent control measures is to attempt to mitigate the proliferation of COVID-19 pandemic and eventually eradicating it by; minimizing the sub-populations of infectious individuals and also minimizing the cost of mitigation strategies, we define the objective functional as

$\begin{array}{r}\mathcal{J}(u_1,{\theta }_3,u_2)={\displaystyle {\int }_0^{T_f}({\mathcal{A}}_1{I}_p(t)+{\mathcal{A}}_2{I}_a(t)+{\mathcal{A}}_3{I}_s(t)}\\ +\frac{1}{2}\left({ℬ}_1{u}_1^2(t)+{ℬ}_2{\theta }_3^2(t)+{ℬ}_3{u}_2^2(t)\right))dt,& (39)\end{array}$

where T_f denotes the terminal time such that t ∈ [0, T_f] while ${\mathcal{A}}_j$ and ${\mathcal{B}}_j$ (j = 1, 2, 3) represent weight constants. The term ${\mathcal{B}}_1{u}_1^2/2$ is the cost associated with adopting non-pharmaceutical intervention measures which include social/physical distancing, wearing face masks, hand washing, and adhering to proper hygiene education among individuals in the subgroups S, V_p, and V, and ${\mathcal{B}}_2{\theta }_3^2/2$ is the cost associated with vaccination of susceptible individuals with a single-dose vaccine. Furthermore, the term ${\mathcal{B}}_3{u}_2^2/2$ is the cost associated with COVID-19 screening-management for both individuals in I_a and I_s compartments. The choice of a quadratic cost function is motivated by the relevant literature in the sequel on optimal control problems [for instance, see Okyere et al. [74], Ghosh et al. [75], Bonyah et al. [76], Purwati et al. [77], Lenhart and Workman [78] and Olaniyi et al. [79]]. Thus, we now proceed to find a ${\mathcal{U}}^{*}=(u_1^{*},{\theta }_3^{*},u_2^{*})$ which fulfills

$\begin{array}{ll}\mathcal{J}({\mathcal{U}}^{*})=min\left\{\mathcal{J}(u_1,{\theta }_3,u_2)\mid u_1,{\theta }_3,u_2\in \mathcal{U}\right\},& (40)\end{array}$

where $\mathcal{U}$ is a Lebesgue measurable control set which has a lower bound as zero and an upper bound 1 for t ∈ [0, T_f].

4.1 Characterization of the optimal controls

Following Pontryagin's maximum principle originally described in Pontryagin [80], we establish the necessary conditions that an optimal control COVID-19 model system (38) needs to fulfill. The technique outlined by Pontryagin et al. [80] enables conversion of the model system (38) in synergy with the objective functional (39) into a problem of minimizing pointwise a Hamiltonian $\mathcal{H}$ with respect to the time-dependent control measures; u₁(t), θ₃(t), u₂(t). Consequently, we define the Hamiltonian of model system (38) as

$\begin{array}{ll}\mathcal{H}(t,y,\mathcal{U})=\mathcal{P}(t,y,\mathcal{U})+{\displaystyle \sum _{j=1}^8}{\lambda }_j{q}_j,& (41)\end{array}$

where y = (S, V_p, V, E, I_p, I_a, I_s, R) represents COVID-19 model state variables, $\mathcal{P}(t,y,\mathcal{U})$ is a Lagrangian which represents the integrand of the objective functional, and λ_j represents the corresponding adjoint variables for the model states S, V_p, V, E, I_p, I_a, I_s, R. Moreover, q_j is the right-hand term of the model system (38). Thus, the Hamiltonian $\mathcal{H}$ is stated explicitly as

$\begin{array}{l}\mathscr{H}={\mathcal{A}}_1{I}_p+{\mathcal{A}}_2{I}_a+{\mathcal{A}}_3{I}_s+\frac{1}{2}({ℬ}_1{u}_1^2+{ℬ}_2{\theta }_3^2+{ℬ}_3{u}_2^2)\\ +{\lambda }_1[\wedge +{\omega }_1{V}_p+{\omega }_{2}V-(\mu +(1-u_1)\beta (I_p\\ +{\eta }_1{I}_a+{\eta }_2{I}_s)+(1-h_1){\theta }_1+(1-h_3){\theta }_3)S]\\ +{\lambda }_2[(1-h_1){\theta }_{1}S-(\mu +{\omega }_1+(1-h_2){\theta }_2\\ +(1-u_1)(1-\sigma )\beta (I_p+{\eta }_1{I}_a+{\eta }_2{I}_s))V_p]\\ +{\lambda }_3[(1-h_2){\theta }_2{V}_p+(1-h_3){\theta }_{3}S-(\mu +{\omega }_2\\ +(1-u_1)(1-\psi )\beta (I_p+{\eta }_1{I}_a+{\eta }_2{I}_s))V]\\ +{\lambda }_4[(1-u_1)(\beta (I_p+{\eta }_1{I}_a+{\eta }_2{I}_s)(S\\ +(1-\sigma )V_p+(1-\psi )V)-(\mu +\varphi )E]\\ +{\lambda }_5[\varphi E-(\mu +\alpha )I_p]\\ +{\lambda }_6[(1-v)\alpha I_p-(\mu +{\gamma }_a+c{u}_2)I_a]\\ +{\lambda }_7[v\alpha I_p-(\mu +d+{\gamma }_s+c{u}_2)I_s]\\ +{\lambda }_8[({\gamma }_a+c{u}_2)I_a+({\gamma }_s+c{u}_2)I_s-\mu R].\end{array}$

We now establish the following results:

Theorem 6. Assuming a solution of the optimal control problem for the objective functional $\mathcal{J}$ over the control set V is obtained, then the adjoint variables λ₁, λ₂, ⋯  , λ₈ fulfill the adjoint equations

$\begin{array}{l}\frac{d{\lambda }_j}{dt}=-\frac{\partial \mathcal{H}}{\partial y}\end{array}$

with λ_j(T_f) = 0 as transversality condition and y = (S, V_p, V, E, I_p, I_a, I_s, R). Moreover, optimal controls are given as

$\begin{array}{l}{u}_1^{*}=min\{1,max\{0,{\Delta }_1\}\},\\ {\theta }^{*}=min\{1,max\{0,{\Delta }_2\}\},\\ u_2^{*}=min\{1,max\{0,{\Delta }_3\}\},\end{array}$

where

$\begin{array}{l}{\Delta }_1=\frac{\begin{array}{c}\left[({\lambda }_4-{\lambda }_1)S+(1-\sigma )({\lambda }_4-{\lambda }_2)V_p+(1-\psi )({\lambda }_4-{\lambda }_3)V\right]\\ \beta (I_p+{\eta }_1{I}_a+{\eta }_2{I}_s)\end{array}}{{\mathcal{B}}_1},\\ {\Delta }_2=\frac{({\lambda }_1-{\lambda }_3)(1-h_3)S}{{\mathcal{B}}_2},\\ {\Delta }_2=\frac{({\lambda }_6-{\lambda }_8)c{I}_a+({\lambda }_7-{\lambda }_8)c{I}_s}{{\mathcal{B}}_3}.\end{array}$

Proof: Let ${\mathcal{U}}^{*}=(u_1^{*},{\theta }^{*},u_2^{*}).$ By Pontryagin's maximum principle, the state variables of the model (38) and adjoint variables satisfy the following relations which are obtained from the Hamiltonian $\mathcal{H}:$

$\begin{array}{l}\frac{d{\lambda }_1}{dt}=-\frac{\partial \mathcal{H}}{\partial S},{\lambda }_1(T_f)=0,\frac{d{\lambda }_2}{dt}=-\frac{\partial \mathcal{H}}{\partial V_p},{\lambda }_2(T_f)=0,\\ \frac{d{\lambda }_3}{dt}=-\frac{\partial \mathcal{H}}{\partial V},{\lambda }_3(T_f)=0,\frac{d{\lambda }_4}{dt}=-\frac{\partial \mathcal{H}}{\partial E},{\lambda }_4(T_f)=0,\\ \frac{d{\lambda }_5}{dt}=-\frac{\partial \mathcal{H}}{\partial I_p},{\lambda }_5(T_f)=0,\frac{d{\lambda }_6}{dt}=-\frac{\partial \mathcal{H}}{\partial I_a},{\lambda }_6(T_f)=0,\\ \frac{d{\lambda }_7}{dt}=-\frac{\partial \mathcal{H}}{\partial I_s},{\lambda }_7(T_f)=0,\frac{d{\lambda }_8}{dt}=-\frac{\partial \mathcal{H}}{\partial R},{\lambda }_8(T_f)=0.\end{array}$

That is, the adjoint equations can now be expressed as

$\begin{array}{ll}\frac{d{\lambda }_1}{dt}=({\lambda }_1-{\lambda }_4)(1-u_1)\beta (I_p+{\eta }_1{I}_a+{\eta }_2{I}_s)& (42)\end{array}$

$\begin{array}{l}+({\lambda }_1-{\lambda }_2)(1-h_1){\theta }_1+({\lambda }_1-{\lambda }_3)(1-h_3){\theta }_3+\mu {\lambda }_1,\end{array}$

$\begin{array}{ll}\frac{d{\lambda }_2}{dt}=({\lambda }_2-{\lambda }_1){\omega }_1+({\lambda }_2-{\lambda }_3)(1-h_2){\theta }_2& (43)\end{array}$

$\begin{array}{l}++({\lambda }_2-{\lambda }_4)(1-u_1)(1-\sigma )\beta (I_p+{\eta }_1{I}_a+{\eta }_2{I}_s)+\mu {\lambda }_2,\end{array}$

$\begin{array}{ll}\frac{d{\lambda }_3}{dt}=({\lambda }_3-{\lambda }_1){\omega }_2+({\lambda }_3-{\lambda }_4)(1-u_1)(1-\psi )& (44)\end{array}$

$\begin{array}{l}+\beta (I_p+{\eta }_1{I}_a+{\eta }_2{I}_s)+\mu {\lambda }_3,\end{array}$

$\begin{array}{ll}\frac{d{\lambda }_4}{dt}=({\lambda }_4-{\lambda }_5)\varphi +\mu {\lambda }_4,& (45)\end{array}$

$\begin{array}{ll}\frac{d{\lambda }_5}{dt}=-{\mathcal{A}}_1+({\lambda }_1-{\lambda }_4)(1-u_1)\beta S& (46)\end{array}$

$\begin{array}{l}+({\lambda }_2-{\lambda }_4)(1-u_1)(1-\sigma )\beta V_p\\ +({\lambda }_3-{\lambda }_4)(1-u_1)(1-\psi )\beta V+({\lambda }_5-{\lambda }_6)\alpha \\ +({\lambda }_6-{\lambda }_7)v\alpha +\mu {\lambda }_5,\\ \frac{d{\lambda }_6}{dt}=-{\mathcal{A}}_2+({\lambda }_1-{\lambda }_4)(1-u_1)\beta {\eta }_{1}S\\ +({\lambda }_2-{\lambda }_4)(1-u_1)(1-\sigma )\beta {\eta }_1{V}_p\\ +({\lambda }_3-{\lambda }_4)(1-u_1)(1-\psi )\beta {\eta }_{1}V\\ +({\lambda }_6-{\lambda }_8)({\gamma }_a+c{u}_2)+\mu {\lambda }_6,\\ \frac{d{\lambda }_7}{dt}=-{\mathcal{A}}_3+({\lambda }_1-{\lambda }_4)(1-u_1)\beta {\eta }_{2}S\\ +({\lambda }_2-{\lambda }_4)(1-u_1)(1-\sigma ){\beta }_2{V}_p\\ +({\lambda }_3-{\lambda }_4)(1-u_1)(1-\psi )\beta {\eta }_{2}V\\ +({\lambda }_7-{\lambda }_8)({\gamma }_s+c{u}_2)+(\mu +d){\lambda }_7,\\ \frac{d{\lambda }_8}{dt}=\mu {\lambda }_8,\end{array}$

with the transversality conditions:

$\begin{array}{r}{\lambda }_1(T_f)={\lambda }_2(T_f)={\lambda }_3(T_f)={\lambda }_4(T_f)={\lambda }_5(T_f)\\ ={\lambda }_6(T_f)={\lambda }_7(T_f)={\lambda }_8(T_f)=0.& (47)\end{array}$

Further in the interior of the set $\mathcal{U}$ where the controls are bounded in 0 ≤ u₁, θ₃, u₂ ≤ 1, we have the following optimality conditions being satisfied:

$\begin{array}{l}0=\frac{\partial \mathcal{H}}{\partial u_1}={\mathcal{B}}_1{u}_1+{\lambda }_1\beta (I_p+{\eta }_1{I}_a+{\eta }_2{I}_s)S\\ +{\lambda }_2(1-\sigma )\beta (I_p+{\eta }_1{I}_a+{\eta }_2{I}_s)V_p\\ +{\lambda }_3(1-\psi )\beta (I_p+{\eta }_1{I}_a+{\eta }_2{I}_s)V-{\lambda }_4\beta (I_p+{\eta }_1{I}_a+{\eta }_2{I}_s)S\\ -{\lambda }_4(1-\sigma )\beta (I_p+{\eta }_1{I}_a+{\eta }_2{I}_s)V_p\\ -{\lambda }_4(1-\psi )\beta (I_p+{\eta }_1{I}_a+{\eta }_2{I}_s)V,\\ 0=\frac{\partial \mathcal{H}}{\partial {\theta }_3}={\mathcal{B}}_2{\theta }_3-{\lambda }_1(1-h_3)S+{\lambda }_3(1-h_3)S,\\ 0=\frac{\partial \mathcal{H}}{\partial u_2}={\mathcal{B}}_3{u}_2-{\lambda }_{6}c{I}_a-{\lambda }_{7}c{I}_s+{\lambda }_{8}c{I}_a+{\lambda }_{8}c{I}_s,\end{array}$

such that

$\begin{array}{l}{u}_1^{*}=\frac{\begin{array}{l}[({\lambda }_4-{\lambda }_1)S+(1-\sigma )({\lambda }_4-{\lambda }_2)V_p\hfill \\ +(1-\psi )({\lambda }_4-{\lambda }_3)V]\beta (I_p+{\eta }_1{I}_a+{\eta }_2{I}_s)\hfill \end{array}}{{ℬ}_1},\\ {\theta }^{*}=\frac{({\lambda }_1-{\lambda }_3)(1-h_3)S}{{ℬ}_2},u_2^{*}=\frac{({\lambda }_6-{\lambda }_8)c{I}_a+({\lambda }_7-{\lambda }_8)c{I}_s}{{ℬ}_3}.\end{array}$

Given the time-dependent controls are bounded below by zero and above by one, we summarize the characterization as

${\mathcal{U}}^{*}=\{\begin{array}{c}0\text{ if }{k}_j^{*}\le 0\hfill \\ k_j^{*}\text{ if }0<k_j^{*}<1\hfill \\ 1\text{ if }{k}_j^{*}\ge 1\hfill \end{array}$

where j = 1, 2, 3; ${\mathcal{U}}^{*}=(u_1^{*},{\theta }_3^{*},u_2^{*})$ and

$\begin{array}{l}{k}_1^{*}=\frac{\begin{array}{c}\left[({\lambda }_4-{\lambda }_1)S+(1-\sigma )({\lambda }_4-{\lambda }_2)V_p+(1-\psi )({\lambda }_4-{\lambda }_3)V\right]\\ \beta (I_p+{\eta }_1{I}_a+{\eta }_2{I}_s)\end{array}}{{\mathcal{B}}_1},\end{array}$

$\begin{array}{ll}{k}_2^{*}=\frac{({\lambda }_1-{\lambda }_3)(1-h_3)S}{{\mathcal{B}}_2},& (48)\end{array}$

$\begin{array}{l}{k}_3^{*}=\frac{({\lambda }_6-{\lambda }_8)c{I}_a+({\lambda }_7-{\lambda }_8)c{I}_s}{{\mathcal{B}}_3}.\end{array}$

5 Numerical simulation and discussion

To validate theoretical findings of model system (1), we perform numerical simulations using MATLAB so as to gain more insight on equilibrium dynamics as well as on time profiles. Analytical analysis of the persistent COVID-19 equilibrium indicates that the model exhibits the phenomenon of backward bifurcation when certain conditions are satisfied as shown in Theorem 5, that is, if the double-dose and single-dose vaccines being administered are not 100% effective against suppressing COVID-19 persistent, then there exist a set of parameters that if they fulfill the condition that Γ₁ > Γ₂, and then bistability phenomenon arises. Hence, plotting the solution of cubic equation (20) (i.e., force of infection at equilibrium) as a function of the effective contact rate β results to Figure 2A which clearly depicts the phenomenon of backward bifurcation where COVID-19 persists even when the control reproduction number is below unity. The epidemiological implication of the occurrence of backward bifurcation signals that it will not be sufficient to decrease ${\mathcal{R}}_c$ below unity. Figure 2B which is obtained by varying vaccine efficacy among fully vaccinated individuals while all other parameters remain fixed as in Table 1 shows that vaccine with high efficacy leads to reduction of control reproduction numbers ${\mathcal{R}}_c$, ${\mathcal{R}}_p$, ${\mathcal{R}}_a$, and ${\mathcal{R}}_s.$ Figure 2B also corroborates with theoretical findings in equation (11) where it is shown that infectious pre-symptomatic individuals contribute significantly in the increase of control reproduction number in comparison with truly asymptomatic and symptomatic infectious individuals. In fact, ${\mathcal{R}}_p$ curve is higher than ${\mathcal{R}}_a$ and ${\mathcal{R}}_s$ curves for any given plausible set of parameter values.

FIGURE 2

Figure 2. (A) Illustration of the backward bifurcation phenomenon when vaccine is assumed to be imperfect. Parameters used are same as those in Table 1 except h₁ = 0.1, h₂ = 0.85, h₃ = 0.01, θ₁ = 0.015, θ₂ = 0.68, θ₃ = 0.075, σ = 0.5, v = 0.1, γ_a = 0.000165, γ_s = 0.00055, d = 0.001, ∧ = 0.0001which correspond to Γ₁ = 0.0056 > Γ₂ = 0.0049. Semi-logarithmic scale is used for a better view of bifurcation curves. The red dotted curve represents the unstable curve, while blue solid curve represents the stable curve. (B) Impact of vaccine efficacy ψ on control reproduction numbers. Semi-logarithmic scale is used for a clear view.

5.1 Single-dose vs. double-dose vaccination or both vaccination strategies

The COVID-19 model proposed adopted two vaccination strategies, that is, vaccinating susceptible individuals with either single-dose or double-dose vaccines. The parameters that relate to vaccination rates are θ₁, θ₂, and θ₃ where θ₁ and θ₂ represent first dose and second dose vaccination rates for a double-dose vaccine while θ₃ is the vaccination rate for a single-dose vaccines. While both strategies remain important in suppression of COVID-19 over a long time scale, it is worth investigating which strategy has more positive impact in mitigating COVID-19. Hence, fixing all other parameters constant as indicated in Table 1 while switching on and off vaccination rates, we obtain Figures 3A–D. Figure 3 reveals that single-dose vaccination strategy leads to a significant decline in COVID-19 prevalence in comparison with a double-dose vaccination strategy and baseline scenario (no vaccination). It is important to stress that these simulation results were generated when vaccine hesitancy is assumed to be equal among susceptible and partially vaccinated individuals (i.e., h₁ = h₂ = h₃ = 0.35). Hence, in a setting where vaccine reluctance (vaccine hesitancy) is more likely to occur, single-dose vaccination strategy would be more recommendable than a double dose-vaccine. This is due to the fact individuals who receive first dose of the double-dose vaccine may choose not to receive second dose as a result of vaccine hesitancy. The advantage of administering a single-dose vaccine is that vaccine hesitancy among partially vaccinated individuals is circumvented. Furthermore, it is important to note that both single-dose and double-dose vaccination strategies lead to a delay in COVID-19 peak when compared with the baseline scenario (no vaccine) indicated by a solid black curve in each figure. Interestingly, when comparing delay in COVID-19 peak between single-dose and double-dose vaccination strategies, there is no significant difference. In contrast, if both vaccination strategies are concurrently implemented, there is considerable difference in both COVID-19 prevalence and delaying in COVID-19 peak (when compared with baseline scenario and either single-dose or double-dose epidemic curve). That is, for the case of combined vaccination strategies, the epidemic curve is more flattened (see Figure 3). Thus, administering both single-dose and double-dose vaccines simultaneously has a positive impact in diminishing COVID-19 proliferation. Furthermore, a more flattened epidemic curve implies minimal burdening of healthcare infrastructure, given there will be no sudden surge of severely ill COVID-19 patients who need hospitalization.

FIGURE 3

Figure 3. Illustrate impact of single dose and double dose vaccination strategies on (A) the exposed individuals (B) pre-symptomatic infectious individuals (C) truly asymptomatic infectious individuals and (D) symptomatic infectious individuals. Parameters used remain as shown in Table 1 except, θ₁, θ₂, and θ₃ which are shown in the figure while initial conditions used are S(0) = 60000000, V_p(0) = 0, V(0) = 0, E(0) = 10954, I_p(0) = 10954, I_a(0) = 7322, I_s(0) = 7545, R(0) = 203968.

5.2 Effect of unwillingness to receive vaccines (vaccine hesitancy)

It is worth noting that despite COVID-19 vaccines being administered in almost every country in the world, the unwillingness of the general populace to accept COVID-19 jabs impend intervention measures put in place to end the COVID-19 pandemic. To elucidate on the impact of vaccine hesitancy on the COVID-19 transmission dynamics, we fix vaccination rates among susceptible and partially vaccinated individuals to a constant value (i.e., θ₁ = θ₂ = θ₃ = 0.1) and vary the fraction of individuals unwilling to receive COVID-19 jabs. That is, h₁, h₂, h₃ ∈ [0, 0.25, 0.50, 0.75, 0.90] while all other parameters remain as defined in Table 1. The first row figures of Figure 4 present the direct consequence of unwillingness to receive double-dose vaccine (by susceptible and partially vaccinated individuals) on the sizes of the infectious populations I_p, I_a, and I_s. The second row figures of Figure 4 show effect of unwillingness to receive single-dose vaccine (by the susceptible populations) on the sizes of infectious subgroups I_p, I_a, and I_s. The last row of Figure 4 shows the combined impact of vaccine hesitancy against single-dose and double-dose vaccines on the infectious cohorts I_p, I_a, and I_s. Figures 4A–C depict that an increase in vaccine hesitancy against single-dose vaccine has an adverse effect on COVID-19 prevalence. This is due to an increase in infectious population sizes (i.e., pre-symptomatic, truly asymptomatic and symptomatic) as vaccine hesitancy increases. An increase in infectious population sizes translates to an increase in likelihood of coming into contact with an infectious person which ultimately leads to more positive cases being confirmed in any given setting. Vaccine hesitancy against double-dose vaccine triggers an early appearance of COVID-19 peak among pre-symptomatic population than on both truly asymptomatic and symptomatic individuals. This observation is epidemiologically crucial, given in several epidemiological models the contribution of pre-symptomatic individuals in transmission of COVID-19 was ignored.

FIGURE 4

Figure 4. Illustrate the effect of vaccine hesitancy against double-dose vaccine, single-dose and both vaccines. (A) Show the effect of vaccine hesitancy against double-dose on pre-symptomatic infectious individuals. (B) Show the effect of vaccine hesitancy against double-dose on infectious truly asymptomatic individuals. (C) Illustrate the effect of vaccine hesitancy against double-dose on infectious symptomatic individuals. (D) Depict the effect of vaccine hesitancy against single-dose on pre-symptomatic infectious individuals. (E) Show effect of vaccine hesitancy against single-dose on infectious truly asymptomatic individuals. (F) Represents effect of vaccine hesitancy against single-dose on infectious symptomatic individuals. (G) Show the impact of vaccine hesitancy against both single and double dose on pre-symptomatic infectious individuals. (H) Show the impact of vaccine hesitancy against both single-dose and double-dose vaccines on infectious truly asymptomatic individuals. (I) Show the impact of vaccine hesitancy against both single dose and double dose vaccines on infectious symptomatic individuals. The blue solid epidemic curve represents the baseline scenario where there is no vaccine hesitancy.

Figures 4D–F present a scenario where there is vaccine hesitancy against single-dose vaccine (by susceptible individuals). It is observed that unwillingness to receive single-dose vaccine by susceptible individuals is more detrimental than vaccine hesitancy against double-dose vaccine. This is due to the pronounced early appearance of COVID-19 peaks among all the three infectious populations (I_p, I_a, I_s) and a wider increase in COVID-19 prevalence when compared to a case of vaccine hesitancy against double dose (see Figures 4A–C). Again, an early peak of COVID-19 epidemic curve is more notable among pre-symptomatic population than on both truly symptomatic and symptomatic populations. The combined impact of vaccine hesitancy against both single-dose and double-dose vaccines (by both susceptible and partially vaccinated cohorts) is presented in Figures 4G–I. It is clearly visible from Figures 4G–I that vaccine hesitancy against both single-dose and double-dose vaccines is disastrous as there is a rapid increase in COVID-19 prevalence by a wider margin when compared with the baseline curve (see the blue solid curve in each figure of Figure 4). Furthermore, a rise in vaccine hesitancy against the two vaccines being administered (single dose and double dose) triggers an early peak of COVID-19 epidemic curve which is again more pronounced when compared with vaccine hesitancy against each vaccination strategy (either single or double dose).

5.3 Simulations of the optimality system

To visualize the impacts of each of the three time-dependent controls and the effects of combining any two of the controls on the dynamics of COVID-19 transmission [81], we simulate the optimality system comprising the state system (38) and the adjoint system (45). Keeping in mind that the optimality system is a two-point boundary problem with both initial and transversality conditions (2) and (47), respectively, we therefore make use of fourth-order forward-backward Runge-Kutta method. We first solve the state system (38) forward in time with the initial condition (2) after initial guessed values for the controls. Then, we solve the adjoint system (45) backward in time with the transversality conditions (47) and update the controls using a convex combination of both the initial guessed values and current values of the control characterizations (48). We continue this iteration until the absolute error between the current and previous solutions of the optimality system is negligibly small [78, 82, 83]. The simulation is performed using MATLAB over the time interval [0, 120] days with the following initial conditions: S(0) = 60000000, V_p(0) = 0, V(0) = 0, E(0) = 10954, I_p(0) = 10954, I_a(0) = 7322, I_s(0) = 7545, R(0) = 203968. We use the same parameter values given in Table 1 together with the weight and rate constants ${\mathcal{A}}_m=1$, (m = 1, 2, 3), ${\mathcal{B}}_n=10$, (n = 1, 2, 3), and c = 1.

As shown in Figure 5, implementation of the non-pharmaceutical time-dependent control u₁(t) reduces the populations of pre-symptomatic, asymptomatic, and symptomatic individuals. The control profile shows that the non-pharmaceutical measure should be adhered to at maximum value for almost the implementation period to flatten the infectious curves. In Figure 6, we see that administration of the optimal vaccination control θ₃(t) for susceptible individuals with a single-dose vaccine helps in reducing the sizes of individuals in the subgroups I_p, I_a, and I_s when compared with the case without vaccination control.

FIGURE 5

Figure 5. Depict the impact of the non-pharmaceutical control u₁(t) on; (A) the infectious pre-symptomatic population. (B) The infectious truly asymptomatic population. (C) The infectious symptomatic population. (D) Represents the control profile of time-dependent control u₁(t).

FIGURE 6

Figure 6. Depicts the impact of the time-dependent optimal vaccination control θ₃(t) on: (A) the infectious pre-symptomatic population. (B) The infectious truly asymptomatic population. (C) The infectious symptomatic population. (D) Represents the control of time-dependent optimal vaccination control θ₃(t).

Since the time-dependent management control u₂(t) is only targeted at the populations of asymptomatic and symptomatic individuals, Figure 7 therefore displays the behaviors of the subgroups I_a and I_s in the presence and absence of u₂(t). It can be seen that the sizes of both asymptomatic and symptomatic individuals with control are lower when compared with the sizes without control. The control profile also indicates that management of infectious cases in the population should be sustained at maximum throughout the implementation period. Moreover, control profiles showing the combination of u₁(t) and θ₃(t), and the combination of u₁(t) and u₂(t) are displayed in Figure 8. In both profiles, we see that the non-pharmaceutical measure u₁(t) is maintained at maximum for almost 100% of the implementation period, while each of vaccination and management controls is maintained at maximum for 30% and 72% of the total implementation period, respectively. This implies that whenever the non-pharmaceutical measure is fully adhered to, less effort will be needed by either vaccination or management control to achieve optimal reduction of COVID-19 spread. Lastly, the influence of combination of both θ₃(t) and u₂(t) on the populations of individuals in the subgroups I_p(t), I_a(t), and I_s(t) is shown in Figure 9. We see that both vaccination and management controls are to be sustained at maximum throughout the required implementation period to reduce the sizes of infectious individuals in the population.

FIGURE 7

Figure 7. Illustrate the impact of the time-dependent screening-management control measure u₂(t) on: (A) the infectious truly asymptomatic population. (B) The infectious symptomatic population. (C) Represents the control profile of time-dependent screening-management control measure u₂(t).

FIGURE 8

Figure 8. Represents control profiles involving combinations of non-pharmaceutical measure with each of the optimal vaccination and management controls. (A) Represents control profiles for both optimal non-pharmaceutical control measure and optimal vaccination measure. (B) Represents control profiles for both optimal non-pharmaceutical control measure and screening-management control.

FIGURE 9

Figure 9. Demonstrate impact of combining the time-dependent optimal vaccination θ₃(t) and management u₂(t) controls on: (A) the infectious pre-symptomatic population (B) The infectious truly asymptomatic population (C) The infectious symptomatic population. (D) Represents the control profiles for both optimal vaccination and optimal screening-management of infectious individuals.

5.4 Cost-effectiveness assessment

Out of the single and several combinations of controls considered in this study, it is important to identify the main control strategy which minimizes COVID-19 transmission optimally at the lowest cost of implementation when available resources are limited. To do this, we make use of the two known methods, namely, average cost-effectiveness ratio (ACER) and incremental cost-effectiveness ratio (ICER) [84–86]. ACER is measured by dividing the cost of implementing a control strategy by the total benefits of such strategy in terms of cases averted. It is given by

$\begin{array}{ll}\text{ACER}=\frac{\text{Cost of control strategy}}{\text{Benefits of the control strategy}}& (49)\end{array}$

Unlike ACER, we use ICER to compare effectiveness of any two competing control strategies, and it is measured by finding the ratio of the difference between the costs of the two strategies to the difference between their health benefits. The ICER is calculated using

$\begin{array}{ll}\text{ICER}=\frac{\text{Difference in total costs}}{\text{Difference in total cases averted}}& (50)\end{array}$

As a result of the numerical simulations of the optimal control problem, we now present the values of the cases averted by each of the control strategies in an increasing order with their corresponding costs of implementation. We have used 49 and 50 to calculate ACER and ICER, respectively, as shown in Table 2.

TABLE 2

Table 2. Increasing order of COVID-19 cases averted by control strategies.

It can be observed that u₁(t) has the least ACER value in comparison with the other control strategies. This suggests that the non-pharmaceutical control is the most cost-effective strategy that can be implemented when resources are limited. To further confirm this result, we first carry out ICER analysis between θ₃(t) and u₂(t), and see that ICER(θ₃(t)) is greater than ICER(u₂(t)). This implies that implementing vaccination control is costlier and less-effective when compared with the optimal screening-management control. Hence, the control strategy θ₃(t) is excluded from the list of control strategies and we concentrate on the remaining as indicated in the Table 3.

TABLE 3

Table 3. ICER analysis without optimal vaccination control θ₃(t).

Comparison between the single control u₂(t) and the combination of θ₃(t) and u₂(t) in the Table 3 shows that ICER(θ₃(t), u₂(t)) is greater than ICER(u₂(t)), implying that implementing the combination of both vaccination and management controls is more expensive than the single implementation of optimal management control. Thus, we remove the combination θ₃(t) and u₂(t) from the list of the available control strategies and continue to analyze ICER for the remaining control strategies as presented in Table 4. We see that the non-pharmaceutical control u₁(t) strongly dominates the management control u₂(t) since ICER(u₁(t)) is less than ICER(u₂(t)). Therefore, we discard the control u₂(t) and re-compute ICER for the remaining strategies.

TABLE 4

Table 4. ICER analysis without the single θ₃(t) and the double (θ₃(t), u₂(t)).

Following similar procedure, we remove the combination of u₁(t) and θ₃(t) since it is strongly dominated by the single control strategy u₁(t) as seen in the Table 5. Hence, we are left with two control strategies as indicated in Table 6. Consequently, the combination of u₁(t) and u₂(t) is removed from the analysis since ICER(u₁(t), u₂(t)) is greater than ICER(u₁(t)). Hence, single implementation of the non-pharmaceutical control strategy is the most cost-effective of all the strategies. This result is in total agreement with the ACER result obtained earlier.

TABLE 5

Table 5. ICER analysis without single θ₃(t), double (θ₃(t), u₂(t)), and single u₂(t).

TABLE 6

Table 6. ICER between single u₁(t) and double (u₁(t), u₂(t)).

6 Conclusion

This study considered a COVID-19 epidemiological framework that modeled transition from one healthy status to another using a deterministic non-linear ordinary differential equations. Given vaccines are widely accepted by policymakers as an important intervention measure for curbing COVID-19 proliferation, we incorporated into the proposed model non-mandatory vaccines. Contrary to previous work [33, 35, 41–44, 46] on COVID-19, we considered a setting where susceptible individuals have a choice of either being vaccinated with a single- or double-dose vaccine. Furthermore, a fraction of individuals who initially choose double-dose vaccine refuse the second dose and therefore remain partially vaccinated (these individuals can be infected with COVID-19 upon exposure with infectious individuals). This differs with [33, 46] studies that assumed vaccinated individuals acquired sufficient vaccine-induced immunity, and therefore, they could not be infected.

The assumption that vaccinated individuals could not be infected (see [33, 46]) hindered the important epidemiological phenomenon of backward bifurcation which has been observed in this study. For instance, findings from equilibrium analysis of the model suggest that if both single-dose and double-dose vaccines administered to the general populace are perfect and permanent then intervention measures that aim to decrease control reproduction number below one will be sufficient in curbing COVID-19 (see Theorem 5). However, if vaccines administered are imperfect and non-permanent, there exists a parameter space where the phenomenon of backward bifurcation may arise leading to persistence of COVID-19 even when the control reproduction number is below one. This observation hints that policymakers (in particular governments across the world) and medical practitioners should advocate production of vaccines that render optimal efficacy (as well as vaccines that provide a livelong vaccine-induced immunity) so that breakthrough infections are prevented at all cost. The studies done by Buonomo et al. [33], Oduro et al. [35], Buckner et al. [41], Choi and Shim [42], Mukandavire et al. [43], Deng et al. [44], and Peter et al. [46] did not exhibit any bistability phenomenon. Both analytical and numerical results on the three control reproduction numbers indicate that the control reproduction number associated with the number of infectious individuals generated as a result of contact with pre-symptomatic individuals is always larger than the other two control reproduction numbers associated with symptomatic and truly asymptomatic infectious individuals, respectively. This finding is important as far as COVID-19 transmission dynamics is being demystified both at epidemiological and biological level. The revelation regarding extent of contribution of pre-symptomatic individuals can inform policymakers and medical practitioners on the importance of implementing early screening and quarantine of individuals exposed to COVID-19 to avert spiraling of COVID-19 cases. The studies conducted in Tiwari et al. [39], Rai et al. [40], Kumar et al. [48], Majumder et al. [87], srivastav et al. [88], Oduro et al. [35], Buckner et al. [41], Choi and Shim [42], Mukandavire et al. [43], Deng et al. [44], and Peter et al. [46] never captured the epidemiological implication of pre-symptomatic individuals due to the fact that they did not incorporate pre-symptomatic cohort. Also during the time of their studies, there was no universal consensus on whether pre-symptomatic individuals contributed in the spread of COVID-19 as it is now known [8].

The time profile findings reveal a profound difference in reduction of COVID-19 prevalence when either single-dose vaccine or double-dose vaccine is administered among susceptible population. Administering single-dose vaccines leads to a significant decline of COVID-19 prevalence than double-dose vaccines, especially in a setting where vaccine hesitancy is present. This epidemiologically implies that in a setting where vaccine hesitancy is rampant, a single-dose vaccine should be recommended. This is due to the fact that single-dose vaccine circumvents the second dose vaccine hesitancy likely to be witnessed among partially vaccinated individuals who may refuse second COVID-19 jab due to reasons such as mild (severe) side effects experienced with the first dose. This is despite the existing evidence from medical practitioners that the vaccine benefits considerably outweigh the risks of side effects associated with COVID-19 vaccines being administered globally. Nevertheless, in a scenario where both vaccines are simultaneously administered, there is a significant decline in COVID-19 prevalence, regardless of the level of vaccine hesitancy. It is important to stress that the findings regarding epidemiological impact of combined effect of single-dose and double-dose vaccines have not been documented in any literature (known so far) possibly due to the fact that no existing study has analyzed a model that incorporated both vaccines (single-dose and double-dose vaccines) and vaccine hesitancy simultaneously. Even the recent studies [33, 35, 41–44, 46] on COVID-19 considered only a single-dose or a double-dose vaccines separately.

Investigation on the impact of unwillingness to receive either single-dose or double-dose vaccine reveals that vaccine hesitancy triggers an early peak of COVID-19 outbreak. Sudden surge of COVID-19 prevalence implies the available medical facilities will be overburdened by severely ill COVID-19 patients who are in need of urgent medical attention. In addition, if vaccine hesitancy against both vaccines (that is single dose and double dose) is rampant in the general populace, then COVID-19 epidemic curve peaks faster than vaccine hesitancy against either a single-dose or double-dose vaccine. This implies that countermeasures put in place to control COVID-19 proliferation need to go beyond providing free COVID-19 jabs but also demystify the benefits associated with COVID-19 vaccines so as to raise the level of vaccine acceptance in the general population.

The extension of the COVID-19 model to an optimal control problem by incorporating three time-dependent controls also shed some epidemiological insights on how optimal non-pharmaceutical and pharmaceutical intervention measures may influence COVID-19 transmission dynamics. The three time-dependent controls which are broadly categorized as either non-pharmaceutical (NPIS) and pharmaceutical included preventive control measures (such as social/physical distancing, wearing face mask, hand washing, and adherent to education on proper hygiene) among susceptible and vaccinated cohorts, vaccination of susceptible individuals with a single dose-vaccine, and screening-management of infectious symptomatic individuals and truly asymptomatic individuals. Previous studies [33, 35, 41–44, 46] on COVID-19 did not consider these three time-dependent control measures as we have considered in this study; hence, findings obtained in the study are unique and different from their findings. Optimal simulation results show that if susceptible and both vaccinated cohorts (either partially vaccinated or fully vaccinated) adhere to COVID-19 non-pharmaceutical intervention measures, then COVID-19 pandemic can be easily eradicated. Thus, in the midst of a COVID-19 pandemic advising the public to adhere to non-pharmaceutical preventive measures even after vaccination can have positive impact through reduction of COVID-19 prevalence to almost zero. Furthermore, cost-effectiveness analysis has shown that non-pharmaceutical intervention control is the most effective and cheapest when compared with pharmaceutical controls (management-screening controls of symptomatic and truly asymptomatic and vaccination control). Optimal management of infectious symptomatic and truly asymptomatic individuals significantly reduces COVID-19 prevalence but does not lead to a delayed COVID-19 peak. On the one hand, optimal vaccination of susceptible individuals has 2-fold benefits, that is, reduction of COVID-19 prevalence as well as delaying the peak of COVID-19 outbreak. This allows policymakers and medical practitioners to have sufficient time to plan and manage the existing COVID-19 cases.

7 Limitations

It is imperative to stress that there are few limitations to this study. First, data on COVID-19 vaccine hesitancy are limited, and the research that analyzed individual level of vaccine hesitancy using data obtained in US was conducted before widespread roll-out of COVID-19 vaccines [89]. Hence, the data may be unreliable in a setting where COVID-19 mitigation measures incorporate vaccines. On the one hand, social-economic status and demographic factors were not considered as were beyond the scope of this study. However, consideration of such factors may impact COVID-19 vaccine hesitancy. Moreover, individual perception toward each type of vaccine available varies which could influence level of vaccine hesitancy. This aspect is not captured in this study as we only investigated the general impact of vaccine hesitancy toward either single- or double-dose vaccines without singling a specific type of vaccine.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

IW: Conceptualization, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – original draft, Writing – review & editing. SO: Conceptualization, Formal analysis, Investigation, Methodology, Project administration, Software, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing. RL: Conceptualization, Funding acquisition, Investigation, Project administration, Resources, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing. KO: Conceptualization, Investigation, Methodology, Project administration, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing.

Funding

The author(s) declare that no financial support was received for the research, authorship, and/or publication of this article.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fams.2023.1292443/full#supplementary-material

References

1. World Health Organization. WHO COVID-19 Preparedness and Response Progress Report-1 February to 30 June 2020. Geneva: WHO. (2020).

Google Scholar

2. WHO. Coronavirus Disease (COVID-19): Situation Report. Geneva: WHO. (2020) p. 198.

Google Scholar

3. Chang X, Wang J, Liu M, Jin Z, Han D. Study on an SIHRS model of COVID-19 pandemic with impulse and time delay under media coverage. IEEE Access. (2021) 9:49387–97. doi: 10.1109/ACCESS.2021.3064632

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Akkilic AN, Sabir Z, Raja MAZ, Bulut H. Numerical treatment on the new fractional-order SIDARTHE COVID-19 pandemic differential model via neural networks. Eur Phys J Plus. (2022) 137:334. doi: 10.1140/epjp/s13360-022-02525-w

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Chatterjee AN, Basir FA, Ahmad B, Alsaedi A, A. fractional-order compartmental model of vaccination for COVID-19 with the fear factor. Mathematics. (2022) 10:1451. doi: 10.3390/math10091451

CrossRef Full Text | Google Scholar

6. Anggriani N, Beay LK. Modeling of COVID-19 spread with self-isolation at home and hospitalized classes. Results Phys. (2022) 36:105378. doi: 10.1016/j.rinp.2022.105378

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Chen TM, Rui J, Wang QP, Zhao ZY, Cui JA, Yin L, et al. mathematical model for simulating the phase-based transmissibility of a novel coronavirus. Infect Dis Poverty. (2020) 9:1–8. doi: 10.1186/s40249-020-00640-3

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Alleman TW, Vergeynst J, De Visscher L, Rollier M, Torfs E, Nopens I, et al. Assessing the effects of non-pharmaceutical interventions on SARS-CoV-2 transmission in Belgium by means of an extended SEIQRD model and public mobility data. Epidemics. (2021) 37:100505. doi: 10.1016/j.epidem.2021.100505

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Tyrrell DAJ, Bynoe ML. Cultivation of viruses from a high proportion of patients with colds. Lancet. (1966) 1:76–7. doi: 10.1016/S0140-6736(66)92364-6

PubMed Abstract | CrossRef Full Text | Google Scholar

10. WHO Coronavirus (COVID-19) Dashboard. (2023). Available online at: https://covid19.who.int/ (accessed November 7, 2023).

Google Scholar

11. Bayati M, Noroozi R, Ghanbari-Jahromi M, Jalali FS. Inequality in the distribution of COVID-19 vaccine: a systematic review. Int J Equity Health. (2022) 21:1–9. doi: 10.1186/s12939-022-01729-x

CrossRef Full Text | Google Scholar

12. Baker CM, Chades I, McVernon J, Robinson AP, Bondell H. Optimal allocation of PCR tests to minimise disease transmission through contact tracing and quarantine. Epidemics. (2021) 37:100503. doi: 10.1016/j.epidem.2021.100503

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Yamey G. Schäferhoff M, Hatchett R, Pate M, Zhao F, McDade KK. Ensuring global access to COVID-19 vaccines. Lancet. (2020) 395:1405–6. doi: 10.1016/S0140-6736(20)30763-7

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Zhang SX, Wang Y, Rauch A, Wei F. Unprecedented disruption of lives and work: Health, distress and life satisfaction of working adults in China one month into the COVID-19 outbreak. Psychiatry Res. (2020) 288:112958. doi: 10.1016/j.psychres.2020.112958

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Riaz MMA, Ahmad U, Mohan A, dos Santos Costa AC, Khan H, Babar MS, et al. Global impact of vaccine nationalism during COVID-19 pandemic. Trop Med Health. (2021) 49:1–4. doi: 10.1186/s41182-021-00394-0

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Le TT, Andreadakis Z, Kumar A, Román RG, Tollefsen S, Saville M, et al. The COVID-19 vaccine development landscape. Nat Rev Drug Discov. (2020) 19:305–6. doi: 10.1038/d41573-020-00073-5

CrossRef Full Text | Google Scholar

17. ÇakmaklıC, Demiralp S, Kalemli-”Ozcan e, Yeşiltaş S, Yıldırım MA. The Economic Case for Global Vaccinations: An Epidemiological Model With International Production Networks. Cambridge: National Bureau of Economic Research. (2021).

Google Scholar

18. Stephenson J. Unequal Access to COVID-19 Vaccines Leaves Less-Wealthy Countries More Vulnerable, Poses Threat to Global Immunity. Philadelphia, PA: American Medical Association. (2021).

PubMed Abstract | Google Scholar

19. Tupy IS, Silva FF, Diniz GFC, Montenegro RL, de Queiroz Stein A, Ferraz D. Resilient Regions in Brazil: Unfolding the Effects of COVID-19 From a Socioeconomic Perspective. Thousand Oaks, CA: International Regional Science Review (2022). p. 01600176221145878.

Google Scholar

20. Bollyky TJ, Gostin LO, Hamburg MA. The equitable distribution of COVID-19 therapeutics and vaccines. JAMA. (2020) 323:2462–3. doi: 10.1001/jama.2020.6641

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Nkengasong JN, Ndembi N, Tshangela A, Raji T. COVID-19 vaccines: how to ensure Africa has access. Nature. (2020) 586:197–9. doi: 10.1038/d41586-020-02774-8

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Gostin LO, Karim SA, Mason Meier B. Facilitating access to a COVID-19 vaccine through global health law. J Law, Med Ethics. (2020) 48:622–6. doi: 10.1177/1073110520958892

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Bassi LL. Allocating COVID-19 vaccines globally: An urgent need. In: JAMA Health Forum. Philadelphia, PA: American Medical Association. (2021) p. 1–2.

PubMed Abstract | Google Scholar

24. World Health Organization. Vaccine Equity Its Only Impossible Until Its Done (2022). Available online at: https://www.who.int/campaigns/vaccine-equity (accessed November 24, 2023).

Google Scholar

25. Dooling K. COVID-19 Vaccine Prioritization: Work Group Considerations. (2020).

Google Scholar

26. Joint Committee on Vaccination and Immunisation: Interim advice on priority groups for COVID-19 vaccination (2020). Available online at: https://www.gov.uk/coronavirus (accessed November 7, 2023).

Google Scholar

27. Chen JT, Krieger N. Revealing the unequal burden of COVID-19 by income, race, ethnicity, and household crowding: US county versus zip code analyses. J Public Health Manage Pract. (2021) 27:S43–56. doi: 10.1097/PHH.0000000000001263

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Ogedegbe G, Ravenell J, Adhikari S, Butler M, Cook T, Francois F, et al. Assessment of racial ethnic disparities in hospitalization and mortality in patients with COVID-19 in New York City. JAMA Network Open. (2020) 3:e2026881–e2026881. doi: 10.1001/jamanetworkopen.2020.26881

CrossRef Full Text | Google Scholar

29. Vahidy FS, Nicolas JC, Meeks JR, Khan O, Pan A, Jones SL, et al. Racial and ethnic disparities in SARS-CoV-2 pandemic: analysis of a COVID-19 observational registry for a diverse US metropolitan population. BMJ Open. (2020) 10:e039849. doi: 10.1136/bmjopen-2020-039849

PubMed Abstract | CrossRef Full Text | Google Scholar

30. McClung N, Chamberland M, Kinlaw K, Matthew DB, Wallace M, Bell BP, et al. The Advisory Committee on Immunization Practices' ethical principles for allocating initial supplies of COVID-19 vaccine—United States, 2020. Morbid Mortal Wkly Rep. (2020) 69:1782. doi: 10.15585/mmwr.mm6947e3

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Kukreti S, Rifai A, Padmalatha S, Lin CY Yu T, Ko WC, et al. Willingness to obtain COVID-19 vaccination in general population: a systematic review and meta-analysis. J Glob Health. (2022) 12:50. doi: 10.7189/jogh.12.05006

CrossRef Full Text | Google Scholar

32. Peretti-Watel P, Seror V, Cortaredona S, Launay O, Raude J, Verger P, et al. A future vaccination campaign against COVID-19 at risk of vaccine hesitancy and politicisation. Lancet Infect Dis. (2020) 20:769–70. doi: 10.1016/S1473-3099(20)30426-6

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Buonomo B, Della Marca R. d'Onofrio A, Groppi M. A behavioural modelling approach to assess the impact of COVID-19 vaccine hesitancy J Theo Biol. (2022) 534:110973. doi: 10.1016/j.jtbi.2021.110973

CrossRef Full Text | Google Scholar

34. Neumann-Böhme S, Varghese NE, Sabat I, Barros PP, Brouwer W, van Exel J, et al. Once we have it, will we use it? A European survey on willingness to be vaccinated against COVID-19. Eur J Health Econ. (2020) 21:977–82. doi: 10.1007/s10198-020-01208-6

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Oduro B, Miloua A, Apenteng OO, Osei PP. COVID-19 vaccination hesitancy model: The impact of vaccine education on controlling the outbreak in the United States. medRxiv [Preprint]. (2021) 2021–05. doi: 10.1101/2021.05.21.21257612

CrossRef Full Text | Google Scholar

36. Kirzinger A, Sparks G, Hamel L, Lopes L, Kearney A, Stokes M, et al. KFF COVID-19 vaccine monitor: July 2021. Kaiser Family Foundation. (2021).

Google Scholar

37. Ngonghala CN, Iboi E, Eikenberry S, Scotch M, MacIntyre CR, Bonds MH, et al. Mathematical assessment of the impact of non-pharmaceutical interventions on curtailing the 2019 novel coronavirus. Math Biosci. (2020) 325:108364. doi: 10.1016/j.mbs.2020.108364

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Wangari IM, Sewe S, Kimathi G, Wainaina M, Kitetu V, Kaluki W. Mathematical modelling of COVID-19 transmission in Kenya: A model with reinfection transmission mechanism. Comput Math Methods Med. (2021) 2021:1–18. doi: 10.1155/2021/5384481

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Tiwari PK, Rai RK, Khajanchi S, Gupta RK, Misra AK. Dynamics of coronavirus pandemic: Effects of community awareness and global information campaigns. Eur Phys J Plus. (2021) 136:994. doi: 10.1140/epjp/s13360-021-01997-6

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Rai RK, Khajanchi S, Tiwari PK, Venturino E, Misra AK. Impact of social media advertisements on the transmission dynamics of COVID-19 pandemic in India. J Appl Math Comp. (2022) 1–26. doi: 10.1007/s12190-021-01507-y

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Buckner JH, Chowell G, Springborn MR. Dynamic prioritization of COVID-19 vaccines when social distancing is limited for essential workers. Proc Nat Acad Sci. (2021) 118:e2025786118. doi: 10.1073/pnas.2025786118

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Choi W, Shim E. Optimal strategies for vaccination and social distancing in a game-theoretic epidemiologic model. J Theor Biol. (2020) 505:110422. doi: 10.1016/j.jtbi.2020.110422

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Mukandavire Z, Nyabadza F, Malunguza NJ, Cuadros DF, Shiri T, Musuka G. Quantifying early COVID-19 outbreak transmission in South Africa and exploring vaccine efficacy scenarios. PLoS ONE. (2020) 15:e0236003. doi: 10.1371/journal.pone.0236003

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Deng J, Tang S, Shu H. Joint impacts of media, vaccination and treatment on an epidemic Filippov model with application to COVID-19. J Theo Biol. (2021) 523:110698. doi: 10.1016/j.jtbi.2021.110698

PubMed Abstract | CrossRef Full Text | Google Scholar

45. d'Onofrio A, Manfredi P. Information-related changes in contact patterns may trigger oscillations in the endemic prevalence of infectious diseases. Journal of Theoretical Biology. 2009;256(3):473–478. doi: 10.1016/j.jtbi.2008.10.005

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Peter OJ, Panigoro HS, Abidemi A, Ojo MM, Oguntolu FA. Mathematical model of COVID-19 pandemic with double dose vaccination. Acta Biotheor. (2023) 71:9. doi: 10.1007/s10441-023-09460-y

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Johnson J. Johnson & Johnson COVID-19 Vaccine Authorized by US FDA For Emergency Use-First Single-Shot Vaccine in Fight Against Global Pandemic (2021). Available online at: https://www.janssen.com/COVID19 (accessed November 8, 2023).

Google Scholar

48. Kumar Rai R, Kumar Tiwari P, Khajanchi S. Modeling the influence of vaccination coverage on the dynamics of COVID-19 pandemic with the effect of environmental contamination. Math Methods Appl Sci. (2023) 46:12425–53. doi: 10.1002/mma.9185

CrossRef Full Text | Google Scholar

49. Cavanaugh AM, Spicer KB, Thoroughman D, Glick C, Winter K. Reduced risk of reinfection with SARS-CoV-2 after COVID-19 vaccination-Kentucky, May-June 2021. Morbid Mortal Wkly Rep. (2021) 70:1081. doi: 10.15585/mmwr.mm7032e1

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Britton A, Slifka KMJ, Edens C, Nanduri SA, Bart SM, Shang N, et al. Effectiveness of the Pfizer-BioNTech COVID-19 vaccine among residents of two skilled nursing facilities experiencing COVID-19 outbreaks-Connecticut, December 2020-February 2021. Morbid Mortal Wkly Rep. (2021) 70:396. doi: 10.15585/mmwr.mm7011e3

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Bender JK, Brandl M, Höhle M, Buchholz U, Zeitlmann N. Analysis of asymptomatic and presymptomatic transmission in SARS-CoV-2 outbreak, Germany, 2020. Emerg Infect Dis. (2021) 27:1159. doi: 10.3201/eid2704.204576

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Taboe HB, Asare-Baah M, Iboi EA, Ngonghala CN. Critical assessment of the impact of vaccine-type and immunity on the burden of COVID-19. Math Biosci. (2023) 360:108981. doi: 10.1016/j.mbs.2023.108981

CrossRef Full Text | Google Scholar

53. Lipsitch M, Krammer F, Regev-Yochay G, Lustig Y, Balicer RD. SARS-CoV-2 breakthrough infections in vaccinated individuals: Measurement, causes and impact. Nat Rev Immunol. (2022) 22:57–65. doi: 10.1038/s41577-021-00662-4

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Curley B. How Long Does Immunity From COVID-19 Vaccination Last? Healthline (2021).

Google Scholar

55. Seow J, Graham C, Merrick B, Acors S, Pickering S, Steel KJA, et al. Longitudinal observation and decline of neutralizing antibody responses in the three months following SARS-CoV-2 infection in humans. Nat Microbiol. (2020) 5:1598–607. doi: 10.1038/s41564-020-00813-8

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Dan JM, Mateus J, Kato Y, Hastie KM, Yu ED, Faliti CE, et al. Immunological memory to SARS-CoV-2 assessed for up to 8 months after infection. Science. (2021) 371:eabf4063. doi: 10.1126/science.abf4063

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Gatto M, Bertuzzo E, Mari L, Miccoli S, Carraro L, Casagrandi R, et al. Spread and dynamics of the COVID-19 epidemic in Italy: Effects of emergency containment measures. Proc Nat Acad Sci. (2020) 117:10484–91. doi: 10.1073/pnas.2004978117

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Zhou C. Evaluating new evidence in the early dynamics of the novel coronavirus COVID-19 outbreak in Wuhan, China with real time domestic traffic and potential asymptomatic transmissions. medRxiv [Preprint]. (2020) 2020-02. doi: 10.1101/2020.02.15.20023440

CrossRef Full Text | Google Scholar

59. Linton NM, Kobayashi T, Yang Y, Hayashi K, Akhmetzhanov AR, Jung Sm, et al. Incubation period and other epidemiological characteristics of 2019 novel coronavirus infections with right truncation: A statistical analysis of publicly available case data. J Clin Med. (2020) 9:538. doi: 10.3390/jcm9020538

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Sun H, Qiu Y, Yan H, Huang Y, Zhu Y, Chen SX. Tracking and predicting COVID-19 epidemic in China mainland. medRxiv [Preprint]. (2020) 2020-02. doi: 10.1101/2020.02.17.20024257

CrossRef Full Text | Google Scholar

61. Zou L, Ruan F, Huang M, Liang L, Huang H, Hong Z, et al. SARS-CoV-2 viral load in upper respiratory specimens of infected patients. N Engl J Med. (2020) 382:1177–9. doi: 10.1056/NEJMc2001737

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Anderson RM, Heesterbeek H, Klinkenberg D, Hollingsworth TD. How will country-based mitigation measures influence the course of the COVID-19 epidemic? Lancet. (2020) 395:931–4. doi: 10.1016/S0140-6736(20)30567-5

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Alene M, Yismaw L, Assemie MA, Ketema DB, Mengist B, Kassie B, et al. Magnitude of asymptomatic COVID-19 cases throughout the course of infection: A systematic review and meta-analysis. PLoS ONE. (2021) 16:e0249090. doi: 10.1371/journal.pone.0249090

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Chagla Z. The BNT162b2 (BioNTech/Pfizer) vaccine had 95% efficacy against COVID-19 7 days after the 2nd dose. Ann Intern Med. (2021) 174:JC15. doi: 10.7326/ACPJ202102160-015

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Van den Driessche P, Watmough J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Math Biosci. (2002) 180:29–48. doi: 10.1016/S0025-5564(02)00108-6

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Castillo-Chavez C, Blower S, van den Driessche P, Kirschner D, Yakubu AA. Mathematical Approaches for Emerging and Reemerging Infectious Diseases: Models, Methods, and Theory. vol 126. Berlin: Springer Science & Business Media (2002).

Google Scholar

67. Anderson B, Jackson J, Sitharam M. Descartes' rule of signs revisited. Am Math Monthly. (1998) 105:447–51. doi: 10.1080/00029890.1998.12004907

CrossRef Full Text | Google Scholar

68. La Salle JP. The stability of dynamical systems. SIAM. (1976). doi: 10.1137/1.9781611970432

CrossRef Full Text | Google Scholar

69. Thongtha A, Modnak C. Optimal COVID-19 epidemic strategy with vaccination control and infection prevention measures in Thailand. Infect Dis Model. (2022) 7:835–55. doi: 10.1016/j.idm.2022.11.002

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Castillo-Chavez C, Song B. Dynamical models of tuberculosis and their applications. Math Biosci Eng. (2004) 1:361–404. doi: 10.3934/mbe.2004.1.361

PubMed Abstract | CrossRef Full Text | Google Scholar

71. Wangari IM, Stone L. Backward bifurcation and hysteresis in models of recurrent tuberculosis. PLoS One. (2018) 13:e0194256. doi: 10.1371/journal.pone.0194256

PubMed Abstract | CrossRef Full Text | Google Scholar

72. Wangari IM, Davis S, Stone L. Backward bifurcation in epidemic models: Problems arising with aggregated bifurcation parameters. Appl Math Model. (2016) 40:1669–75. doi: 10.1016/j.apm.2015.07.022

CrossRef Full Text | Google Scholar

73. Wangari IM. Backward Bifurcation and Reinfection in Mathematical Models of Tuberculosis. Melbourne: RMIT University (2017).

Google Scholar

74. Okyere E, Olaniyi S, Bonyah E. Analysis of Zika virus dynamics with sexual transmission route using multiple optimal controls. Sci Afr. (2020) 9:e00532. doi: 10.1016/j.sciaf.2020.e00532

CrossRef Full Text | Google Scholar

75. Ghosh M, Olaniyi S, Obabiyi OS. Mathematical analysis of reinfection and relapse in malaria dynamics. Appl Math Comput. (2020) 373:125044. doi: 10.1016/j.amc.2020.125044

CrossRef Full Text | Google Scholar

76. Bonyah E, Khan MA, Okosun KO. Gómez-Aguilar J. On the co-infection of dengue fever and Zika virus Optim Cont Appl Method. (2019) 40:394–421. doi: 10.1002/oca.2483

CrossRef Full Text | Google Scholar

77. Purwati UD, Riyudha F, Tasman H, et al. Optimal control of a discrete age-structured model for tuberculosis transmission. Heliyon. (2020) 6. doi: 10.1016/j.heliyon.2019.e03030

PubMed Abstract | CrossRef Full Text | Google Scholar

78. Lenhart S, Workman JT. Optimal Control Applied to Biological Models. Boca Raton, FL: CRC Press (2007). doi: 10.1201/9781420011418

CrossRef Full Text | Google Scholar

79. Olaniyi S, Okosun K, Adesanya S, Lebelo R. Modelling malaria dynamics with partial immunity and protected travellers: Optimal control and cost-effectiveness analysis. J Biol Dyn. (2020) 14:90–115. doi: 10.1080/17513758.2020.1722265

PubMed Abstract | CrossRef Full Text | Google Scholar

80. Pontryagin LS. Mathematical Theory of Optimal Processes. Wiley (1962).

Google Scholar

81. Olaniyi S, Akanni J, Adepoju O. Optimal control and cost-effectiveness analysis of an illicit drug use population dynamics. J Appl Nonlinear Dynam. (2023) 12:133–46. doi: 10.5890/JAND.2023.03.010

CrossRef Full Text | Google Scholar

82. Olaniyi S, Ajala OA, Abimbade SF. Optimal control analysis of a mathematical model for recurrent malaria dynamics. Operat Res Forum (Springer). (2023) 4:14. doi: 10.1007/s43069-023-00197-5

CrossRef Full Text | Google Scholar

83. Olaniyi S, Obabiyi O, Okosun K, Oladipo A, Adewale S. Mathematical modelling and optimal cost-effective control of COVID-19 transmission dynamics. Eur Phys J Plus). (2020) 135:938. doi: 10.1140/epjp/s13360-020-00954-z

PubMed Abstract | CrossRef Full Text | Google Scholar

84. Asamoah JKK, Okyere E, Abidemi A, Moore SE, Sun GQ, Jin Z, et al. Optimal control and comprehensive cost-effectiveness analysis for COVID-19. Results Phys. (2022) 33:105177. doi: 10.1016/j.rinp.2022.105177

PubMed Abstract | CrossRef Full Text | Google Scholar

85. Olaniyi S, MukamuriM, Okosun K, Adepoju O. Mathematical analysis of a social hierarchy-structured model for malaria transmission dynamics. Results Phys. (2022) 34:104991. doi: 10.1016/j.rinp.2021.104991

PubMed Abstract | CrossRef Full Text | Google Scholar

86. Olaniyi S, Abimbade SF, Ajala OA, Chuma FM. Efficiency and economic analysis of intervention strategies for recurrent malaria transmission. Qual Quant. (2023) 2023:1–19. doi: 10.1007/s11135-023-01664-1

CrossRef Full Text | Google Scholar

87. Majumder M, Tiwari PK, Pal S. Impact of saturated treatments on HIV-TB dual epidemic as a consequence of COVID-19: Optimal control with awareness and treatment. Nonlinear Dyn. (2022) 109:143–76. doi: 10.1007/s11071-022-07395-6

PubMed Abstract | CrossRef Full Text | Google Scholar

88. Srivastav AK, Tiwari PK, Srivastava PK, Ghosh M, Kang Y, A. mathematical model for the impacts of face mask, hospitalization and quarantine on the dynamics of COVID-19 in India: Deterministic versus stochastic. Math Biosci Eng. (2021) 18:182–213. doi: 10.3934/mbe.2021010

PubMed Abstract | CrossRef Full Text | Google Scholar

89. Gerretsen P, Kim J, Caravaggio F, Quilty L, Sanches M, Wells S, et al. Individual determinants of COVID-19 vaccine hesitancy. PLoS ONE. (2021) 16:e0258462. doi: 10.1371/journal.pone.0258462

PubMed Abstract | CrossRef Full Text | Google Scholar