Cost–benefit analysis of the COVID-19 vaccination model incorporating different infectivity reductions

The spread and control of coronavirus disease 2019 (COVID-19) present a worldwide economic and medical burden to public health. It is imperative to probe the effect of vaccination and infectivity reductions in minimizing the impact of COVID-19. Therefore, we analyze a mathematical model incorporating different infectivity reductions. This work provides the most economical and effective control methods for reducing the impact of COVID-19. Using data from Ghana as a sample size, we study the sensitivity of the parameters to estimate the contributions of the transmission routes to the effective reproduction number $R_e.$ We also devise optimal interventions with cost–benefit analysis that aim to maximize outcomes while minimizing COVID-19 incidences by deploying cost-effectiveness and optimization techniques. The outcomes of this work contribute to a better understanding of COVID-19 epidemiology and provide insights into implementing interventions needed to minimize the COVID-19 burden in similar settings worldwide.

1 Introduction

Coronavirus disease 2019 (COVID-19) is one of the infectious diseases that has caused a global medical and financial burden. Due to its high transmission rate, the World Health Organization (WHO) officially announced the prevalence of it as a global pandemic on 8th March 2020. The coronavirus spreads directly or indirectly from one infectious living host cell to another by replicating and causing infectious disease in the host (humans and/or animals) [1]. Therefore, it is imperative to investigate how the virus spreads in the community and devise strategies to halt the disease’s transmission.

Globally, COVID-19 has caused over 6 million deaths, and by March 2022, approximately 480 million incidences had been reported, according to [2]. However, the COVID-19 reinfection scenario using the SIR model shows that transmission dynamics could arise as a result of immunity waning, even in cases where the force of reinfection is relatively weak [3]. As mentioned in [4], the COVID-19 model, which entails infection through objects contaminated with SARS-CoV-2, is suggested to be made public knowledge.

COVID-19 spreads quickly, threatening global health and igniting a pandemic. This pandemic has had a major effect on other sectors, especially the socioeconomic sector [5]. As of 30 April 2022, there were 512,466,045 reported incidences and 6,257,512 COVID-19 fatalities worldwide, according to the Worldometer. With 6,046,467 COVID-19 cases reported overall, Indonesia ranks seventh in Asia and 18th globally [6].

These findings highlight the need to investigate the spread of COVID-19 in order to minimize its transmission. Using a deterministic model is one way that mathematics is crucial to simulating the epidemic phenomenon of the disease’s spread. The analysis of the COVID-19 model with declining immunity has advanced significantly. For the latter, the natural immunity period is defined by the vaccine efficacy level, which dictates when to start the mass vaccination strategy based on models involving symptomatic and asymptomatic infected populations, as explained in [6].

Mathematical models have been an important tool in epidemiology since the 18th century and can be used to determine the level of spread of infectious disease, aside from medical and biological research studies and strategies for controlling diseases. A modified susceptible–exposed–infectious–removed (SEIR) model with vaccination, quarantine, and isolated SVEQIMR is used to analyze the spread of COVID-19 disease in [7]. In this article, compartmental models were used to demonstrate the pattern of vaccination and control strategies for infectious diseases.

Furthermore, the compartmental model is used to study the behavior of COVID-19 after the introduction of vaccines in [8]. In this article, the population is divided into distinct compartments to demonstrate the control efficacy of infectious diseases in India. Again, a parameterized nonlinear SEIHR model to analyze the transmission of coronavirus disease in Indonesia using the compartmental model is analyzed in [9].

The transmission dynamics, basic reproduction number, and control measures of COVID-19 are analyzed using the compartmental model in [10]. In this article, effective public health interventions were proposed to control the spread of the disease in Wuhan, China. The cost-effective analysis, global stability, and control strategy for the spread of COVID-19 in Ghana are analyzed using the compartmental model in [11]. In this article, control measures were outlined by formulating the SEAIRV model to investigate human–environment–human transmission. However, these studies do not capture the effects of reducing the infectivity rate on the disease’s transmission.

The aforementioned issues have driven this study to explore the effects of vaccination and infectivity reductions on the spread of COVID-19, devise optimal control interventions, and analyze the cost–benefits of implementing the interventions. The subsequent sections are as follows: in Section 2, we present the design of the epidemiological model together with the definition of the parameters and variables in the model. In Section 3, we analyze the positivity of the model’s solutions, computation of the models’ reproduction numbers, and stability of the model. In Section 4, we present the model parameter estimation and sensitivity analysis. In Section 5, we present the numerical simulations and model analysis. In Section 6, we present the optimal control strategies and analyze the cost–benefit analysis. We finally present the concluding remarks of this study in Section 7.

2 $\mathit{S}\mathit{V}\mathit{E}\mathit{Q}\mathit{I}\mathit{M}\mathit{R}$ model formulation

The flowchart below illustrates the mechanism of the model. We denote $N\left(t\right)$ as the total population divided into seven different compartments: susceptible individuals at a given time $S\left(t\right)$, vaccinated individuals at a given time $V\left(t\right)$, exposed individuals at a given time $E\left(t\right)$, individuals under quarantine at a given time $Q\left(t\right)$, infected individuals at a given time $I\left(t\right),$ isolated individuals at a given time $M\left(t\right)$, and individuals recovered at a given time $R\left(t\right)$. We assume that the vaccinated individuals become susceptible again due to the vaccine’s inefficacy. Again, a portion of susceptible individuals undergo self-quarantine, while others enter the exposed class. Some of the exposed individuals become infected and then recover naturally without any special treatment. The remaining proportion of individuals either go through quarantine or become infected. It is also assumed that the individuals in the susceptible and vaccinated compartments would come into direct or indirect contact with the individuals in the exposed, quarantined, and isolated compartments. Therefore, we introduce infectivity reductions $r_e,r_q$, and $r_m$ attributed to the exposed, quarantined, and isolated compartments, respectively. The introduction of these different infectivity reductions is to minimize transmission from the potential carriers of the virus. Again, we assume that the quarantined and isolated individuals are not part of the active population at the time of incidence, so the total active population for the standard incidence at any given time is $N_1=N\left(t\right)-Q\left(t\right)-M\left(t\right).$ Therefore, the forces for infection of susceptible and vaccinated classes are given as ${\rho }_1=\frac{{\beta }_s\left(I+r_{e}E+r_{q}Q+r_{m}M\right)}{N_1}$ and ${\rho }_2=\frac{{\beta }_v\left(I+r_{e}E+r_{q}Q+r_{m}M\right)}{N_1}$, respectively. Recovered individuals become susceptible again once they come into contact with the possible carriers of the COVID-19 virus. The description of the model’s parameters is given in Table 1. The optimal control interventions ${\kappa }_p\left(t\right),{\kappa }_v\left(t\right)$, and ${\kappa }_i\left(t\right)$ are explained in the subsequent section.

Table 1

Table 1. Definition of model parameters.

The total population $N\left(t\right)$ in Eq. (1) is defined based on Figure 1 as follows:

$N\left(t\right)=S\left(t\right)+V\left(t\right)+E\left(t\right)+Q\left(t\right)+I\left(t\right)+M\left(t\right)+R\left(t\right).(1)$

Figure 1

Figure 1. SVEQIMR model flowchart.

The following nonlinear ordinary differential equations illustrate the model.

$\left\{\begin{array}{l}\frac{\mathit{d}\mathit{S}}{\mathit{d}\mathit{t}}=\mathit{\Lambda }+\mathit{\sigma }\mathit{V}+{\mathit{\delta }}_1\mathit{Q}+\mathit{\tau }\mathit{R}-{\mathit{\rho }}_1\mathit{S}-{\mathit{\omega }}_1\mathit{S},\\ \frac{\mathit{d}\mathit{V}}{\mathit{d}\mathit{t}}=\mathit{\vartheta }\mathit{S}-{\mathit{\rho }}_2\mathit{V}-{\mathit{\omega }}_2\mathit{V},\\ \frac{\mathit{d}\mathit{E}}{\mathit{d}\mathit{t}}={\mathit{\rho }}_1\mathit{S}+{\mathit{\rho }}_2\mathit{V}-{\mathit{\omega }}_3\mathit{E},\\ \frac{\mathit{d}\mathit{Q}}{\mathit{d}\mathit{t}}=\mathit{ŋ}\mathit{S}+\mathit{\phi }\mathit{E}-{\mathit{\omega }}_4\mathit{Q},\\ \frac{\mathit{d}\mathit{I}}{\mathit{d}\mathit{t}}=\mathit{\gamma }\mathit{E}-{\mathit{\omega }}_5\mathit{I},\\ \frac{\mathit{d}\mathit{M}}{\mathit{d}\mathit{t}}={\mathit{\delta }}_2\mathit{Q}+\mathit{\alpha }\mathit{I}-{\mathit{\omega }}_6\mathit{M},\\ \frac{\mathit{d}\mathit{R}}{\mathit{d}\mathit{t}}={\mathit{\epsilon }}_{\mathit{e}}\mathit{E}+{\mathit{\epsilon }}_{\mathit{i}}\mathit{I}+{\mathit{\epsilon }}_{\mathit{m}}\mathit{M}-{\mathit{\omega }}_7\mathit{R}.\end{array}\right.(2)$

where $\begin{array}{l}{\rho }_1=\frac{{\beta }_s\left(I+r_{e}E+r_{q}Q+r_{m}M\right)}{N_1},{\rho }_2=\frac{{\beta }_v\left(I+r_{e}E+r_{q}Q+r_{m}M\right)}{N_1},{\omega }_1=\mathrm{\mu }+\vartheta +ŋ,{\omega }_2=\mathrm{\mu }+\mathrm{\sigma },\\ {\omega }_3=\mathrm{\mu }+\phi +\gamma +{\epsilon }_e,{\omega }_4=\mathrm{\mu }+{\delta }_1+{\delta }_2,{\omega }_5=\mathrm{\mu }+\alpha +{\epsilon }_i\end{array}$, ${\omega }_6=\mathrm{\mu }+{\epsilon }_m$, and ${\omega }_7=\mathrm{\mu }+\tau .$

With initial conditions $S\left(0\right)\ge 0,V\left(0\right)\ge 0,E\left(0\right)\ge 0,Q\left(0\right)\ge 0,I\left(0\right)\ge 0,M\left(0\right)\ge 0,\text{and\hspace{0.17em}}R\left(0\right)\ge 0,$ all initial parameters of model (2) are nonnegative.

3 $\mathit{S}\mathit{V}\mathit{E}\mathit{Q}\mathit{I}\mathit{M}\mathit{R}$ model mathematical analysis

The detailed analysis of the model is conducted mathematically to show the positivity of solutions and their boundedness within a specific domain.

3.1 Positivity and boundedness of solutions

The variables and parameters of model (2) are nonnegative since it is based on a population with COVID-19.

Theorem 1. Define $W\left(t\right)=S\left(t\right),V\left(t\right),E\left(t\right),Q\left(t\right),I\left(t\right),M\left(t\right),R\left(t\right)$; if $W\left(0\right)\ge 0,$ then $W\left(t\right)\ge 0,$ and its solutions and initial values are nonnegative for $t>0$ and bounded in the region ${\mathrm{R}}_{+}^7.$

Proof: Let us consider the following instance where there exists an initial time $t_i$ such that

$\min \left\{W\left(t_i\right)\right\}>0$ and $\min \left\{W\left(t\right)\right\}>0$ for all $t\in \left[0\right.,\left.t_i\right).$

Here, $W\left(t\right)=S\left(t\right),V\left(t\right),E\left(t\right),Q\left(t\right),I\left(t\right),M\left(t\right),R\left(t\right).$ Without the loss of generalization, $\min \left\{W\left(t_i\right)\right\}=S\left(t_i\right)$.

Therefore, $S\left(t_i\right)=0,V\left(t_i\right)>0,Q\left(t_i\right)>0,R\left(t_i\right)>0$, and $S\left(t\right)>0,$ for all $t\in \left[0\right.,\left.t_i\right).$ However,

$\frac{dS\left(t_i\right)}{dt}=\mathrm{\Lambda }+\sigma V\left(t_i\right)+{\delta }_{1}Q\left(t_i\right)+\tau R\left(t_i\right)>0,\mathrm{\Lambda }\ge 0,S\left(t_i\right)>S\left(0\right)\ge 0.$

This contradicts the claim $S\left(t_i\right)=0.$ Therefore, $S\left(t\right)>0$ for all $t\ge 0.$ This shows that all the solutions are positive for $t\ge 0$ in all other cases.

Theorem 2. To prove the boundedness of model (2), define a positive invariant set as $K=\left\{\left(S_0,V_0,E_0,Q_0,I_0,M_0,R_0\right)\in R_{+}^7:N\left(t\right)\le \frac{\mathrm{\Lambda }}{\mathrm{\mu }}\right\}$ and attract positive solutions.

Proof: Considering Eq. 1.

The rate of change in the total population is given as

$N^{\prime }\left(t\right)=S^{\prime }\left(t\right)+V^{\prime }\left(t\right)+Q^{\prime }\left(t\right)+I^{\prime }\left(t\right)+M^{\prime }\left(t\right)+R^{\prime }\left(t\right),$

$N^{\prime }\left(t\right)=\mathrm{\Lambda }-\mu \left(S\left(t\right)+V\left(t\right)+E\left(t\right)+Q\left(t\right)+I\left(t\right)+M\left(t\right)+R\left(t\right)\right),$

$N^{\prime }=\mathrm{\Lambda }–\mu N.(3)$

From Eq. (3), it follows that,

$N\left(\mathrm{t}\right)\le \frac{\mathrm{\Lambda }}{\mu }+N\left(0\right)e^{-\mu t}.$

Then, $0<N\left(\mathrm{t}\right)\le \frac{\mathrm{\Lambda }}{\mu }.$ Thus, $N\left(\mathrm{t}\right)$ is bounded, and all solutions in the $K$ approach enter or remain in $K.$ If $t\to \infty ,0\le N\left(\mathrm{t}\right)$ shows that $N\left(t\right)$ is a set of positive invariant and is in the region $R_{+}^7.$ This theorem proves the existence of COVID-19 at a given time in an area that was not infected with COVID-19 disease, and all the initial state variables are positive. This completes the proof.

3.2 Disease-free equilibrium and effective reproduction number

The disease equilibrium point, ${\xi }_0$, of the model is achieved by equating system (2) to 0. The result is defined as follows:

$\left\{\begin{array}{l}\mathit{\Lambda }+\mathit{\sigma }\mathit{V}+{\mathit{\delta }}_1\mathit{Q}+\mathit{\tau }\mathit{R}-{\mathit{\rho }}_1\mathit{S}-{\mathit{\omega }}_1\mathit{S}=0,\\ \mathit{\vartheta }\mathit{S}-{\mathit{\rho }}_2\mathit{V}-{\mathit{\omega }}_2\mathit{V}=0,\\ {\mathit{\rho }}_1\mathit{S}+{\mathit{\rho }}_2\mathit{V}-{\mathit{\omega }}_3\mathit{E}=0,\\ \mathit{ŋ}\mathit{S}+\mathit{\phi }\mathit{E}-{\mathit{\omega }}_4\mathit{Q}=0,\\ \mathit{\gamma }\mathit{E}-{\mathit{\omega }}_5\mathit{I}=0,\\ {\mathit{\delta }}_2\mathit{Q}+\mathit{\alpha }\mathit{I}-{\mathit{\omega }}_6\mathit{M}=0,\\ {\mathit{\epsilon }}_{\mathit{e}}\mathit{E}+{\mathit{\epsilon }}_{\mathit{i}}\mathit{I}+{\mathit{\epsilon }}_{\mathit{m}}\mathit{M}-{\mathit{\omega }}_7\mathit{R}=0.\end{array}\right.(4)$

At disease-free equilibrium (DFE), we set $E=Q=I=M=R=0$, and the following results are obtained from Eq. (4):

$\mathrm{\Lambda }+\sigma V-{\omega }_{1}S=0,$

$\vartheta S={\omega }_{2}V,$

$S_0=\frac{\mathrm{\Lambda }\left(\mathrm{\mu }+\mathrm{\sigma }\right)}{\left(\mathrm{\mu }+\vartheta +ŋ\right)\left(\mathrm{\mu }+\mathrm{\sigma }\right)-\mathrm{\sigma }\vartheta },$

$V_0=\frac{\mathrm{\Lambda }\vartheta }{\left(\mathrm{\mu }+\vartheta +ŋ\right)\left(\mathrm{\mu }+\mathrm{\sigma }\right)-\mathrm{\sigma }\vartheta }.$

At disease-free equilibrium,

${\xi }_0=\left(\frac{\mathrm{\Lambda }\left(\mathrm{\mu }+\mathrm{\sigma }\right)}{\left(\mathrm{\mu }+\vartheta +ŋ\right)\left(\mathrm{\mu }+\mathrm{\sigma }\right)-\mathrm{\sigma }\vartheta },\frac{\mathrm{\Lambda }\vartheta }{\left(\mathrm{\mu }+\vartheta +ŋ\right)\left(\mathrm{\mu }+\mathrm{\sigma }\right)-\mathrm{\sigma }\vartheta },0,0,0,0,0\right).(5)$

The corresponding Jacobian matrix of the system (2) evaluated at ${\xi }_0=\left(S_0,V_0,0,0,0,0,0\right)$ (see, Eq. (5)) to obtain the disease-free equilibrium Jacobian matrix $J_{{\xi }_0}$ is given as follows:

$J_{{\xi }_0}=\left(\begin{array}{ccccccc}-{\omega }_1& \mathrm{\sigma }& -\frac{{\beta }_s{S}_0}{N_0}{r}_e& {\delta }_1-\frac{{\beta }_s{S}_0}{N_0}{r}_q& -\frac{{\beta }_s{S}_0}{N_0}& -\frac{{\beta }_s{S}_0}{N_0}{r}_m& 0\\ \vartheta & -{\omega }_2& -\frac{{\beta }_v{V}_0}{N_0}{r}_e& -\frac{{\beta }_v{V}_0}{N_0}{r}_q& -\frac{{\beta }_v{V}_0}{N_0}& -\frac{{\beta }_v{V}_0}{N_0}{r}_m& 0\\ 0& 0& \frac{{\beta }_s{S}_0+{\beta }_v{V}_0}{N_0}{r}_e-{\omega }_3& \frac{{\beta }_s{S}_0+{\beta }_v{V}_0}{N_0}{r}_q& \frac{{\beta }_s{S}_0+{\beta }_v{V}_0}{N_0}& \frac{{\beta }_s{S}_0+{\beta }_v{V}_0}{N_0}{r}_m& 0\\ ŋ& 0& \phi & -{\omega }_4& 0& 0& 0\\ 0& 0& \gamma & 0& -{\omega }_5& 0& 0\\ 0& 0& 0& {\delta }_2& \alpha & -{\omega }_6& 0\\ 0& 0& {\epsilon }_e& 0& {\epsilon }_i& {\epsilon }_m& -\mathrm{\mu }\end{array}\right),(6)$

where $N_0=N-Q-M$ using ${\xi }_0$; hence, $N_0=S_0+V_0.$

Using $J_{{\xi }_0}$, Eq. (6) the following matrices can be deduced to evaluate the effective reproduction number $R_e$. We define the next-generation matrix $G=F{V}^{-1}$ as the square matrix, which consists of matrix $F$ representing new infections and matrix V.

$F=\left(\begin{array}{cccc}\frac{{\beta }_s{S}_0+{\beta }_v{V}_0}{N_0}{r}_e& \frac{{\beta }_s{S}_0+{\beta }_v{V}_0}{N_0}& \frac{{\beta }_s{S}_0+{\beta }_v{V}_0}{N_0}{r}_q& \frac{{\beta }_s{S}_0+{\beta }_v{V}_0}{N_0}{r}_m\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right),V=\left(\begin{array}{cccc}{\omega }_3& 0& 0& 0\\ -\phi & {\omega }_4& 0& 0\\ -\gamma & 0& {\omega }_5& 0\\ 0& -{\delta }_2& -\alpha & {\omega }_6\end{array}\right),$

$V^{-1}=\left(\begin{array}{cccc}\frac{1}{{\omega }_3}& 0& 0& 0\\ \frac{\phi }{{\omega }_3{\omega }_4}& \frac{1}{{\omega }_4}& 0& 0\\ \frac{\gamma }{{\omega }_3{\omega }_5}& 0& \frac{1}{{\omega }_5}& 0\\ \frac{\alpha \gamma {\omega }_4+\phi {\delta }_2{\omega }_5}{{\omega }_3{\omega }_4{\omega }_5{\omega }_6}& \frac{{\delta }_2}{{\omega }_4{\omega }_6}& \frac{\alpha }{{\omega }_5{\omega }_6}& \frac{1}{{\omega }_6}\end{array}\right),G=F{V}^{-1}=\left(\begin{array}{cccc}{g}_{11}& g_{12}& g_{13}& g_{14}\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right),(7)$

where $\begin{array}{l}{g}_{11}=\frac{{\beta }_s{S}_0+{\beta }_v{V}_0}{N_0{\omega }_3}{r}_e+\frac{\left({\beta }_s{S}_0+{\beta }_v{V}_0\right)\gamma }{N_0{\omega }_3{\omega }_5}+\frac{\left({\beta }_s{S}_0+{\beta }_v{V}_0\right)\phi }{N_0{\omega }_3{\omega }_4}{r}_q+\frac{\left({\beta }_s{S}_0+{\beta }_v{V}_0\right)\left({\omega }_4\alpha \gamma +{\omega }_5\phi {\delta }_2\right)}{N_0{\omega }_3{\omega }_4{\omega }_5{\omega }_6}{r}_m,g_{12}=\frac{{\beta }_s{S}_0+{\beta }_v{V}_0}{{N_0\omega }_4}{r}_q+\\ \frac{\left({\beta }_s{S}_0+{\beta }_v{V}_0\right){\delta }_2}{N_0{\omega }_4{\omega }_6}{r}_m,g_{13}=\frac{{\beta }_s{S}_0+{\beta }_v{V}_0}{N_0{\omega }_5}+\frac{\left({\beta }_s{S}_0+{\beta }_v{V}_0\right)\alpha }{N_0{\omega }_5{\omega }_6}{r}_m,g_{14}=\frac{{\beta }_s{S}_0+{\beta }_v{V}_0}{N_0{\omega }_6}{r}_m.\end{array}$

The effective reproduction number, $R_e$, of the model is evaluated as the spectral radius of matrix $G$; that is, $\mathcal{P}\left(F{V}^{-1}\right)$ given as

$R_e=R_{eE}+R_{eI}+R_{eQ}+R_{eM},$

where $R_{eE}=\frac{{\beta }_s{S}_0+{\beta }_v{V}_0}{{N_0\omega }_3}{r}_e,R_{eI}=\frac{\left({\beta }_s{S}_0+{\beta }_v{V}_0\right)\gamma }{{N_0\omega }_3{\omega }_5},R_{eQ}=\frac{\left({\beta }_s{S}_0+{\beta }_v{V}_0\right)\phi }{N_0{\omega }_3{\omega }_4}{r}_q,$ and $R_{eM}=\frac{\left({\beta }_s{S}_0+{\beta }_v{V}_0\right)\left({\omega }_4\alpha \gamma +{\omega }_5\phi {\delta }_2\right)}{{N_0\omega }_3{\omega }_4{\omega }_5{\omega }_6}{r}_m$.

Substituting $S_0,V_0$, and $N_0$ into Eq. (7), gives

$R_e=\left(\frac{{\beta }_s\left(\mathrm{\mu }+\mathrm{\sigma }\right)+{\beta }_v\vartheta }{\left(\mathrm{\mu }+\mathrm{\sigma }+\vartheta \right)}\right)\left(\frac{r_e}{{\omega }_3}+\frac{\gamma }{{\omega }_3{\omega }_5}+\frac{\phi r_q}{{\omega }_3{\omega }_4}+\frac{\left({\omega }_4\alpha \gamma +{\omega }_5\phi {\delta }_2\right)r_m}{{\omega }_3{\omega }_4{\omega }_5{\omega }_6}\right).(8)$

Therefore, there is a unique equilibrium in the model, which implies that $R_e$ in Eq. (8) is unique. The effective reproduction number, $R_e,$ indicates the number of secondary infections that one infected person can produce if they come into contact with people living in a safe zone. $R_e$ can be summarized as follows: the exposed class contributed to a significant secondary infection, which is the first term of $R_e,$ where a proportion of $\frac{\gamma }{{\omega }_3}$ individuals entered the infected class. The infected class generated a significant secondary infection, which is the second term of $R_e.$ The quarantined class contributed to a significant secondary infection, which is the third term of $R_e.$ The isolated class contributed to a significant secondary infection, which is the last term of $R_e.$ We consider the dynamics of vaccination and infectivity reductions on $R_e$ as follows in Eqs 9–11 respectively.

• In the absence of vaccination, the reproduction number, $R_a^v,$ is given as

$R_a^v={\beta }_s\left(\frac{r_e}{{\omega }_3}+\frac{\gamma }{{\omega }_3{\omega }_5}+\frac{\phi r_q}{{\omega }_3{\omega }_4}+\frac{\left({\omega }_4\alpha \gamma +{\omega }_5\phi {\delta }_2\right)r_m}{{\omega }_3{\omega }_4{\omega }_5{\omega }_6}\right).(9)$

• In the absence of infectivity reductions, the reproduction number, $R_a^i,$ is given as

$R_a^i=\left(\frac{{\beta }_s\left(\mathrm{\mu }+\mathrm{\sigma }\right)+{\beta }_v\vartheta }{\left(\mathrm{\mu }+\mathrm{\sigma }+\vartheta \right)}\right)\left(\frac{1}{{\omega }_3}+\frac{\gamma }{{\omega }_3{\omega }_5}+\frac{\phi }{{\omega }_3{\omega }_4}+\frac{\left({\omega }_4\alpha \gamma +{\omega }_5\phi {\delta }_2\right)}{{\omega }_3{\omega }_4{\omega }_5{\omega }_6}\right).(10)$

• In the absence of vaccination and infectivity reductions, the reproduction number, $R_a^n,$ is given as

$R_a^n={\beta }_s\left(\frac{1}{{\omega }_3}+\frac{\gamma }{{\omega }_3{\omega }_5}+\frac{\phi }{{\omega }_3{\omega }_4}+\frac{\left({\omega }_4\alpha \gamma +{\omega }_5\phi {\delta }_2\right)}{{\omega }_3{\omega }_4{\omega }_5{\omega }_6}\right).(11)$

3.3 Stability analysis

Next, since the system has unique equilibrium points, we check its stability.

3.4 Local stability of disease-free equilibrium ${\mathit{\xi }}_0$

The Jacobian matrix, $J_{{\xi }_0}$, is given in Equation 7. It can be seen from (6) that ${\lambda }_1=-{\omega }_1,{\lambda }_2=-{\omega }_2,{\text{and}\lambda }_3=-{\omega }_7$ are the three eigenvalues of $J_{{\xi }_0}$. The other eigenvalues are derived from the reduced matrix in Eq. (12):

$\left(J_{1{\xi }_0}-\lambda I_4\right)=\left(\begin{array}{cccc}A-{\omega }_3-\lambda & A_1& A_2& A_3\\ \phi & -{\omega }_4-\lambda & 0& 0\\ \gamma & 0& -{\omega }_5-\lambda & 0\\ 0& {\delta }_2& \alpha & -{\omega }_6-\lambda \end{array}\right),(12)$

where $A=\frac{{\beta }_s{S}_0+{\beta }_v{V}_0}{N_0}{r}_e,A_1=\frac{{\beta }_s{S}_0+{\beta }_v{V}_0}{N_0}{r}_q,A_2=\frac{{\beta }_s{S}_0+{\beta }_v{V}_0}{N_0},$ and $A_3=\frac{{\beta }_s{S}_0+{\beta }_v{V}_0}{N_0}{r}_m.$

Therefore, the remaining eigenvalues are the roots of the following characteristic polynomial:

${\mathrm{\Gamma }}_4\left(\lambda \right)={\lambda }^4+H_3{\lambda }^3+H_2{\lambda }^2+H_1\lambda +H_0=0,(13)$

where

$H_3={\omega }_3+{\omega }_4+{\omega }_5+{\omega }_6-A,$

$H_2={\omega }_3\left({\omega }_4+{\omega }_5+{\omega }_6\right)+{\omega }_4\left({\omega }_5+{\omega }_6\right)+{\omega }_5{\omega }_6-A\left({\omega }_4+{\omega }_5+{\omega }_6\right)-A_2\gamma ,$

$H_1={\omega }_3\left({\omega }_4{\omega }_5+{\omega }_4{\omega }_6+{\omega }_5{\omega }_6\right)+{\omega }_4{\omega }_5{\omega }_6-A\left({\omega }_4{\omega }_5+{\omega }_4{\omega }_6+{\omega }_5{\omega }_6\right)-B{\delta }_{rt}\left(a_3+a_4\right)-A_1\phi \left({\omega }_5+{\omega }_6\right)-A_2\gamma \left({\omega }_4+{\omega }_6\right)-A_3\left(\alpha \gamma +\phi {\delta }_2\right),$

$H_0={\omega }_3{\omega }_4{\omega }_5{\omega }_6-A{\omega }_4{\omega }_5{\omega }_6-A_1\phi {\omega }_5{\omega }_6-A_2\gamma {\omega }_4{\omega }_6-A_3\left(\alpha \gamma {\omega }_4+\phi {\delta }_2{\omega }_5\right),$

$H_0=1-\left[\frac{A}{{\omega }_3}+\frac{A_2\gamma }{{\omega }_3{\omega }_5}+\frac{A_1\phi }{{\omega }_4{\omega }_5}+\frac{A_3\left(\alpha \gamma {\omega }_4+\phi {\delta }_2{\omega }_5\right)}{{\omega }_3{\omega }_4{\omega }_5{\omega }_6}\right],$

$H_0=1-\left(R_{eE}+R_{eI}+R_{eQ}+R_{eM}\right),$

$H_0=1-R_e.$

Therefore, using the Routh–Hurwitz stability conditions, the roots of ${\mathrm{\Gamma }}_4\left(\lambda \right)$ in Eq. (13) of the reduced matrix have negative real parts if the following conditions hold: $H_3>0,H_2>0,H_1>0,$ and $H_0>0.$ It is obvious that the condition holds if $R_e<1$. This proves that in model (2), disease-free equilibrium, ${\xi }_0$, is locally stable if $R_e<1$ and unstable if $R_e>1$.

3.5 Existence of endemic equilibrium of the model

Let us consider $P^{*}=\left(S^{*},V^{*},E^{*},Q^{*},I^{*},M^{*},R^{*}\right)$ as the endemic equilibrium for system (2) and equate the derivative to zero (0). The following results are obtained:

$\left\{\begin{array}{l}\mathit{\Lambda }+\mathit{\sigma }{\mathit{V}}^{*}+{\mathit{\delta }}_1{\mathit{Q}}^{*}+\mathit{\tau }{\mathit{R}}^{*}-\mathit{Z}{\mathit{\beta }}_{\mathit{s}}{\mathit{S}}^{*}-{\mathit{\omega }}_1{\mathit{S}}^{*}=0,\\ \mathit{\vartheta }{\mathit{S}}^{*}-\mathit{Z}{\mathit{\beta }}_{\mathit{v}}{\mathit{V}}^{*}-{\mathit{\omega }}_2{\mathit{V}}^{*}=0,\\ \mathit{Z}\left({\mathit{\beta }}_{\mathit{s}}{\mathit{S}}^{*}+{{\mathit{\beta }}_{\mathit{v}}\mathit{V}}^{*}\right)-{\mathit{\omega }}_3{\mathit{E}}^{*}=0,\\ \mathit{ŋ}{\mathit{S}}^{*}+\mathit{\phi }{\mathit{E}}^{*}-{\mathit{\omega }}_4{\mathit{Q}}^{*}=0,\\ \mathit{\gamma }{\mathit{E}}^{*}-{\mathit{\omega }}_5{\mathit{I}}^{*}=0,\\ {\mathit{\delta }}_2{\mathit{Q}}^{*}+\mathit{\alpha }{\mathit{I}}^{*}-{\mathit{\omega }}_6{\mathit{M}}^{*}=0,\\ {\mathit{\epsilon }}_{\mathit{e}}{\mathit{E}}^{*}+{\mathit{\epsilon }}_{\mathit{i}}{\mathit{I}}^{*}+{\mathit{\epsilon }}_{\mathit{m}}{\mathit{M}}^{*}-{\mathit{\omega }}_7{\mathit{R}}^{*}=0,\end{array}\right.(14)$

where $Z=\frac{\left(I^{*}+r_e{E}^{*}+r_q{Q}^{*}+r_m{M}^{*}\right)}{N_1^{*}},$ for $N_1^{*}=N^{*}-Q^{*}-M^{*}.$

Then, $P^{*}=\left(S^{*},V^{*},E^{*},Q^{*},I^{*},M^{*},R^{*}\right),$ and solving Eq. 14 simultaneously gives the following

$\left\{\begin{array}{l}{\mathit{S}}^{*}=\frac{\mathit{\vartheta }{\mathit{\omega }}_4{\mathit{\omega }}_6\left(\mathit{Z}{\mathit{\beta }}_{\mathit{v}}+{\mathit{\omega }}_2\right)\left(\mathit{\Lambda }+{\mathit{y}}_1{\mathit{I}}^{*}\right)}{{\mathit{y}}_4},\\ {\mathit{V}}^{*}=\frac{\mathit{\vartheta }{\mathit{\omega }}_4{\mathit{\omega }}_6\left(\mathit{\Lambda }+{\mathit{y}}_1{\mathit{I}}^{*}\right)}{{\mathit{y}}_4},\\ {\mathit{E}}^{*}=\left(\frac{{\mathit{\omega }}_5}{\mathit{\gamma }}\right){\mathit{I}}^{*},\\ {\mathit{Q}}^{*}=\frac{1}{{\mathit{\omega }}_4}\left(\frac{\mathit{\gamma }\mathit{ŋ}\mathit{\vartheta }{\mathit{\delta }}_2{\mathit{\omega }}_4{\mathit{\omega }}_6\left(\mathit{Z}{\mathit{\beta }}_{\mathit{v}}+{\mathit{\omega }}_2\right)\left(\mathit{\Lambda }+{\mathit{y}}_1{\mathit{I}}^{*}\right)+\mathit{\phi }{\mathit{\omega }}_5{\mathit{y}}_3{\mathit{I}}^{*}}{\mathit{\gamma }{\mathit{y}}_3}\right),\\ {\mathit{M}}^{*}=\frac{1}{{\mathit{\omega }}_6}\left(\frac{\mathit{\gamma }\mathit{ŋ}\mathit{\vartheta }{\mathit{\omega }}_4{\mathit{\omega }}_6\left(\mathit{Z}{\mathit{\beta }}_{\mathit{v}}+{\mathit{\omega }}_2\right)\left(\mathit{\Lambda }+{\mathit{y}}_1{\mathit{I}}^{*}\right)+\left(\mathit{\phi }{\mathit{\omega }}_5+\mathit{\alpha }\mathit{\gamma }\right){\mathit{y}}_3{\mathit{I}}^{*}}{\mathit{\gamma }{\mathit{y}}_4}\right),\\ {\mathit{R}}^{*}=\frac{1}{{\mathit{\omega }}_7}\left(\frac{\left({{\mathit{\epsilon }}_{\mathit{e}}\mathit{\omega }}_5+\mathit{\gamma }{\mathit{\epsilon }}_{\mathit{i}}\right){\mathit{I}}^{*}}{\mathit{\gamma }}+\frac{{\mathit{\epsilon }}_{\mathit{m}}}{{\mathit{\omega }}_6}\left(\frac{\mathit{\gamma }\mathit{ŋ}\mathit{\vartheta }{\mathit{\delta }}_2{\mathit{\omega }}_4{\mathit{\omega }}_6\left(\mathit{Z}{\mathit{\beta }}_{\mathit{v}}+{\mathit{\omega }}_2\right)\left(\mathit{\Lambda }+{\mathit{y}}_1{\mathit{I}}^{*}\right)+\left(\mathit{\phi }{\mathit{\omega }}_5+\mathit{\alpha }\mathit{\gamma }\right){\mathit{y}}_3{\mathit{I}}^{*}}{\mathit{\gamma }{\mathit{y}}_4}\right)\right).\end{array}\right.(15)$

Substituting the expressions for $S^{*},V^{*},\text{and\hspace{0.17em}}{E}^{*}$ in Eq. 15 into the fifth equation of Eq. 14 and simplifying, we obtain the following equation for $I^{*}$:

$\left(\frac{{\omega }_4{\omega }_6\left(\mathrm{\Lambda }+y_1{I}^{*}\right)\left(Z{\beta }_v+{\omega }_2+\vartheta \right)}{y_4}\right)-\left(\frac{{{{\omega }_3\omega }_{5}I}^{*}}{\gamma }\right)=0.$

Simplifying the above equation gives Eq. (16)

$I^{*}=\frac{\mathrm{\Lambda }\gamma {\omega }_4{\omega }_6\left(Z{\beta }_v+{\omega }_2+\vartheta \right)}{{{\omega }_3\omega }_4{y}_4-\gamma {\omega }_4{\omega }_6{y}_1\left(Z{\beta }_v+{\omega }_2+\vartheta \right)},(16)$

where $\begin{array}{l}{y}_1=\frac{\phi {\delta }_1{\omega }_5}{\gamma {\omega }_4}+\frac{\tau }{{\omega }_7}\left(\frac{{{\epsilon }_e\omega }_5}{\gamma }+{\epsilon }_i+\frac{{\epsilon }_m}{{\omega }_6}\left(\frac{\phi {\delta }_2{\omega }_5}{\gamma }+\alpha \right)\right),y_2=\sigma \vartheta {\omega }_4{\omega }_6+{\epsilon }_mŋ{\delta }_2\left(Z{\beta }_v+{\omega }_2\right),\\ y_3=\sigma \vartheta {\omega }_4{\omega }_6+{\epsilon }_mŋ{\delta }_2\left(Z{\beta }_v+{\omega }_2\right)-{\omega }_4{\omega }_6\left(Z{\beta }_v\left(\mu +ŋ\right)-Z^2{\beta }_s{\beta }_v-{\omega }_1{\omega }_2\right),{\text{and}y}_4=\frac{Z}{K{y}_3}\left(R_e\right)-1.\end{array}$

Note that, in evaluating $S^{*}$ and $V^{*},$

$y_4=\left(Z{\beta }_s{+\omega }_2)({\beta }_v+{\omega }_1\right){\omega }_4{\omega }_6-y_2,$

$y_4=\left(Z^2{\beta }_s{\beta }_v+Z{\beta }_s{\omega }_2+Z{\beta }_v{\omega }_1+{\omega }_1{\omega }_2\right){\omega }_4{\omega }_6-y_2.$

Substituting ${\omega }_1$ and ${\omega }_2$ and simplifying gives

$y_4=Z{\beta }_s\left(\mathrm{\mu }+\sigma \right)+Z{\beta }_v\vartheta +Z{\beta }_v\left(\mathrm{\mu }+ŋ\right)+Z^2{\beta }_s{\beta }_v+\left(\mathrm{\mu }+\vartheta +ŋ\right)\left(\mathrm{\mu }+\vartheta \right)-y_2,$

$y_4=Z{\beta }_s\left(\mathrm{\mu }+\vartheta \right)+Z{\beta }_v\vartheta -y_3.$

Let $K=\frac{1}{\left(\mathrm{\mu }+\mathrm{\sigma }+\vartheta \right)}\left(\frac{r_e}{{\omega }_3}+\frac{\gamma }{{\omega }_3{\omega }_5}+\frac{\phi r_q}{{\omega }_3{\omega }_4}+\frac{\left({\omega }_4\alpha \gamma +{\omega }_5\phi {\delta }_2\right)r_m}{{\omega }_3{\omega }_4{\omega }_5{\omega }_6}\right)$.

It follows that

$y_4=ZK\left({\beta }_s\left(\mathrm{\mu }+\sigma \right)+{\beta }_v\vartheta \right)-K{y}_3,$

$y_4=\frac{Z\left({\beta }_s\left(\mathrm{\mu }+\sigma \right)+{\beta }_v\vartheta \right)}{\left(\mathrm{\mu }+\mathrm{\sigma }+\vartheta \right)}\left(\frac{r_e}{{\omega }_3}+\frac{\gamma }{{\omega }_3{\omega }_5}+\frac{\phi r_q}{{\omega }_3{\omega }_4}+\frac{\left({\omega }_4\alpha \gamma +{\omega }_5\phi {\delta }_2\right)r_m}{{\omega }_3{\omega }_4{\omega }_5{\omega }_6}\right)-K{y}_3,$

$y_4=\frac{Z}{K{y}_3}\left(R_e\right)-1.(17)$

The expressions for $S^{*}$ and $V^{*}$ can be rewritten as follows

$S^{*}=\frac{\vartheta {\omega }_4{\omega }_6\left(Z{\beta }_v+{\omega }_2\right)\left(\mathrm{\Lambda }+y_1{I}^{*}\right)}{\frac{Z}{K{y}_3}\left(R_e\right)-1},$

$V^{*}=\frac{\vartheta {\omega }_4{\omega }_6\left(\mathrm{\Lambda }+y_1{I}^{*}\right)}{\frac{Z}{K{y}_3}\left(R_e\right)-1}.$

Let us consider the denominators of $S^{*}$ and $V^{*}$ in the following cases:

(a) If $R_e>1$, then $\frac{Z}{K{y}_3}\left(R_e\right)-1=\frac{Z}{K{y}_3}\left(R_e\right)-1>0.$ This results in the endemic equilibrium point $K^{*}=\left(S^{*},V^{*},E^{*},Q^{*},I^{*},M^{*},R^{*}\right)$ since the expressions for $S^{*}$ and $V^{*}$ are nonnegative, which implies $I^{*}>0$.

(b) If $R_e<1$, then since $\frac{Z}{K{y}_3}\left(R_e\right)-1=\frac{Z}{K{y}_3}\left(R_e\right)-1<0,$ which implies $I^{*}<0$, endemic equilibrium does not exist.

Using Eq. (17), for $I^{*}>0$ when $R_e>1,$ it implies that the endemic equilibrium $P^{*}$ has a positive and unique equilibrium point when $R_e>1.$

3.6 Global stability of endemic equilibrium

Let us consider a globally positively definite and unbounded function $L\left(x\right)$ with a globally negative time derivative. For $L\left(x\right)<0$ for all $x\ne x^{*},$ then at equilibrium, $x^{*}$ is globally stable for the autonomous system $x^{\prime }=f\left(x\right)$ and $L\left(x\right)$ is the Lyapunov function.

Theorem 3. The system (14) has an endemic equilibrium $P^{*}=\left(S^{*},V^{*},E^{*},Q^{*},I^{*},M^{*}\right)$ that satisfies $S^{*}>0,V^{*}>0,E^{*}>0,Q^{*}>0,I^{*}>0,\text{and\hspace{0.17em}}{M}^{*}>0$. If $R_e>1,$ the endemic equilibrium point globally asymptotically stable in a positive region $R_{+}^7$ on $H$ with

$\left\{\begin{array}{l}\mathit{\Lambda }=\mathit{Z}{\mathit{\beta }}_{\mathit{s}}{\mathit{S}}^{*}+{\mathit{\omega }}_1{\mathit{S}}^{*}-\mathit{\sigma }{\mathit{V}}^{*}-{\mathit{\delta }}_1{\mathit{Q}}^{*},\\ {\mathit{\omega }}_2=\frac{\mathit{\vartheta }{\mathit{S}}^{*}-\mathit{K}{\mathit{\beta }}_{\mathit{v}}{\mathit{V}}^{*}}{{\mathit{V}}^{*}},\\ {\mathit{\omega }}_3=\frac{\mathit{Z}\left({\mathit{\beta }}_{\mathit{s}}{\mathit{S}}^{*}+{{\mathit{\beta }}_{\mathit{v}}\mathit{V}}^{*}\right)}{{\mathit{E}}^{*}},\\ {\mathit{\omega }}_4=\frac{\mathit{ŋ}{\mathit{S}}^{*}+\mathit{\phi }{\mathit{E}}^{*}}{{\mathit{Q}}^{*}},\\ {\mathit{\omega }}_5=\frac{\mathit{\gamma }{\mathit{E}}^{*}}{{\mathit{I}}^{*}},\\ {\mathit{\omega }}_6=\frac{{\mathit{\delta }}_2{\mathit{Q}}^{*}+\mathit{\alpha }{\mathit{I}}^{*}}{{\mathit{M}}^{*}}.\end{array}\right.(18)$

Proof. We define the Lyapunov function $L$ and its endemic equilibrium as follows:

$\begin{array}{c}L=k_i\left(\left(S-S^{*}-S^{*}\mathit{ln}\frac{S}{S^{*}}\right)|+|\left(V-V^{*}-V^{*}\mathit{ln}\frac{V}{V^{*}}\right)\right.\\ +\left(E-E^{*}-E^{*}\mathit{ln}\frac{E}{E^{*}}\right)+\left(Q-Q^{*}-Q^{*}\mathit{ln}\frac{Q}{Q^{*}}\right)\\ +\left.\left(I-I^{*}-I^{*}\mathit{ln}\frac{I}{I^{*}}\right)|+|\left(M-M^{*}-M^{*}\mathit{ln}\frac{M}{M^{*}}\right)\right),\end{array}(19)$

from Eq. (19), we have $k_i>0\text{\hspace{0.17em}and\hspace{0.17em}}i=1,2,3,4,5,6.$

$\frac{dL}{dt}=k_i\left(\left(1-\frac{S^{*}}{S}\right)S^{\prime }+\left(1-\frac{V^{*}}{V}\right)V^{\prime }+\left(1-\frac{E^{*}}{E}\right)E^{\prime }+\left(1-\frac{Q^{*}}{Q}\right)Q^{\prime }+\left(1-\frac{I^{*}}{I}\right)I^{\prime }+\left(1-\frac{M^{*}}{M}\right)M^{\prime }\right).(20)$

Substituting (2) and (18) into (20) gives

$\begin{array}{c}\frac{dL}{dt}=k_1\left(1-\frac{S^{*}}{S}\right)\left(Z{\beta }_s{S}^{*}+{\omega }_1{S}^{*}-\sigma V^{*}-{\delta }_1{Q}^{*}+\sigma V+{\delta }_{1}Q-{\rho }_{1}S-{\omega }_{1}S\right)\\ +k_2\left(1-\frac{V^{*}}{V}\right)\left(\vartheta S-Z{\beta }_{v}V-\left(\frac{\vartheta S^{*}-Z{\beta }_v{V}^{*}}{V^{*}}\right)V\right)\\ +k_3\left(1-\frac{E^{*}}{E}\right)\left(Z{\beta }_{s}S+Z{\beta }_{v}V-\left(\frac{Z\left({\beta }_s{S}^{*}+{{\beta }_{v}V}^{*}\right)}{E^{*}}\right)E\right)\\ +k_4\left(1-\frac{Q^{*}}{Q}\right)\left(ŋS+\phi E-\frac{ŋS^{*}+\phi E^{*}}{Q^{*}}\right)Q+k_5\left(1-\frac{I^{*}}{I}\right)\\ \left(\gamma E-\frac{\gamma E^{*}}{I^{*}}\right)I+k_6\left(1-\frac{M^{*}}{M}\right)\left({\delta }_{2}Q+\alpha I-\frac{{\delta }_2{Q}^{*}+\alpha I^{*}}{M^{*}}\right)M.\end{array}(21)$

Solving the (21) value gives

$\begin{array}{c}\frac{dL}{dt}=k_1{\omega }_1\left(2-\frac{S}{S^{*}}-\frac{S^{*}}{S}\right)S^{*}+k_{1}Z{\beta }_s{S}^{*}\left(3-\frac{S^{*}}{S}-\frac{E}{E^{*}}-\frac{S{E}^{*}}{S^{*}E}\right)\\ +k_4ŋS^{*}\left(3-\frac{S}{S^{*}}-\frac{Q}{Q^{*}}-\frac{S^{*}{Q}^{*}}{SQ}\right)+k_1\sigma V^{*}\left(3-\frac{S}{S^{*}}-\frac{V^{*}}{V}-\frac{S^{*}V}{S{V}^{*}}\right)\\ +k_1{\delta }_1{Q}^{*}\left(3-\frac{S}{S^{*}}-\frac{Q^{*}}{Q}-\frac{S^{*}Q}{S{Q}^{*}}\right)\\ +k_2{k}_{3}Z{\beta }_v{V}^{*}\left(3-\frac{V^{*}}{V}-\frac{E}{E^{*}}-\frac{V{E}^{*}}{V^{*}E}\right)+k_4\phi E^{*}\left(3-\frac{E^{*}}{E}-\frac{Q}{Q^{*}}-\frac{E{Q}^{*}}{E^{*}Q}\right)\\ +k_6\alpha I^{*}\left(3-\frac{I^{*}}{I}-\frac{M}{M^{*}}-\frac{I{M}^{*}}{I^{*}M}\right)+k_5\gamma E^{*}\left(3-\frac{E^{*}}{E}-\frac{I}{I^{*}}-\frac{E{I}^{*}}{E^{*}I}\right)\\ +k_6{\delta }_2{Q}^{*}\left(3-\frac{Q^{*}}{Q}-\frac{M}{M^{*}}-\frac{Q{M}^{*}}{Q^{*}M}\right).\end{array}(22)$

Let $f_1=\frac{S}{S^{*}},f_2=\frac{V}{V^{*}},f_3=\frac{E}{E^{*}},f_4=\frac{Q}{Q^{*}},f_5=\frac{I}{I^{*}},\text{and\hspace{0.17em}}{f}_6=\frac{M}{M^{*}}.$

$k_1,k_2,k_3,k_4,k_5,\text{and\hspace{0.17em}}{k}_6$ are obtained by setting the coefficients of $f_1,f_2,f_3,f_4,f_1{f}_4\text{\hspace{0.17em}and\hspace{0.17em}}{f}_2{f}_4$ equal to 0, which after solving gives

$k_1=k_3=1,k_2=\frac{\sigma V^{*}}{\vartheta S^{*}-Z{\beta }_v{V}^{*}},k_4=\frac{{\delta }_1{Q}^{*}}{ŋS^{*}+\phi E^{*}},k_5=\frac{Z\left({\beta }_s{S}^{*}+{{\beta }_{v}V}^{*}\right)}{\gamma E^{*}},\text{and\hspace{0.17em}}{k}_6=\frac{{\delta }_1{Q}^{*}}{{\delta }_2{Q}^{*}+\alpha I^{*}}.$

Hence, the result is

$\begin{array}{c}\frac{dL}{dt}={\omega }_1\left(2-f_1-\frac{1}{f_1}\right)S^{*}+Z{\beta }_s{S}^{*}\left(3-\frac{1}{f_1}-f_3-\frac{f_1}{f_3}\right)\\ +{\sigma V}^{*}\left(3-f_1-\frac{1}{f_2}-\frac{f_2}{f_1}\right)+{\delta }_1{Q}^{*}\left(3-f_1-\frac{1}{f_4}-\frac{f_4}{f_1}\right)\\ +k_{2}Z{\beta }_v{V}^{*}\left(3-\frac{1}{f_2}-f_3-\frac{f_2}{f_3}\right)+k_4\phi E^{*}\left(3-\frac{1}{f_3}-f_4-\frac{f_3}{f_4}\right)\\ +k_5\gamma E^{*}\left(3-\frac{1}{f_3}-f_5-\frac{f_3}{f_5}\right)+k_6\alpha I^{*}\left(3-\frac{1}{f_5}-f_6-\frac{f_5}{f_6}\right)\\ +k_6{\delta }_2{Q}^{*}\left(3-\frac{1}{f_5}-f_6-\frac{f_5}{f_6}\right).\end{array}(23)$

From the above Eqs (22, 23), it can be realized that if the arithmetic mean is greater than or equal to their geometric mean, then $\frac{dL}{dt}\le 0.\frac{dL}{dt}=0$ only holds if $f_1{=f}_2=f_3=f_4=f_5=f_6.$ It implies that ${S=S}^{*},{V=V}^{*},E=E^{*},Q=Q^{*},\text{and\hspace{0.5em}}I=I^{*}=M=M^{*}\text{in\hspace{0.17em}}H$; therefore, the largest invariant set is $\left\{S,V,E,Q,I,M\in H:\frac{dL}{dt}=0\right\}.$ The endemic equilibrium $P^{*}$ is globally asymptotically stable in the positive region $R_{+}^7$ when $R_e>1$ based on the Lyapunov–LaSalle stability theorem.

4 Model parameterization and sensitivity analysis

We evaluate the reproduction number to find the transmission routes of the pandemic and the proportion of individuals. A sensitivity analysis was carried out to evaluate the contribution of each parameter to the reproduction number.

4.1 Estimation of model parameters

We evaluate the reproduction numbers using the values in Table 3. The effective reproduction number $R_e=0.85575$ and a proportion of 0.19998 of the exposed individuals move to the infected class at any given time. The contribution from each transmission route is shown in the table below.

Table 2 reveals that the infectious and exposed classes contribute significantly to the transmission of the disease. This implies that there must be admissible measures to halt the transmission of the disease.

Table 2

Table 2. Contribution of transmission routes to the effective reproduction number ${\mathit{R}}_{\mathit{e}}$.

4.2 Sensitivity analysis

We illustrate the arithmetic behavior of model (2) parameters’ sensitivity. The sensitivity of the parameters in $R_e$ for the model is defined in Eq. 24 as

${Ƥ}_p^{R_e}=\frac{\partial R_e}{\partial p}.\frac{p}{R_0}\approx \frac{\%∆R_e}{\%∆p}.(24)$

$\begin{array}{c}{Ƥ}_{\mathrm{\Lambda }}^{R_e}=1,{Ƥ}_{{\mathrm{\beta }}_s}^{R_e}=\frac{{\mathrm{\beta }}_s\left(\mu +\sigma \right)}{a_1},{Ƥ}_{{\mathrm{\beta }}_v}^{R_e}=\frac{{\mathrm{\beta }}_v\mathrm{\vartheta }}{a_1},{Ƥ}_{r_e}^{R_e}=R_{eE}{a}_2,{Ƥ}_{r_q}^{R_e}=R_{eQ}{a}_2,\\ {Ƥ}_{r_m}^{R_e}=R_{eM}{a}_2,{Ƥ}_{{\epsilon }_e}^{R_e}=-{\epsilon }_e\left(\frac{R_{eE}}{{\omega }_3}+\frac{R_{eI}}{{\omega }_3}+\frac{R_{eQ}}{{\omega }_3}\right),{Ƥ}_{{\epsilon }_i}^{R_e}=-{\epsilon }_i\left(\frac{R_{eI}}{{\omega }_5}+\frac{R_{eQ}}{{\omega }_5}-\frac{\phi {\delta }_2{r}_m}{R_{eM}}\right),\\ {Ƥ}_{{\epsilon }_m}^{R_e}=-\left(\frac{R_{eM{\epsilon }_m}}{{\omega }_6}\right)\left(R_{eE}+R_{eI}+R_{eM}+\frac{\phi r_q}{{\omega }_3{\omega }_4}\right),{Ƥ}_{{\delta }_1}^{R_e}=-{\delta }_1\left(\frac{R_{eQ}}{{\omega }_4}+\frac{R_{eM}}{{\omega }_4}-\frac{\mathrm{\alpha }\gamma r_m}{a_2}\right),\\ {Ƥ}_{{\delta }_2}^{R_e}=-{\delta }_2\left(\frac{R_{eQ}}{{\omega }_4}+\frac{R_{eM}}{{\omega }_4}-\frac{\left(\mathrm{\alpha }\gamma +\phi {\omega }_5\right)r_m}{a_2}\right),{Ƥ}_{\mathrm{\alpha }}^{R_e}=-\mathrm{\alpha }\left(\frac{R_{eI}}{{\omega }_5}+\frac{R_{0M}}{{\omega }_5}-\frac{R_{eM}}{a_2}\right),\\ {Ƥ}_{\gamma }^{R_e}=-\gamma \left(\frac{R_{eE}}{{\omega }_5}+\frac{R_{eI}}{{{\omega }_3\omega }_5}+\frac{R_{eQ}}{{\omega }_3^2{\omega }_5}+R_{eM}-\frac{1}{{{\mathrm{\beta }}_s{\omega }_{2}a}_1}\right),{Ƥ}_{\mathrm{\mu }}^{R_e}=-\mu a_1\left(\frac{R_{eE}}{{\omega }_3}+\frac{R_{0I}}{{\omega }_3}+\frac{\gamma }{{\omega }_3{\omega }_5^2}+\frac{\phi r_q}{{\omega }_3^2{\omega }_4}+\frac{\phi r_q}{{\omega }_3{\omega }_4^2}-{\mathrm{\beta }}_s{a}_2\right),\\ {Ƥ}_{\sigma }^{R_e}=-\left({\mathrm{\beta }}_s\sigma \right)\left(R_{eE}+R_{eI}+R_{eQ}+\frac{R_{eM}}{a_3^2}-\frac{R_{eM}}{a_3^2}\right),{Ƥ}_{\mathrm{\vartheta }}^{R_e}=-\left({\mathrm{\beta }}_v\mathrm{\vartheta }\right)\left(R_{eE}+R_{eI}+R_{eQ}+\frac{R_{eM}}{a_3}-\frac{1}{a_2}\right),\\ \text{and\hspace{0.17em}}{Ƥ}_{\phi }^{R_e}=-\phi \left(\frac{R_{eE}}{{\omega }_3}+\frac{R_{eI}}{{\omega }_3}+\frac{R_{eQ}}{{\omega }_3^2{\omega }_4}+\frac{R_{eM}}{{\omega }_3{a}_3}-\frac{r_q}{{{\omega }_3\omega }_5}-\frac{{\delta }_2{r}_m}{{{\omega }_{3}a}_2}\right),\end{array}$

where $a_1={\mathrm{\beta }}_s{\omega }_2+{\mathrm{\beta }}_v\mathrm{\vartheta },a_2=R_{eE}+R_{eI}+R_{eQ}+R_{eM},\text{and}{a}_3=\mu +\sigma +\vartheta .$

From the above analysis, it is observed that $R_e$ will increase proportionally as the following parameters $\left({\mathrm{\beta }}_s,{\mathrm{\beta }}_v,r_e,r_q,r_m\right)$ increase; on the other hand, $R_0$ decreases proportionally as the following parameters $\left(\mathrm{\alpha },\gamma ,\sigma ,\phi ,\mathrm{\mu },\mathrm{\vartheta },{\epsilon }_e,{\epsilon }_i,{\epsilon }_m,{\delta }_1,{\delta }_2\right)$ increase since sensitivity indices are sign determined. This can be used to show the numerical importance of the various parameters in $R_e$. Table 3 illustrates the numerical results of the sensitivity index of the various parameters in $R_e$.

Table 3

Table 3. Values of the sensitivity index of parameters in ${\mathit{R}}_{\mathit{e}}$.

Based on the numerical illustrations shown in Table 3, the effective reproduction number $R_e$ will change as the parameters change. Sensitivity analysis is interpreted based on the sign associated with the particular parameter. The effective reproduction number $R_e$ will decrease as the parameter values with the negative sign increase, while it increases when they decrease. Considering that ${\mathrm{\beta }}_s=0.8192,$ it indicates that ${\mathrm{\beta }}_s$ will increase $R_e$ by 81% whenever there is a 1% increment in the transmission rate of the susceptible individuals. This depicts the explanation for all the parameter values with a positive sign. On the other hand, ${\epsilon }_i=-0.1546$ indicates that ${\epsilon }_i$ will decrease $R_e$ by 15% whenever there is a 1% increment in the recovery rate of the infected individuals. This depicts the explanation for all the parameter values with a negative sign.

Now, we examine the relationship between the parameters of effective reproduction number $R_e$ by checking the partial rank correlation coefficients (PRCCs) of the parameters. The following figure illustrates the behavioral pattern of each parameter in the transmission dynamics of the disease.

From Figure 2, one could realize that $\left({\mathrm{\beta }}_s,{\mathrm{\beta }}_v,r_e,r_q,\gamma ,\phi \right)$ have a high positive effect on $R_e,$ while $\left(\mathrm{\mu },\mathrm{\vartheta },{\epsilon }_e\right)$ have a high negative effect on $R_e.$ One could realize that all these parameters are associated with vaccination and interventions to mitigate the infectivity rate in the population.

Figure 2

Figure 2. Partial rank correlation coefficient (PRCC) of effective reproduction number ${\mathit{R}}_{\mathit{e}}$ parameters.

4.3 Effects of vaccination and infectivity reductions on the effective reproduction number ${\mathit{R}}_{\mathit{e}}$

This subsection explores the dynamics of $R_e$ with respect to vaccination and infectivity reduction parameters $\left(\vartheta ,\sigma ,r_e,r_q,r_m\right)$. The following figure demonstrates the changing effects of vaccination and infectivity reduction parameters on $R_e.$

Figure 3A demonstrates the influence of vaccination rate $\vartheta$ and vaccine inefficacy $\sigma .$ One can realize that the effective reproduction number $R_e$ increases as the vaccination rate $\vartheta$ decreases, while $R_e$ decreases as vaccine inefficacy $\sigma$ decreases. Figure 3B–g show the effects of vaccination rate $\vartheta$, vaccine inefficacy $\sigma$, and infectivity reductions. It is observed that to mitigate the disease’s transmission, there should be effective control measures to minimize the infectivity levels and maximize the vaccination rate to halt the disease’s transmission. Figures 3H–J show the effect of the infectivity reductions on $R_e.$ It could be observed that an increase in $r_e$ at any given time results in an increase in $R_e$ (Figure 3H; Figure 3I). It can be seen that $r_e$ and $r_q$ have a significant effect on $R_e$, as demonstrated in Figure 2. Based on the above graphical representation, it is important to minimize the infectivity levels of the exposed class $r_e$ and quarantined class $r_q$ to halt the disease’s transmission.

Figure 3

Figure 3. Effect of vaccination and infectivity reductions on the effective reproduction number ${\mathit{R}}_{\mathit{e}}.$ (A) a plot of σ and θ, (B) r_e and θ, (C) r_q and θ, (D) r_m and θ, (E) r_e and σ, (F) r_q and σ, (G) r_m and σ, (H) r_q and r_e (I) r_m and r_e (J) r_m and r_q. The meaning of the respective parameters is in Table 1.

5 Numerical simulations and discussions

In this section, the real data are observed to depict the situation of the transmission. The real data are compared with the model’s solution accuracy. The proposed model is applied to explore the transmission of COVID-19 disease in Ghana using data from the WHO [17]. This includes the daily number of infections and the total number of daily infections after the introduction of vaccines. The data are analyzed and compared to other literature works for the simulations. The daily infections and the total number of infections are simulated to analyze the behavioral pattern of the transmission of the disease. In particular, we illustrate the changing effects of the vaccination rate $\vartheta$, vaccine inefficacy $\sigma$ , and infectivity reductions on the size of infectious individuals.

5.1 Application of the $\mathit{S}\mathit{V}\mathit{E}\mathit{Q}\mathit{I}\mathit{M}\mathit{R}$ model

Once the model is formulated, it is necessary to compare it with data to check its validity. Here, we want to check the accuracy and authenticity of the model by verifying the extent to which the model can represent the real situation, as described in [18]. The simulation consists of the application of the data from Ghana to illustrate the transmission of COVID-19 for the period 1st January 2022 to 1st March 2022 when the individuals were vaccinated. The following state variables are considered using data from Ghana for the period 1st January 2022 to 1st March 2022 [17]. The estimated total population of Ghana is $31732129$ [19]; hence, $N\left(0\right)=31732129,$ and the assumed initial values are as follows: $S\left(0\right)=200000,V\left(0\right)=120000,E\left(0\right)=150000,I\left(0\right)=997,Q\left(0\right)=1000,M\left(0\right)=800,\text{and\hspace{0.17em}}R\left(0\right)=500.$ All the parameters used for the simulations are shown in Table 3.

The SVEQIMR model (fitted) depicts the pattern of the real situation (real data), as shown in Figure 4A, which represents the plot of the daily number of infections, and Figure 4B represents the total number of infections for the period of 1st January 2022 to 1st March 2022. The plot of the model results indicates the pattern of the real situation. That is, the fitted model and real situation agree with each other and illustrate the transmission dynamics of the spread of the pandemic.

Figure 4

Figure 4. Plot of (A) the daily cases and (B) the total number of cases of the fitted model verses the real data from 1st January to 1st March 2022.

5.2 Effects of vaccination and infectivity reductions on COVID-19 incidences

Here, we analyze the influence of vaccination rate and vaccine inefficacy on the spread of COVID-19 disease using the same dataset. Let us consider the following parameters: vaccination rates $\vartheta =\left(0,0.15,0.25,0.35,0.45\right)$ and vaccine inefficacy $\sigma =\left(0,0.02,0.05,0.07,0.09\right)$ for the infected class.

Again, we demonstrate the changing effects of the different infectivity reductions on the behavioral pattern of the transmission of COVID-19 disease by considering the following arbitrary values $\left(0.05,0.10,0.15,0.20,0.25\right)$ for the infectivity reduction of exposed individuals $r_e$, quarantined individuals $r_q$, and isolated individuals $r_m$. The results are presented in Figure 5.

Figure 5

Figure 5. Plot of (A) the impact of vaccination rate $\mathit{\vartheta }$ on COVID-19 daily cases (B) and the impact of vaccine inefficacy $\mathit{\sigma }$, (C) infectivity reduction of exposed individuals ${\mathit{r}}_{\mathit{e}}$, (D) infectivity reduction of quarantined individuals ${\mathit{r}}_{\mathit{q}}$, and (E) infectivity reduction of isolated individuals ${\mathit{r}}_{\mathit{m}}$ on COVID-19 daily cases.

From Figure 5A, it is observed that an increase in the vaccination rate $\vartheta$ results in a significant decrease in the daily number of reported cases. On the contrary, an increase in vaccine inefficacy $\sigma$ results in an increase in the daily number of reported cases, as shown in Figure 5B. It could also be observed that the disease dies out gradually with time as vaccine inefficacy $\sigma$ approaches zero (0); however, a high vaccination rate $\vartheta$ reduces the number of infections. Therefore, there is a significant effect of different values of vaccination rate $\vartheta$ and vaccine inefficacy $\sigma$ on the disease’s spread. The above illustration shows that vaccination against the spread of the coronavirus is very important and must be adhered to.

From Figure 5C, it is realized that an increase in the infectivity reduction of exposed individuals $r_e$ results in a significant increase and higher peaks in the daily reported cases. Furthermore, an increase in the infectivity reduction of the quarantined class $r_q$ results in a moderate increase in the daily reported cases compared to the peaks of $r_e$, as shown in Figure 5D. The peaks of daily reported cases are lower than those of the exposed class because the individuals in this group are confined and monitored. In Figure 5E, there is a significant decline in the daily reported cases as the infectivity reduction of the isolated individuals $r_m$ decreases, which has fewer peaks than the others. From the above graphical representation, the infectivity reduction of the exposed individuals $r_e$ should be controlled carefully because the higher the infectivity reduction, the higher the number of daily reported cases and vice versa. This graphical representation can be confirmed by the numerical illustration of the contribution of the transmission route in Table 2. The exposed compartment contributed significantly to the effective reproduction number $R_e$, with a proportion of 0.19998 of the individuals moving to the infected class. This means that approximately 20% of individuals in the exposed compartment become infected at a given time. The above diagram depicts the order of contribution of the transmission route to the effective reproduction number $R_e.$

6 Optimal control problem and cost–benefit analysis

We modify Equation 2 with the following optimal control variables: ${\kappa }_p\left(t\right),$ representing the public awareness of the prevalence of COVID-19 and related infections where ${\kappa }_p\left(t\right)\in \left[0,1\right]$, which reduces the forces of infection ${\rho }_1$ and ${\rho }_2$, by $1-{\kappa }_p\left(t\right)$; ${\kappa }_v\left(t\right),$ denoting the control effort to intensify COVID-19 vaccination; and ${\kappa }_i\left(t\right),$ denoting the control effort for infectivity reduction, which reduces infectivity by $1-{\kappa }_i\left(t\right)$. All these efforts denote the control interventions in minimizing the transmission of COVID-19 and its reinfections. The modified Eq. (2) is given in Eq. (25) as follows:

$\left\{\begin{array}{l}\frac{\mathit{d}\mathit{S}}{\mathit{d}\mathit{t}}=\frac{\mathit{d}\mathit{S}}{\mathit{d}\mathit{t}}=\mathit{\Lambda }+\mathit{\sigma }\mathit{V}+{\mathit{\delta }}_1\mathit{Q}+\mathit{\tau }\mathit{R}-\left(1-{\mathit{\kappa }}_{\mathit{p}}\right){\mathit{\rho }}_1\mathit{S}-\left(\mathrm{\mu }+{\mathit{\kappa }}_{\mathit{v}}\mathit{\vartheta }+\mathit{ŋ}\right)\mathit{S},\\ \frac{\mathit{d}\mathit{V}}{\mathit{d}\mathit{t}}={\mathit{\kappa }}_{\mathit{v}}\mathit{\vartheta }\mathit{S}-\left(1-{\mathit{\kappa }}_{\mathit{p}}\right){\mathit{\rho }}_2\mathit{V}-\left(\mathrm{\mu }+\mathrm{\sigma }\right)\mathit{V},\\ \frac{\mathit{d}\mathit{E}}{\mathit{d}\mathit{t}}=\left(1-{\mathit{\kappa }}_{\mathit{p}}\right)\left({\mathit{\rho }}_1\mathit{S}+{\mathit{\rho }}_2\mathit{V}\right)-\left(\mathrm{\mu }+\left(1-{\mathit{\kappa }}_{\mathit{i}}\right)\mathit{\phi }+\left(1-{\mathit{\kappa }}_{\mathit{i}}\right)\mathit{\gamma }+{\mathit{\epsilon }}_{\mathit{e}}\right)\mathit{E},\\ \frac{\mathit{d}\mathit{Q}}{\mathit{d}\mathit{t}}=\mathit{ŋ}\mathit{S}+\left(1-{\mathit{\kappa }}_{\mathit{i}}\right)\mathit{\phi }\mathit{E}-\left(\mathrm{\mu }+{\mathit{\delta }}_1+{\left(1-{\mathit{\kappa }}_{\mathit{i}}\right)\mathit{\delta }}_2\right)\mathit{Q},\\ \frac{\mathit{d}\mathit{I}}{\mathit{d}\mathit{t}}=\left(1-{\mathit{\kappa }}_{\mathit{i}}\right)\mathit{\gamma }\mathit{E}-\left(\mathrm{\mu }+\mathit{\alpha }+{\mathit{\epsilon }}_{\mathit{i}}\right)\mathit{I},\\ \frac{\mathit{d}\mathit{M}}{\mathit{d}\mathit{t}}=\left(1-{\mathit{\kappa }}_{\mathit{i}}\right){\mathit{\delta }}_2\mathit{Q}+\mathit{\alpha }\mathit{I}-\left(\mathrm{\mu }+{\mathit{\epsilon }}_{\mathit{m}}\right)\mathit{M},\\ \frac{\mathit{d}\mathit{R}}{\mathit{d}\mathit{t}}={\mathit{\epsilon }}_{\mathit{e}}\mathit{E}+{\mathit{\epsilon }}_{\mathit{i}}\mathit{I}+{\mathit{\epsilon }}_{\mathit{m}}\mathit{M}-\left(\mathrm{\mu }+\mathit{\tau }\right)\mathit{R},\end{array}\right.(25)$

where ${\rho }_1$ and ${\rho }_2$ remain as in Section 2 with the initial conditions given in the model (2).

6.1 Objective functional

We now formulate the optimal trajectories that show the effect of the control efforts ${\kappa }_p\left(t\right),{\kappa }_v\left(t\right),\text{and}{\kappa }_i\left(t\right)$ subjected to (36); the objective functional $Q$ is given as

$Q\left({\kappa }_p,{\kappa }_v,{\kappa }_i\right)={\int }_0^{t_f}\left[h_{1}E+h_{2}I+h_{3}Q+h_{4}M+\frac{1}{2}{z}_1{\kappa }_p^2\left(t\right)+\frac{1}{2}{z}_2{\kappa }_v^2\left(t\right)+\frac{1}{2}{z}_3{\kappa }_i^2\left(t\right)\right]dt.(26)$

We focus on minimizing the cost function (26), and the total cost of implementing the optimal control is given as

$Z={\int }_0^{t_f}\left[\frac{1}{2}{z}_1{\kappa }_p^2\left(t\right)+\frac{1}{2}{z}_2{\kappa }_v^2\left(t\right)+\frac{1}{2}{z}_3{\kappa }_i^2\left(t\right)\right]dt.(27)$

The parameters $z_1,z_2,$and $z_3$ in Eq. 27 are the balancing cost factors for ${\kappa }_p\left(t\right),{\kappa }_v\left(t\right),{\text{and}\kappa }_i\left(t\right)$, respectively. All the control efforts ${\kappa }_p\left(t\right),{\kappa }_v\left(t\right),\text{and\hspace{0.17em}}{\kappa }_i\left(t\right)$ are assumed to be bounded by Lebesgue measurable time-dependent functions on the interval $\left[0,t_f\right]$, where $t_f$ is the final time with the control effort set defined in Eq. 28 as

$\mathrm{\Gamma }=\left({\kappa }_p,{\kappa }_v,{\kappa }_i\text{\hspace{0.17em}for\hspace{0.17em}}\left(0\le {\kappa }_p,{\kappa }_v,{\kappa }_i\le 1,\right),0\le t\le t_f\right).(28)$

Now, we establish point-wise Hamiltonian $H$ through Pontryagin’s maximum principle to transform the optimal control system (25) and its associated objective functional (26). The following optimal solution is achieved.

$\begin{array}{ll}\hfill H& =h_{1}E+h_{2}I+h_{3}Q+h_{4}M+\frac{1}{2}{z}_1{\kappa }_p^2\left(t\right)+\frac{1}{2}{z}_2{\kappa }_v^2\left(t\right)+\frac{1}{2}{z}_3{\kappa }_i^2\left(t\right)+{\lambda }_S\left(\mathrm{\Lambda }+\sigma V+{\delta }_{1}Q+\tau R-\left(1-{\kappa }_p\right){\rho }_{1}S-\left(\mathrm{\mu }+{\kappa }_v\vartheta +ŋ\right)S\right)\hfill \\ \hfill & +{\lambda }_V\left({\kappa }_v\vartheta S-\left(1-{\kappa }_p\right){\rho }_{2}V-\left(\mathrm{\mu }+\mathrm{\sigma }\right)V\right)+{\lambda }_E\left(\left(1-{\kappa }_p\right)\left({\rho }_{1}S+{\rho }_{2}V\right)-\left(\mathrm{\mu }+\left(1-{\kappa }_i\right)\phi +\gamma +{\epsilon }_e\right)E\right)\hfill \\ \hfill & +{\lambda }_Q\left(ŋS+\left(1-{\kappa }_i\right)\phi E-\left(\mathrm{\mu }+{\delta }_1+{\left(1-{\kappa }_i\right)\delta }_2\right)Q\right)+{\lambda }_I\left(\gamma E-\left(\mathrm{\mu }+\alpha +{\epsilon }_i\right)I\right)\hfill \\ \hfill & +{\lambda }_M\left(\left(1-{\kappa }_i\right){\delta }_{2}Q+\alpha I-\left(\mathrm{\mu }+{\epsilon }_m\right)M\right)+{\lambda }_R\left({\epsilon }_{e}E+{\epsilon }_{i}I+{\epsilon }_{m}M-\left(\mathrm{\mu }+\tau \right)R\right),\hfill \end{array}(29)$

where ${\lambda }_S,{\lambda }_V,{\lambda }_E,{\lambda }_Q,{\lambda }_I,{\lambda }_M,\text{and\hspace{0.17em}}{\lambda }_R$ in Eq. 29 are the co-state variables with respect to the state variables, $S,V,E,Q,I,M,\text{and\hspace{0.17em}}R.$

Theorem 4. Given ${\kappa }_p^{*}\left(t\right),{\kappa }_v^{*}\left(t\right),\text{and\hspace{0.17em}}{\kappa }_i^{*}\left(t\right)$ as the optimal controls and the corresponding solutions $S^0,V^0,E^0,Q^0,I^0,M^0,\text{and\hspace{0.17em}}{R}^0$ of the system (25), which minimizes $Z\left({\kappa }_p\left(t\right),{\kappa }_v\left(t\right),{\kappa }_i\left(t\right)\right)$ over $\mathrm{\Gamma },$ then there exist co-state variables ${\lambda }_S,{\lambda }_V,{\lambda }_E,{\lambda }_Q,{\lambda }_I,{\lambda }_M,{\text{and}\lambda }_R$ that satisfy

$\frac{d{\lambda }_j}{dt}=-\frac{\partial H}{\partial j},(30)$

with conditions ${\lambda }_j\left(t_f\right)=0,$ where $j=S,V,E,Q,I,M,R.$ Then, the optimality conditions that minimize the Hamiltonian, H, of (29) with respect to the controls are given as

$\left\{\begin{array}{l}{\mathit{\kappa }}_{\mathit{p}}^{*}\left(\mathit{t}\right)=\min \left\{{\mathit{\kappa }}_{\mathit{p}}\max ,\max \left(0,\frac{\left({\mathit{\lambda }}_{\mathit{E}}-{\mathit{\lambda }}_{\mathit{S}}\right){\mathit{\rho }}_1^0{\mathit{S}}^0+\left({\mathit{\lambda }}_{\mathit{E}}-{\mathit{\lambda }}_{\mathit{V}}\right){\mathit{\rho }}_2^0{\mathit{V}}^0}{{\mathit{z}}_1}\right)\right\},\\ {\mathit{\kappa }}_{\mathit{v}}^{*}\left(\mathit{t}\right)=\min {\left\{{\mathit{\kappa }}_{\mathit{v}}\max ,\max \left(0,\frac{\left({\mathit{\lambda }}_{\mathit{S}}-{\mathit{\lambda }}_{\mathit{V}}\right)\mathit{\vartheta }{\mathit{S}}^0}{{\mathit{z}}_2}\right)\right\}}^{\prime }\\ {\mathit{\kappa }}_{\mathit{i}}^{*}\left(\mathit{t}\right)=\min \left\{{\mathit{\kappa }}_{\mathit{i}}\max ,\max \left(0,\frac{\left({\mathit{\lambda }}_{\mathit{Q}}-{\mathit{\lambda }}_{\mathit{E}}\right)\mathit{\phi }{\mathit{E}}^0+\left({\mathit{\lambda }}_{\mathit{I}}-{\mathit{\lambda }}_{\mathit{E}}\right)\mathit{\gamma }{\mathit{E}}^0+\left({\mathit{\lambda }}_{\mathit{M}}-{\mathit{\lambda }}_{\mathit{Q}}\right){\mathit{\delta }}_2{\mathit{Q}}^0}{{\mathit{z}}_3}\right)\right\}.\end{array}\right.(31)$

Proof: We formulate the adjoint equation for the optimal system by taking the partial derivative of Eq. 30 as follows in Eq. (32):

$\left\{\begin{array}{l}\frac{\mathit{d}{\mathit{\lambda }}_{\mathit{S}}}{\mathit{d}\mathit{S}}=\left(1-{\mathit{\kappa }}_{\mathit{p}}\right)\left({\mathit{\lambda }}_{\mathit{S}}-{\mathit{\lambda }}_{\mathit{E}}\right){\mathit{\rho }}_1+\left({\mathit{\lambda }}_{\mathit{V}}-{\mathit{\lambda }}_{\mathit{S}}\right){\mathit{\kappa }}_{\mathit{v}}\mathit{\vartheta }+\left({\mathit{\lambda }}_{\mathit{S}}-{\mathit{\lambda }}_{\mathit{V}}\right)\mathit{\sigma }\mathit{V}+\left({\mathit{\lambda }}_{\mathit{Q}}-{\mathit{\lambda }}_{\mathit{S}}\right)\mathit{ŋ}+\left({\mathit{\lambda }}_{\mathit{S}}-{\mathit{\lambda }}_{\mathit{Q}}\right){\mathit{\delta }}_1\mathit{Q}+\left({\mathit{\lambda }}_{\mathit{S}}-{\mathit{\lambda }}_{\mathit{R}}\right)\mathit{\tau }\mathit{R}+\mathit{\mu }{\mathit{\lambda }}_{\mathit{S}},\\ \frac{\mathit{d}{\mathit{\lambda }}_{\mathit{V}}}{\mathit{d}\mathit{V}}=\left(1-{\mathit{\kappa }}_{\mathit{p}}\right)\left({\mathit{\lambda }}_{\mathit{E}}-{\mathit{\lambda }}_{\mathit{V}}\right){\mathit{\rho }}_2+\left({\mathit{\lambda }}_{\mathit{V}}-{\mathit{\lambda }}_{\mathit{S}}\right){\mathit{\kappa }}_{\mathit{v}}\mathit{\vartheta }\mathit{S}+\left({\mathit{\lambda }}_{\mathit{S}}-{\mathit{\lambda }}_{\mathit{V}}\right)\mathit{\sigma }+\mathit{\mu }{\mathit{\lambda }}_{\mathit{V}},\\ \frac{\mathit{d}{\mathit{\lambda }}_{\mathit{E}}}{\mathit{d}\mathit{E}}=-{\mathit{h}}_1+\left(1-{\mathit{\kappa }}_{\mathit{p}}\right)\left[\left({\mathit{\lambda }}_{\mathit{E}}-{\mathit{\lambda }}_{\mathit{S}}\right){\mathit{\rho }}_1\mathit{S}+\left({\mathit{\lambda }}_{\mathit{E}}-{\mathit{\lambda }}_{\mathit{V}}\right){\mathit{\rho }}_2\mathit{V}\right]+\left({\mathit{\lambda }}_{\mathit{Q}}-{\mathit{\lambda }}_{\mathit{E}}\right)\left(1-{\mathit{\kappa }}_{\mathit{i}}\right)\mathit{\phi }+\left({\mathit{\lambda }}_{\mathit{I}}-{\mathit{\lambda }}_{\mathit{E}}\right)\left(1-{\mathit{\kappa }}_{\mathit{i}}\right)\mathit{\gamma }+\left({\mathit{\lambda }}_{\mathit{R}}-{\mathit{\lambda }}_{\mathit{E}}\right){\mathit{\epsilon }}_{\mathit{e}}+\mathit{\mu }{\mathit{\lambda }}_{\mathit{E}},\\ \frac{\mathit{d}{\mathit{\lambda }}_{\mathit{Q}}}{\mathit{d}\mathit{Q}}=-{\mathit{h}}_2+\left({\mathit{\lambda }}_{\mathit{Q}}-{\mathit{\lambda }}_{\mathit{S}}\right)\mathit{ŋ}\mathit{S}+\left({\mathit{\lambda }}_{\mathit{S}}-{\mathit{\lambda }}_{\mathit{Q}}\right){\mathit{\delta }}_1+\left({\mathit{\lambda }}_{\mathit{Q}}-{\mathit{\lambda }}_{\mathit{E}}\right)\left(1-{\mathit{\kappa }}_{\mathit{i}}\right)\mathit{\phi }\mathit{E}+\left({\mathit{\lambda }}_{\mathit{M}}-{\mathit{\lambda }}_{\mathit{Q}}\right)\left(1-{\mathit{\kappa }}_{\mathit{i}}\right){\mathit{\delta }}_2+\mathit{\mu }{\mathit{\lambda }}_{\mathit{Q}},\\ \frac{\mathit{d}{\mathit{\lambda }}_{\mathit{I}}}{\mathit{d}\mathit{I}}=-{\mathit{h}}_3+\left(1-{\mathit{\kappa }}_{\mathit{i}}\right)\left({\mathit{\lambda }}_{\mathit{I}}-{\mathit{\lambda }}_{\mathit{E}}\right)\mathit{\gamma }\mathit{E}+\left({\mathit{\lambda }}_{\mathit{M}}-{\mathit{\lambda }}_{\mathit{I}}\right)\mathit{\alpha }+\left({\mathit{\lambda }}_{\mathit{R}}-{\mathit{\lambda }}_{\mathit{I}}\right){\mathit{\epsilon }}_{\mathit{i}}+\mathit{\mu }{\mathit{\lambda }}_{\mathit{I}},\\ \frac{\mathit{d}{\mathit{\lambda }}_{\mathit{M}}}{\mathit{d}\mathit{M}}=-{\mathit{h}}_4+\left(1-{\mathit{\kappa }}_{\mathit{i}}\right)\left({\mathit{\lambda }}_{\mathit{M}}-{\mathit{\lambda }}_{\mathit{Q}}\right){\mathit{\delta }}_2\mathit{Q}+\left({\mathit{\lambda }}_{\mathit{M}}-{\mathit{\lambda }}_{\mathit{I}}\right)\mathit{\alpha }\mathit{I}+\left({\mathit{\lambda }}_{\mathit{R}}-{\mathit{\lambda }}_{\mathit{M}}\right){\mathit{\epsilon }}_{\mathit{m}}+\mathit{\mu }{\mathit{\lambda }}_{\mathit{M}},\\ \frac{\mathit{d}{\mathit{\lambda }}_{\mathit{R}}}{\mathit{d}\mathit{R}}=\left({\mathit{\lambda }}_{\mathit{R}}-{\mathit{\lambda }}_{\mathit{E}}\right){\mathit{\epsilon }}_{\mathit{e}}\mathit{E}+\left({\mathit{\lambda }}_{\mathit{R}}-{\mathit{\lambda }}_{\mathit{I}}\right){\mathit{\epsilon }}_{\mathit{i}}\mathit{I}+\left({\mathit{\lambda }}_{\mathit{R}}-{\mathit{\lambda }}_{\mathit{M}}\right){\mathit{\epsilon }}_{\mathit{m}}\mathit{M}+\left({\mathit{\lambda }}_{\mathit{S}}-{\mathit{\lambda }}_{\mathit{R}}\right)\mathit{\tau }+\mathit{\mu }{\mathit{\lambda }}_{\mathit{R}}.\end{array}\right.(32)$

The control set illustrates the co-state system with the optimal conditions.

$\left\{\begin{array}{l}\frac{\partial \mathit{H}}{\partial {\mathit{\kappa }}_{\mathit{p}}}={\mathit{z}}_1{\mathit{\kappa }}_{\mathit{p}}+\left({\mathit{\lambda }}_{\mathit{S}}-{\mathit{\lambda }}_{\mathit{E}}\right){\mathit{\rho }}_1^0{\mathit{S}}^0+\left({\mathit{\lambda }}_{\mathit{V}}-{\mathit{\lambda }}_{\mathit{E}}\right){\mathit{\rho }}_2^0{\mathit{V}}^0,\\ \frac{\partial \mathit{H}}{\partial {\mathit{\kappa }}_{\mathit{v}}}={\mathit{z}}_2{\mathit{\kappa }}_{\mathit{v}}+\left({\mathit{\lambda }}_{\mathit{V}}-{\mathit{\lambda }}_{\mathit{S}}\right)\mathit{\vartheta }{\mathit{S}}^0,\\ \frac{\partial \mathit{H}}{\partial {\mathit{\kappa }}_{\mathit{i}}}={\mathit{z}}_3{\mathit{\kappa }}_{\mathit{i}}+\left({\mathit{\lambda }}_{\mathit{E}}-{\mathit{\lambda }}_{\mathit{Q}}\right)\mathit{\phi }{\mathit{E}}^0+\left({\mathit{\lambda }}_{\mathit{E}}-{\mathit{\lambda }}_{\mathit{I}}\right)\mathit{\gamma }{\mathit{E}}^0+\left({\mathit{\lambda }}_{\mathit{Q}}-{\mathit{\lambda }}_{\mathit{M}}\right){\mathit{\delta }}_2{\mathit{Q}}^0.\end{array}\right.(33)$

We solve for ${\kappa }_p\left(t\right),{\kappa }_v\left(t\right)$, and ${\kappa }_i\left(t\right)\text{\hspace{0.17em}as\hspace{0.17em}}{\kappa }_p^{*}\left(t\right),{\kappa }_v^{*}\left(t\right)$, and ${\kappa }_i^{*}\left(t\right)$ of Eq. (33), and the results confirms the expression in Eq. (31) are as follows:

$\left\{\begin{array}{l}{\mathit{\kappa }}_{\mathit{p}}^{*}\left(\mathit{t}\right)=\min \left\{{\mathit{\kappa }}_{\mathit{p}}\max ,\max \left(0,\frac{\left({\mathit{\lambda }}_{\mathit{E}}-{\mathit{\lambda }}_{\mathit{S}}\right){\mathit{\rho }}_1^0{\mathit{S}}^0+\left({\mathit{\lambda }}_{\mathit{E}}-{\mathit{\lambda }}_{\mathit{V}}\right){\mathit{\rho }}_2^0{\mathit{V}}^0}{{\mathit{z}}_1}\right)\right\},\\ {\mathit{\kappa }}_{\mathit{v}}^{*}\left(\mathit{t}\right)=\min {\left\{{\mathit{\kappa }}_{\mathit{v}}\max ,\max \left(0,\frac{\left({\mathit{\lambda }}_{\mathit{S}}-{\mathit{\lambda }}_{\mathit{V}}\right)\mathit{\vartheta }{\mathit{S}}^0}{{\mathit{z}}_2}\right)\right\}}^{\prime }\\ {\mathit{\kappa }}_{\mathit{i}}^{*}\left(\mathit{t}\right)=\min \left\{{\mathit{\kappa }}_{\mathit{i}}\max ,\max \left(0,\frac{\left({\mathit{\lambda }}_{\mathit{Q}}-{\mathit{\lambda }}_{\mathit{E}}\right)\mathit{\phi }{\mathit{E}}^0+\left({\mathit{\lambda }}_{\mathit{I}}-{\mathit{\lambda }}_{\mathit{E}}\right)\mathit{\gamma }{\mathit{E}}^0+\left({\mathit{\lambda }}_{\mathit{M}}-{\mathit{\lambda }}_{\mathit{Q}}\right){\mathit{\delta }}_2{\mathit{Q}}^0}{{\mathit{z}}_3}\right)\right\}.\end{array}\right.(34)$

Therefore, using the bounds of the controls ${\kappa }_p\left(t\right),{\kappa }_v\left(t\right)$, and ${\kappa }_i\left(t\right),$ the control efforts are in the compact form given by the optimal condition of the system in Eq. (34); hence, the proof is complete.

6.2 Optimal control strategies

Here, our aim is to determine the number of infections after deploying the optimal control interventions. We explore the effects of implementing the interventions; therefore, the optimality system (36) is solved forward in time and the adjoint system backward in time with the corresponding lower and upper bounds of the controls. We used the population of Ghana to study the behavioral pattern of COVID-19. The estimated total population of Ghana is $31732129$ [19]; hence, $N\left(0\right)=31732129,$ and the assumed initial values are as follows: $S\left(0\right)=200000,V\left(0\right)=120000,E\left(0\right)=150000,I\left(0\right)=997,Q\left(0\right)=1000,M\left(0\right)=800,\text{and\hspace{0.17em}}R\left(0\right)=500,$ together with Λ = 1,364, $\tau =0.2$, and $ŋ=0.21$ and parameter values illustrated in Table 3. The balance costs associated with the objective functional are assumed to be $z_1=5,z_2=10,\text{and\hspace{0.17em}}{z}_3=20,$ and weight $h_i=100$, where $i=1,2,3,4.$ The lower bound (LB) and upper bound (UB) are assumed to be ${LB}_1=0,{UB}_1=1,{LB}_2=0,{UB}_2=1,{LB}_3=0,\text{and\hspace{0.17em}}{UB}_3=1.$ The results are illustrated according to the strategies to implement the control strategies.

6.2.1 Strategy 1: implementation of public awareness (${\mathit{\kappa }}_{\mathit{p}}$)

The optimal solutions illustrated in Figure 6 account for the observations when the control effort ${\kappa }_p$ is applied accordingly.

Figure 6

Figure 6. Optimal solutions of implementing strategy 1. Plot (A) blue dotted line is the optimal solution for implementing strategy 1. (B) Optimal control profile for strategy 1.

The optimal solutions illustrated above depict the following observations when public education is only applied:

(a) Figure 6A represents the effect of the control effort ${\kappa }_p$ on the infectious individuals. It implies that the number of individuals will decrease if the control effort is optimally implemented in halting the disease’s transmission. Conversely, it will increase significantly.

(b) Figure 6B represents the profile of the control effort for public awareness of COVID-19. It implies that education on COVID-19 should reach 50% of the population throughout the implementation to halt COVID-19 transmission.

6.2.2 Strategy 2: implementation of vaccination (${\kappa }_v$)

The optimal solutions illustrated in Figure 7 account for the observations when the control effort ${\kappa }_v$ is applied accordingly.

Figure 7

Figure 7. Optimal solutions of implementing strategy 2. Plot (A) blue dotted line is the optimal solution for implementing strategy 2. (B) Optimal control profile for strategy 2.

The optimal solutions illustrated above depict the following observations when vaccination is only applied:

(a) Figure 7A represents the effect of the control effort ${\kappa }_v$ on the infectious individuals. It implies that the number of individuals will decrease if the control effort is optimally implemented in halting the disease’s transmission. Conversely, it will increase significantly.

(b) Figure 7B represents the profile of the control effort for vaccination to prevent COVID-19. It implies that approximately 25% of the population should be vaccinated within 80 days and intensified further up to 75% in the subsequent days throughout the implementation to halt COVID-19 transmission.

6.2.3 Strategy 3: implementation of infectivity treatment (${\mathit{\kappa }}_{\mathit{i}}$)

The optimal solutions illustrated in Figure 8 account for the observations when the control effort ${\kappa }_i$ is applied accordingly.

Figure 8

Figure 8. Optimal solutions of implementing strategy 3. Plot (A) blue dotted line is the optimal solution for implementing strategy 3. (B) Optimal control profile for strategy 3.

The optimal solutions illustrated above depict the following observations when infectivity treatment is only applied:

(a) Figure 8A represents the effect of the control effort ${\kappa }_i$ on the infectious individuals. It implies that the number of individuals will decrease if the control effort is optimally implemented in halting the disease’s transmission. Conversely, it will increase significantly.

(b) Figure 8B represents the profile of the control effort for COVID-19 reinfection. It implies that approximately 75% of the individuals suspected to be carriers of the virus should be treated/monitored throughout the implementation to halt COVID-19 transmission.

6.2.4 Strategy 4: implementation of public awareness and vaccination (${\kappa }_p,{\kappa }_v$)

The optimal solutions illustrated in Figure 9 account for the observations when the control efforts ${\kappa }_p\text{\hspace{0.17em}and\hspace{0.17em}}{\kappa }_v$ are applied accordingly.

Figure 9

Figure 9. Optimal solutions of implementing strategy 4. Plot (A) blue dotted line is the optimal solution for implementing strategy 4. (B) Optimal control profile for strategy 4.

The optimal solutions illustrated above depict the following observations when control efforts for public awareness and vaccination are applied:

(a) Figure 9A represents the effect of the control efforts ${\kappa }_p\text{\hspace{0.17em}and\hspace{0.17em}}{\kappa }_v$ on the infectious individuals. It implies that the number of individuals will decrease if the control effort is optimally implemented in halting the disease’s transmission. Conversely, it will increase significantly.

(b) Figure 9B represents the profile of the control efforts for public awareness and vaccination against COVID-19. It implies that approximately 80% of the population should be educated on COVID-19 within 22 days, which can be relaxed to approximately 25%, and 50% of the population should be vaccinated throughout the implementation of these interventions to halt COVID-19 transmission.

6.2.5 Strategy 5: implementation of public awareness and infectivity treatment (${\kappa }_p,{\kappa }_i$)

The optimal solutions illustrated in Figure 10 account for the observations when all the control efforts ${\kappa }_p\text{\hspace{0.17em}and\hspace{0.17em}}{\kappa }_i$ are applied accordingly.

Figure 10

Figure 10. Optimal solutions of implementing strategy 5. Plot (A) blue dotted line is the optimal solution for implementing strategy 5. (B) Optimal control profiles for strategy 5.

The optimal solutions illustrated above depict the following observations when all the control efforts for public awareness and infectivity treatment are applied:

(a) Figure 10A represents the effect of the control efforts ${\kappa }_p\text{\hspace{1em}and\hspace{0.17em}}{\kappa }_i$ on the infectious individuals. It implies that the number of individuals will decrease to a minimum within 30 days if the control effort is optimally implemented in halting the disease’s transmission. Conversely, it will increase significantly.

(b) Figure 10B represents the profile of the control efforts for public awareness and infectivity treatment. It implies that approximately 80% of the population should be educated on COVID-19, and 50% of suspected carriers of the virus should be treated/monitored throughout the implementation of these interventions to halt COVID-19 transmission.

6.2.6 Strategy 6: implementation of vaccination and infectivity treatment (${\mathit{\kappa }}_{\mathit{v}},{\mathit{\kappa }}_{\mathit{i}}$)

The optimal solutions illustrated in Figure 11 account for the observations when the control efforts ${\kappa }_v\text{\hspace{0.17em}and\hspace{0.17em}}{\kappa }_i$ are applied accordingly.

Figure 11

Figure 11. Optimal solutions of implementing strategy 6. Plot (A) blue dotted line is the optimal solution for implementing strategy 6. (B) Optimal control profiles for strategy 6.

The optimal solutions illustrated above depict the following observations when all the control efforts for vaccination infectivity reduction are applied:

(a) Figure 11A represents the effect of the control efforts ${\kappa }_v\text{\hspace{0.17em}and\hspace{0.17em}}{\kappa }_i$ on the infectious individuals. It implies that the number of individuals will decrease to a minimum within 20 days if the control effort is optimally implemented in halting the disease’s transmission. Conversely, it will increase significantly.

(b) Figure 11B represents the profile of the control efforts for public education and vaccination against COVID-19. It implies that all the interventions $U_1\text{\hspace{0.17em}and\hspace{0.17em}}{U}_2$ should be implemented at levels higher than 25% from the start of implementation throughout the subsequent days to halt COVID-19 transmission.

6.2.7 Strategy 7: implementation of all controls (${\mathit{\kappa }}_{\mathit{p}},{\mathit{\kappa }}_{\mathit{v}},{\mathit{\kappa }}_{\mathit{i}}$)

The optimal solutions illustrated in Figure 12 account for the observations when the control efforts ${\kappa }_p,{\kappa }_v,\text{and\hspace{0.17em}}{\kappa }_i$ are applied accordingly.

Figure 12

Figure 12. Optimal solutions of implementing strategy 7. Plot (A) blue dotted line is the optimal solution for implementing strategy 7. (B) Optimal control profiles for strategy 7.

The optimal solutions illustrated above depict the following observations when all the control efforts are applied:

(a) Figure 12A represents the effect of the control efforts ${\kappa }_p,{\kappa }_v,\text{and\hspace{0.17em}}{\kappa }_i$ on the infectious individuals. It implies that the number of individuals will decrease to a minimum within 10 days if the control effort is optimally implemented in halting the disease’s transmission. Conversely, it will increase significantly.

(b) Figure 12B represents the profile of all control efforts. It implies that approximately 80% of the population should be educated on COVID-19, approximately 50% of the population should be vaccinated, and 50% of suspected carriers of the virus should also be treated/monitored throughout the implementation period of these interventions to halt COVID-19 transmission.

6.3 Cost–benefit analysis

Once the strategies are given, it is imperative to know the cost associated with implementing such intervention(s). Therefore, we explore the cost associated with each control strategy to check their effectiveness. We outline some cost-effectiveness approaches to further understand the control strategies.

We consider two procedures, namely, average cost-effectiveness ratio (ACER) and incremental cost-effectiveness ratio (ICER), which have been explained in [20–23], to carry out epidemiological studies.

6.4 Average cost-effectiveness ratio

We define the ACER of implementing a strategy as

$\text{ACER}=\frac{\text{Overall}\text{cost}\text{generated}\text{by}\text{applying}\text{the}\text{strategy}}{\text{Overall}\text{infection}\text{averted}\text{by}\text{applying}\text{the}\text{strategy}}.(35)$

The overall cost $Z$ stated in (27) would be used to evaluate the total cost that the intervention would generate in Eq. 36. We then compare the ACER values of each strategy, and the one with the lowest value is the most cost-effective, saving costs. Therefore, the cost-effective intervention is considered the strategy with the least ACER value. The expression in Eq. 35 is illustrated as follows.

From Table 4, control strategy 1, which involves the implementation of public education only, has the least value of ACER, indicating cost savings. However, relying solely on this metric is not enough to choose a strategy; hence, we further explore other approaches.

Table 4

Table 4. Strategies’ ACER values with their overall infection averted and cost incurred.

6.5 Incremental cost-effectiveness ratio

We define the ICER of implementing a strategy as

$\text{ICER}=\frac{\text{The}\text{cost}\text{difference}\text{generated}\text{by}\text{strategies}x\text{and}y}{\text{Difference}\text{in}\text{the}\text{overall}\text{infection}\text{averted}\text{in}\text{strategies}x\text{and}y}.(36)$

The total cost function $Z$ stated in (27) would be used to estimate the overall cost that the intervention would generate in Eq. (36). It is worth knowing that the averted total number of infections is the difference between the initial values of $E_x\text{\hspace{0.17em}and\hspace{0.17em}}{I}_x$, without control(s) and with controls. The outcomes are tabulated in increasing order of infection averted.

The ICER in Table 5 is calculated as

$\text{ICER}\left(2\right)=\frac{1.2513\times {10}^3-0}{2.0184\times {10}^7-0}=6.1995\times {10}^{-5},$

$\text{ICER}\left(1\right)=\frac{1.4064\times {10}^3-1.2513\times {10}^3}{1.5269\times {10}^8-2.0184\times {10}^7}=1.1705\times {10}^{-6},$

$\text{ICER}\left(4\right)=\frac{1.5279\times {10}^3-1.4064\times {10}^3}{1.5279\times {10}^8-1.5269\times {10}^8}=0.0012,$

$\text{ICER}\left(3\right)=\frac{5.6206\times {10}^3-1.5279\times {10}^3}{3.9351\times {10}^8-3.9351\times {10}^8}=1.7002\times {10}^{-5},$

$\text{ICER}\left(5\right)=\frac{6.2463\times {10}^3-5.6206\times {10}^3}{3.9388\times {10}^8-2.3231\times {10}^8}=0.0017,$

$\text{ICER}\left(6\right)=\frac{6.8719\times {10}^3-6.2463\times {10}^3}{5.2611\times {10}^8-3.9388\times {10}^8}=4.7312\times {10}^{-6},$

$\text{ICER}\left(7\right)=\frac{7.4975\times {10}^3-6.8719\times {10}^3}{5.2648\times {10}^8-5.2611\times {10}^8}=0.0017.$

Table 5

Table 5. Strategies’ ICER values with their overall infection averted and cost incurred.

Assessing strategies 2 and 1 in Table 5, it is noticed from the ICER that strategy 2 is expensive to deploy in a resource-limited setting; therefore, strategy 2 is removed from the list of possible controls, and the ICER is calculated again. This is presented in Table 6.

Table 6

Table 6. Strategies’ ICER values with their overall infection averted and cost incurred.

Assessing strategies 1 and 4 in Table 6, it is noticed from the ICER that strategy 4 is expensive to deploy in a resource-limited setting; therefore, strategy 4 is removed from the list of possible controls, and the ICER is calculated again. This is presented in Table 7.

Table 7

Table 7. Strategies’ ICER values with their overall infection averted and cost incurred.

Assessing strategies 1 and 3 in Table 7, it is noticed from the ICER that strategy 3 is expensive to deploy in a resource-limited setting; therefore, strategy 3 is removed from the list of possible controls, and the ICER is calculated again. This is presented in Table 8.

Table 8

Table 8. Strategies’ ICER values with their overall infection averted and cost incurred.

Assessing strategies 1 and 5 in Table 8, it is noticed from the ICER that strategy 5 is expensive to deploy in a resource-limited setting; therefore, strategy 5 is removed from the list of possible controls, and the ICER is calculated again. This is presented in Table 9.

Table 9

Table 9. Strategies’ ICER values with their overall infection averted and cost involved.

Assessing strategies 1 and 6 in Table 9, it is noticed from the ICER that strategy 6 is expensive to deploy in a resource-limited setting; therefore, strategy 6 is removed from the list of possible controls, and the ICER is calculated again. This is presented in Table 10.

Table 10

Table 10. Strategies’ ICER values with their overall infection averted and cost incurred.

Finally, assessing strategies 1 and 7 in Table 10, it is noticed from the ICER that strategy 7 is expensive to deploy in a resource-limited setting; therefore, strategy 7 is removed from the list of possible controls. Therefore, we conclude that strategy 1 is the most cost-effective strategy to use among the several strategies under study here. From the above analysis, it is obvious that strategy 1, which involves public education, is the intervention that saves cost.

7 Conclusion

We have presented a work that analyzes the changing effects of vaccination and infectivity reductions on the transmission of COVID-19 using data from Ghana. We have estimated the model’s parameters and analyzed their effects on disease transmission through numerical and graphical illustrations. Again, we have exhibited the threshold dynamics of the effective reproduction number $R_{\mathrm{e}}$ together with the contributions from the transmission routes (Table 2). We have demonstrated the sensitivity of the model’s parameters to study their effects on the effective reproduction number $R_{\mathrm{e}}$ (Table 3; Figure 2; Figure 3).

The aim of this work is to study the effect of vaccination and infectivity reductions in controlling COVID-19 transmission and devise control interventions that save cost to mitigate the transmission; therefore, we have formulated optimal control strategies together with the cost–benefit analysis that consider control measures involving both pharmaceutical and non-pharmaceutical interventions to control COVID-19. We implemented the strategies (Figure 6A–Figure 12A), and it was realized that public education, vaccination, and infectivity reductions to prevent COVID-19 should be intensified and reach approximately 25% of the population from the beginning and intensified in the subsequent days (Figure 6B–Figure 12B).

It is also worth knowing that public education saves cost as per the cost–benefit analysis compared to the other strategies raised in this work. This intervention can minimize COVID-19, as illustrated in Figure 6A. This intervention should reach approximately 50% of the population throughout the period of its implementation in order to realize the results of strategy 1 (Figure 6). Although strategy 1 saves cost, other strategies elaborated in this work can also be applied, but one has to consider the cost involved in implementing the strategy. The cost involved in applying the optimal control strategies is presented in Table 4–Table 10.

The outcomes of the findings imply that both pharmaceutical and non-pharmaceutical measures are very important in controlling the transmission of COVID-19. These control measures should always be vigorously implemented to create public awareness on COVID-19 and its reinfection, as illustrated in Figure 6B–Figure 12B, in order to reduce the effective contact rates and rates of acquiring COVID-19, as illustrated in Figure 3.

Although we have demonstrated the dynamics of COVID-19 transmission with vaccination and different infectivity reductions, this work is focused on the homogeneity of the population, and we hope to extend this study to explore the transmission dynamics of COVID-19 reinfection by considering heterogeneity of the population such as age and sex. We encourage individuals to adhere to personal hygiene and be aware of COVID-19 reinfection.

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 authors.

Author contributions

RA: conceptualization, data curation, formal analysis, investigation, methodology, writing–original draft, and writing–review and editing. ZJ: funding acquisition, supervision, and writing–review and editing. JY: formal analysis, supervision, and writing–review and editing. JA: formal analysis, investigation, and writing–review and editing.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This research was supported by the National Natural Science Foundation of China grants 12231012 and 61873154.

Acknowledgments

The authors appreciate the support of Complex Systems Research Center, Shanxi University.

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.

References

1. World Health Organization. Coronavirus disease 2019(COVID-19), situation report -51, data as reported by 8 March 2020 (2019). Available at: https://www.who.int/emergencies/disease/novel-coronavirus-2019/situation-reports (Accessed January 10, 2024).

Google Scholar

2. Chen Y, Wang Y, Fleming J, Yu Y, Gu Y, Liu C, et al. Active or latent tuberculosis increases susceptibility to COVID-19 and disease severity (2020). Available at: https://www.medrxiv.org/ (Accessed January 10, 2024).

Google Scholar

3. Salman AM, Ahmed I, Mohd MH, Jamiluddin MS, Dheyab MA. Scenario analysis of COVID-19 transmission dynamics in Malaysia with the possibility of reinfection and limited medical resources scenarios. Comput Biol Med (2021) 133:104372. doi:10.1016/j.compbiomed.2021.104372

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Zamir M, Nadeem F, Alqudah MA, Abdeljawad T. Future implications of COVID-19 through Mathematical modeling. Results Phys (2022) 33:105097. doi:10.1016/j.rinp.2021.105097

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Worldometers. Reported cases and deaths by country or territory (2022). Available at: https://www.worldometers.info/coronavirus/#countries (Accessed April 30, 2022).

Google Scholar

6. Zegarra MAA, Infante SD, Carrasco DB, Liceaga DO. COVID-19 optimal vaccination policies: a modeling study on efficacy, natural and vaccine-induced immunity responses. Math Biosci (2021) 337:108614. doi:10.1016/j.mbs.2021.108614

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Adachi A. Grand challenge in human/animal virology: unseen, smallest replicative entities shape thewholeglobe. FrontMicrobiol (2019) 11:431. doi:10.3389/fmicb.2020.00431

PubMed Abstract | CrossRef Full Text | Google Scholar

8. 1mg. Coronavirus No-panic help guide (2020). Available at: https://smef.org.uk/wpcontent/uploads/2020/03/Corona-Ebook.pdf. (Accessed January 10, 2024)

Google Scholar

9. Annas S, Pratama MI, Rifandi M, Sanusi W, Side S, Stability analysis and numerical simulation of SEIR model for pandemic COVID-19 spread in Indonesia. Chaos, Solitons and Fractal (2020) 139:110072. doi:10.1016/j.chaos.2020.110072

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Pang L, Liu S, Zhang X, Tian T, Zhao Z, Transmission dynamics and control strategies of COVID-19 in Wuhan, China. J Biol Syst (2020) 28:543–60. doi:10.1142/s0218339020500096

CrossRef Full Text | Google Scholar

11. Asamoah JKK, Owusu MA, Jin Z, Oduro F, Abidemi A, Gyasi EO, Global stability and cost-effectiveness analysis of COVID-19 considering the impact of the environment: using data from Ghana. Chaos, Solitons and Fractals (2020) 140:110103. doi:10.1016/j.chaos.2020.110103

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Tchoumi S, Diagne M, Rwezaura H, Tchuenche J, Malaria and COVID-19 co-dynamics: a mathematical model and optimal control. Appl Math Model (2021) 99:294–327. e327. doi:10.1016/j.apm.2021.06.016

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Yan X, Zou Y, Optimal and sub-optimal quarantine and isolation control in SARS epidemics. Math Comp Model (2008) 47(1-2):235–45. doi:10.1016/j.mcm.2007.04.003

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Khan MA, Atangana A, Mathematical modeling and analysis of COVID-19: a study of new variant omicron. Physica A: Stat Mech Its Appl (2022) 599:127452. Article 127452. doi:10.1016/j.physa.2022.127452

CrossRef Full Text | Google Scholar

15. Rwezaura H, Diagne M, Omame A, de Espindola A, Tchuenche J, Mathematical modeling and optimal control of SARS-CoV-2 and tuberculosis co-infection: a case study of Indonesia. Model Earth Syst Environ (2022) 8(4):5493–520. e5520. doi:10.1007/s40808-022-01430-6

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Omame A, Abbas M, Onyenegecha C, A fractional-order model for COVID-19 and tuberculosis co-infection using Atangana–Baleanu derivative. Chaos, Solitons and Fractals (2021) 153:111486. Article 111486. doi:10.1016/j.chaos.2021.111486

PubMed Abstract | CrossRef Full Text | Google Scholar

17. World Health Organization. World health organization covid– 19 report (2024). Available at: https://covid19.who.int/region/afro/country/gh. (Accessed January 10, 2024).

Google Scholar

18. Dietz K. The estimation of the basic reproduction number for infectious diseases. Stat Methods Med Res (1993) 2:23–41. doi:10.1177/096228029300200103

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Osei E, Amu H, Kye-Duodu G, Kwabla MP, Danso E, Binka FN, et al. Impact of COVID-19 pandemic on Tuberculosis and HIV services in Ghana: an interrupted time series analysis. PLoS ONE (2023) 18(9):e0291808. doi:10.1371/journal.pone.0291808

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Asamoah JKK, Owusu MA, Jin Z, Oduro FT, Abidemi A, Gyasi EO, Global stability and cost-effectiveness analysis of COVID-19 considering the impact of the environment: using data from Ghana. Chaos, Solitons and Fractals (2020) 140:110103. doi:10.1016/j.chaos.2020.110103

PubMed Abstract | CrossRef Full Text | Google Scholar

21. Agusto F, Leite M, Optimal control and cost-effective analysis of the 2017 meningitis outbreak in Nigeria. Infect Dis Model (2017) 4(2019):161–87. doi:10.1016/j.idm.2019.05.003

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Asamoah JKK, Okyere E, Abidemi A, Moore SE, Sun G, 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

23. Asamoah JKK, Jin Z, Sun G, Non-seasonal and seasonal relapse model for Q fever disease with comprehensive cost-effectiveness analysis. Results Phys (2021) 22:103889. doi:10.1016/j.rinp.2021.103889

CrossRef Full Text | Google Scholar