Predictive Modeling and Control Strategies for the Transmission of Middle East Respiratory Syndrome Coronavirus

Abstract

The Middle East respiratory syndrome coronavirus (MERS-CoV) is a highly infectious respiratory illness that poses a significant threat to public health. Understanding the transmission dynamics of MERS-CoV is crucial for effective control and prevention strategies. In this study, we develop a precise mathematical model to capture the transmission dynamics of MERS-CoV. We incorporate some novel parameters related to birth and mortality rates, which are essential factors influencing the spread of the virus. We obtain epidemiological data from reliable sources to estimate the model parameters. We compute its basic reproduction number (R0). Stability theory is employed to analyze the local and global properties of the model, providing insights into the system’s equilibrium states and their stability. Sensitivity analysis is conducted to identify the most critical parameter affecting the transmission dynamics. Our findings revealed important insights into the transmission dynamics of MERS-CoV. The stability analysis demonstrated the existence of stable equilibrium points, indicating the long-term behavior of the epidemic. Through the evaluation of optimal control strategies, we identify effective intervention measures to mitigate the spread of MERS-CoV. Our simulations demonstrate the impact of time-dependent control variables, such as supportive care and treatment, in reducing the number of infected individuals and controlling the epidemic. The model can serve as a valuable tool for public health authorities in designing effective control and prevention strategies, ultimately reducing the burden of MERS-CoV on global health.

1. Introduction

The first identification of the Middle East respiratory syndrome coronavirus (MERS-CoV), a viral respiratory illness, took place in Saudi Arabia in 2012, as reported by multiple studies, including those conducted by [1,2,3]. MERS-CoV is believed to have originated from an animal source and has been identified in both humans and animals. The transmission of the disease occurs through close contact with an infected individual, in any form. The World Health Organization (2019) has reported a global total of 2519 laboratory-confirmed cases of MERS-CoV infection, with 866 associated deaths. One of the largest outbreaks of MERS-CoV occurred in South Korea in 2015.

Coronaviruses constitute a diverse family of viruses that are known to infect humans, causing respiratory illnesses that can vary in severity from mild cold-like symptoms to severe respiratory syndromes, such as severe acute respiratory syndrome (SARS). MERS-CoV can cause zoonotic infections in humans through direct or indirect contact with camels or camel-related products. Additionally, human-to-human transmission has been reported, particularly in healthcare settings [4,5,6]. Since 2002, three novel coronaviruses have emerged and caused deadly zoonotic diseases in humans. The first was SARS, which emerged in November 2002. The second was MERS, which emerged in April 2012. The most recent and ongoing pandemic is COVID-19, which emerged in December 2019 and has affected millions of people worldwide [7,8,9,10]. The authors in [11,12] modified a mathematical model of COVID-19 by including the quarantine class and measured the disease transmission. Shen et al. studied a mathematical model of COVID-19 by presenting the vaccinated class with an optimal control analysis [13]. The authors in [14] used the data of Saudi Arabia and investigated the transmission dynamics of COVID-19. Tsay et al. analyzed the state estimation and optimal control for the COVID-19 outbreak model in the US [15]. Libotte et al. used the optimal strategy for vaccines in COVID-19 treatment considered for both mono- and multi-objective optimization [16]. Using the optimal control models also gives information about the impact of individual vaccination during an epidemic, together with the key considerations for political and economic decision making.

MERS typically presents with symptoms such as fever, cough, and shortness of breath. It is believed to spread through respiratory secretions, such as through coughing, from an infected person, similar to other coronaviruses. Several studies have investigated the potential role of camel handlers in the transmission of the virus to determine its source of infection. To understand and predict the dynamics of infectious diseases, researchers have developed various models based on biologically feasible parameters [17,18,19]. These models are essential tools for analyzing and forecasting the spread of diseases [20]. Although several case studies have explored the transmission of MERS-CoV, the literature on its transmission dynamics is limited. Using available data, Cauchemez et al. [21] estimated the incubation period and generation time of MERS-CoV, and calculated the reproductive numbers for both animal-to-human and human-to-human transmission. Chowell et al. [22] took a different approach and compared the reproductive numbers of SARS and MERS. Assiri et al. [23] reported one of the largest outbreaks of MERS-CoV, describing the virus as transmissible from human to human. The virus has spread globally through travel-associated cases, with reported incidences in countries including Algeria, Austria, China, Egypt, Italy, Netherlands, Philippines, South Korea, Thailand, the UK, and the US. Groundbreaking research has been published by numerous esteemed researchers, delving into the exploration of various model types, such as: SIR epidemic models [24,25], the discrete-time prey–predator model [26], and the memristor system [27]. Several infectious disease models have been investigated by researchers by using different approaches, which are available in the literature, such as [28,29,30,31,32].

Members of the coronavirus family, MERS-CoV and COVID-19 (caused by SARS-CoV-2) have certain similarities. It is crucial to remember that these viruses are diverse from one another and have distinctive traits and effects on human health. The following examples demonstrate how MERS-CoV can be used to treat various illnesses, including COVID-19. Both COVID-19 and MERS-CoV can cause serious respiratory infections in people. However, their overall effects differ considerably in a number of ways. MERS-CoV and COVID-19 can also be compared to influenza viruses, particularly those that cause severe respiratory infections. Different influenza viruses (A, B, and C) are what cause the illness. Despite the fact that some symptoms and the means of transmission are similar, influenza viruses have unique genetic traits and often create seasonal outbreaks. However, MERS-CoV and SARS-CoV-2 can cause pandemics or sporadic epidemics with ongoing human-to-human transmission. Several researchers have suggested mathematical models by applying different approaches to an infectious disease and studying its dynamics from different angles [33,34,35,36,37].

As far as the novelty is concerned, we study the model presented in [38], by incorporating the natural birth rate and death rate due to MERS-CoV. We modified the mathematical model for MERS-CoV transmission dynamics. This model consists of six groups: susceptible class $\mathcal{S}$, exposed class (or high risk latent) $\mathcal{E}$, symptomatic and infectious class $\mathcal{I}$, infectious but asymptotic class $\mathcal{A}$, hospitalized class $\mathcal{H}$, and recovery class $\mathcal{R}$. After constructing the model, the basic reproductive number is calculated by using the next generation method, and the local and global stability of the equilibrium points are determined. Lyapunov function theory is then utilized to analyze the global behavior of the model. Furthermore, the principles of optimal control theory are employed to reduce the number of infected persons and maximize the recovery rate within a given population.

2. Model Formulation

Here, we study the mathematical formulation of the deterministic model for MERS-CoV, using a set of differential equations. Specifically, the model describes the dynamics of the host population using the following system of equations:

$$\begin{array}{cc}& \dot{\mathcal{S}}(t)=bN-\frac{\phi \mathcal{I}\mathcal{S}}{\mathrm{N}}-\frac{\phi q\mathcal{H}\mathcal{S}}{\mathrm{N}}-{\eta }_0\mathcal{S},\hfill \\ & \dot{\mathcal{E}}(t)=\frac{\phi \mathcal{I}\mathcal{S}}{\mathrm{N}}+\frac{\phi q\mathcal{H}\mathcal{S}}{\mathrm{N}}-(\chi +{\eta }_0)\mathcal{E},\hfill \\ & \dot{\mathcal{I}}(t)=\chi \xi \mathcal{E}-({\vartheta }_a+{\vartheta }_1)\mathcal{I}-({\eta }_0+{\eta }_1)\mathcal{I},\hfill \\ & \dot{\mathcal{A}}(t)=\chi (1-\xi )\mathcal{E}-({\eta }_0+{\eta }_2)\mathcal{A},\hfill \\ & \dot{\mathcal{H}}(t)={\vartheta }_a\mathcal{I}-{\vartheta }_{\vartheta }\mathcal{H}-{\eta }_0\mathcal{H},\hfill \\ & \dot{\mathcal{R}}(t)={\vartheta }_1\mathcal{I}+{\vartheta }_{\vartheta }\mathcal{H}-{\eta }_0\mathcal{R},\hfill \end{array}$$

(1)

with initial conditions

$$\mathcal{S}(0)\ge 0,\mathcal{E}(0)\ge 0,\mathcal{I}(0)\ge 0,\mathcal{A}(0)\ge 0,\mathcal{H}(0)\ge ,\mathcal{R}\ge 0,$$

where the used parameters in the above system are: $b\mathrm{N}$ represents the rate of birth for the host populace, while the transmission rate from human to human per unit time is represented by $\phi$. The parameter q determines the relative transmissibility of hospitalized individuals. $\chi$ represents the rate at which individuals transition from the exposed compartment $\mathcal{E}$ to the infectious compartment $\mathcal{I}$. The proportion of individuals who progress from $\mathcal{E}$ to $\mathcal{I}$ is given by $\xi$, while the remaining ($1-\xi$) progress to class $\mathcal{A}$. The average rate at which symptomatic persons are hospitalized is denoted by ${\vartheta }_a$, while ${\vartheta }_1$ represents the rate of recovery without hospitalization, and ${\vartheta }_{\vartheta }$ represents the rate of recovery of hospitalized patients. The rate of natural death is represented by ${\eta }_0$, while ${\eta }_1$ and ${\eta }_2$ represent deaths due to MERS-CoV.

Assume that the total populace is represented by $\mathrm{N}(t)$ at time t, and satisfies $\mathrm{N}(t)=\mathcal{S}+\mathcal{E}+\mathcal{I}+\mathcal{A}+\mathcal{H}+\mathcal{R}.$

Adding all the equations of system (1), we have

$$\frac{d\mathrm{N}}{dt}=b\mathrm{N}-{\eta }_0\mathcal{S}-{\eta }_0\mathcal{E}-({\eta }_0+{\eta }_1)\mathcal{I}-({\eta }_0+{\eta }_2)\mathcal{A}-{\eta }_0\mathcal{H}-{\eta }_0\mathcal{R}.$$

Therefore, from the above relation for biological applications, the considered system (1) occurred in the closed set as

$$\mathbb{F}=\left\{(\mathcal{S},\mathcal{E},\mathcal{I},\mathcal{A},\mathcal{H},\mathcal{R})\in {\mathcal{R}}_{+}^6,0<,\mathcal{S}+\mathcal{E}+\mathcal{I}+\mathcal{A}+\mathcal{H}+\mathcal{R}\le \frac{b\mathrm{N}}{{\eta }_0}\right\}.$$

3. Basic Reproduction Number

The basic reproduction number determines whether an epidemic will appear or the infection will die out. It represents the expected average number of new infections that will be generated by a single infective person, both directly and indirectly, when introduced into a fully susceptible populace. In this study, we use the approach of Driessche and Watmough [39,40] to calculate the basic reproduction number for the aforementioned system (1).

$$\mathrm{F}=\left[\begin{array}{ccc}0& \frac{\phi {\mathcal{S}}_0}{\mathrm{N}}& \frac{\phi q{\mathcal{S}}_0}{\mathrm{N}}\\ \\ 0& 0& 0\\ \\ 0& 0& 0\end{array}\right], \mathrm{V}=\left[\begin{array}{ccc}\chi +{\eta }_0& 0& 0\\ \\ -\chi \xi & ({\vartheta }_a+{\vartheta }_1)+({\eta }_0+{\eta }_1)& 0\\ \\ 0& -{\vartheta }_a& {\vartheta }_{\vartheta }+{\eta }_0\end{array}\right]$$

(2)

to find

$$\mathrm{F}{\mathrm{V}}^{-1}=\left(\begin{array}{ccc}\frac{\chi \phi \xi {\mathcal{S}}_0({\vartheta }_{\vartheta }+{\eta }_0+q{\vartheta }_a)\mathcal{I}}{\mathrm{N}(\chi +{\eta }_0)({\vartheta }_{\vartheta }+{\eta }_0)[({\vartheta }_a+{\vartheta }_1)+({\eta }_0+{\eta }_1)]}& \frac{\phi {\mathcal{S}}_0({\vartheta }_{\vartheta }+{\eta }_0+q{\vartheta }_a)}{\mathrm{N}({\vartheta }_{\vartheta }+{\eta }_0)[({\vartheta }_a+{\vartheta }_1)+({\eta }_0+{\eta }_1))]}& \frac{\phi {\mathcal{S}}_{0}q}{\mathrm{N}[({\vartheta }_{\vartheta }+{\eta }_0)({\vartheta }_a+{\vartheta }_1)+({\eta }_0+{\eta }_1)]}\\ \\ 0& 0& 0\\ \\ 0& 0& 0\end{array}\right).$$

Thus, the required basic reproduction number ${\mathbf{R}}_0$ is followed by

$${\mathbf{R}}_0=\frac{\chi \phi \xi {\mathcal{S}}_0{Q}_1}{\mathrm{N}(\chi +{\eta }_0)Q_2{Q}_3}.$$

The terms $Q_1$, $Q_2$, and $Q_3$ are defined as follows:

$Q_1=({\vartheta }_{\vartheta }+{\eta }_0+q{\vartheta }_a)$;

$Q_2=({\vartheta }_{\vartheta }+{\eta }_0)$;

$Q_3=({\vartheta }_a+{\vartheta }_1)+({\eta }_0+{\eta }_1)$.

These terms correspond to the susceptible individuals at the disease-free equilibrium (DFE).

Analysis of Sensitivity

Here, we conduct a sensitivity analysis of some of the parameters used in the model. This technique helps us to identify the parameters that have a significant effect on the basic reproduction number (See

Table 1. Sensitivity indices of different parameters.

Notation	Sensitivity Values	Notation	Sensitivity Values
$\chi$	0.0384651	${\vartheta }_1$	−0.058565
${\vartheta }_{\vartheta }$	−0.00000453	${\vartheta }_a$	−0.93704
${\eta }_0$	−0.001464128	${\eta }_1$	−0.041392
$\phi$	0.99999	q	0.0000476
$bN$	0.99999	$\xi$	0.000432

and Figure 1). We use the approach as described by Chintis [41] to calculate the sensitivity index of ${\mathbf{R}}_0$. Specifically, the sensitivity index ${\Delta }_h^{{\mathbf{R}}_0}$ of a parameter h is presented by the formula ${\Delta }_h^{{\mathbf{R}}_0}=\frac{\partial {\mathbf{R}}_0}{\partial k}\frac{h}{{\mathbf{R}}_0}$.

4. Equilibria Points

The aforementioned system (1) has two possible equilibria: one is the disease-free equilibrium (DFE) and second one is the endemic equilibrium (EE). The DFE, denoted by ${\mathcal{F}}_0$, is given by ${\mathcal{F}}_0=\left(\frac{b\mathrm{N}}{{\eta }_0},0,0,0,0,0\right)$. The EE, denoted by ${\mathcal{F}}_1$, is found by setting “$\mathcal{S}={\mathcal{S}}_{*}$, $\mathcal{E}={\mathcal{E}}_{*}$, $\mathcal{I}={\mathcal{I}}_{*}$, $\mathcal{A}={\mathcal{A}}_{*}$, $\mathcal{H}={\mathcal{H}}_{*}$, and $\mathcal{R}={\mathcal{R}}_{*}$, and the LHS of the resulting system to zero”. We obtain the following expression after simplification, ${\mathcal{S}}_{*}$, ${\mathcal{E}}_{*}$, ${\mathcal{I}}_{*}$, ${\mathcal{A}}_{*}$, ${\mathcal{H}}_{*}$, and ${\mathcal{R}}_{*}$ at the EE.

$$\begin{array}{ccc}\hfill & & {\mathcal{S}}_{*}=\frac{\mathrm{N}b\mathrm{N}{Q}_2}{\phi Q_2+\phi q{\vartheta }_a+{\eta }_0\mathrm{N}{Q}_2},\hfill \\ \\ & & {\mathcal{E}}_{*}=\frac{\mathrm{N}b\mathrm{N}{Q}_2+\phi q{\vartheta }_a+\mathrm{N}b\mathrm{N}{Q}_2}{\phi Q_2+\phi q{\vartheta }_a+{\eta }_0\mathrm{N}{Q}_2},\hfill \\ \\ & & {\mathcal{I}}_{*}=\frac{(\chi \xi )(\mathrm{N}b\mathrm{N}{Q}_2+({\mathbf{R}}_0-1)+\mathrm{N}b\mathrm{N}{Q}_2)}{Q_3\phi +\phi q{\vartheta }_a+{\eta }_0\mathrm{N}{Q}_2},\hfill \\ \\ & & {\mathcal{A}}_{*}=\frac{\chi (1-\xi )\mathrm{N}b\mathrm{N}{Q}_2^2+\phi q{\vartheta }_a+\mathrm{N}b\mathrm{N}{Q}_2}{\phi Q_2+\phi q{\vartheta }_a+{\eta }_0\mathrm{N}{Q}_3},\hfill \\ \\ & & {\mathcal{H}}_{*}=\frac{{\vartheta }_a{\mathcal{I}}_{*}}{Q_2},\hfill \\ \\ & & {\mathcal{R}}_{*}=\frac{Q_2{\mathcal{I}}_{*}+Q_1}{Q_2}.\hfill \end{array}$$

(3)

4.1. Local Stability

We show the local asymptotic stability (LAS) of the DFE as well as the EE of the system (1) with the help of the following theorem.

Theorem 1. 

If the basic reproductive number ${\mathbf{R}}_0$ is less than 1, the DFE point is LAS.

Proof. 

To show the local stability of the system, about the point DFE, the Jacobian matrix for the said system (1) is

$$\begin{array}{c}\hfill J_0=\left(\begin{array}{cccccc}-{\eta }_0& 0& -\frac{\phi {\mathcal{S}}_0}{\mathrm{N}}& 0& \frac{\phi q{\mathcal{S}}_0}{\mathrm{N}}& 0\\ 0& -(\chi +{\eta }_0)& \frac{\phi {\mathcal{S}}_0}{\mathrm{N}}& 0& \frac{\phi q{\mathcal{S}}_0}{\mathrm{N}}& 0\\ 0& \chi \xi & -Q_3& 0& 0& 0\\ 0& \chi (1-\xi )& 0& -{\eta }_2& 0& 0\\ 0& 0& {\vartheta }_a& 0& -Q_2& 0\\ 0& 0& {\vartheta }_1& 0& {\vartheta }_{\vartheta }& -{\eta }_0\end{array}\right).\end{array}$$

(4)

By conducting a row operation, reducing the matrix to echelon form, the following Jacobian matrix is obtained

$$\begin{array}{c}\hfill \left(\begin{array}{cccccc}-{\eta }_0& 0& -\frac{\phi {\mathcal{S}}_0}{\mathrm{N}}& 0& \frac{\phi q{\mathcal{S}}_0}{\mathrm{N}}& 0\\ 0& -(\chi +{\eta }_0)& \frac{\phi {\mathcal{S}}_0}{\mathrm{N}}& 0& \frac{\phi q{\mathcal{S}}_0}{\mathrm{N}}& 0\\ 0& 0& \mathbf{A}& \frac{\phi q{\mathcal{S}}_0\chi \xi }{\mathrm{N}}& 0& \frac{\phi q{\mathcal{S}}_0\chi \xi }{\mathrm{N}}\\ 0& 0& 0& \mathbf{B}& \frac{\phi q{\mathcal{S}}_0\chi }{\mathrm{N}}& 0\\ 0& 0& 0& 0& \mathbf{C}& 0\\ 0& 0& 0& 0& 0& \mathbf{D}\end{array}\right),\end{array}$$

(5)

$$\begin{array}{ccc}\hfill \mathbf{A}& =& -Q_3{Q}_2(\chi +{\eta }_0)-\frac{\phi {\mathcal{S}}_0\chi \xi }{\mathrm{N}},\hfill \\ \hfill \mathbf{B}& =& -({\eta }_0+{\eta }_2)(\chi +{\eta }_0)Q_1-\frac{\phi {\mathcal{S}}_0\chi \xi }{\mathrm{N}},\hfill \\ \hfill \mathbf{C}& =& -Q_1{Q}_2{Q}_3({\mathcal{S}}_0+{\eta }_0),\hfill \\ \hfill \mathbf{D}& =& -{\eta }_0{Q}_3(\kappa +{\eta }_0)-\frac{(\phi {\mathcal{S}}_0\chi \xi )}{\mathrm{N}}-[(1-{\mathbf{R}}_0)(N(\chi +{\eta }_0)Q_2{Q}_3)\phi q{\mathcal{S}}_0\chi \xi ].\hfill \end{array}$$

According to [42], when ${\mathbf{R}}_0<1$, the matrices $\mathbf{A}$, $\mathbf{B}$, $\mathbf{C}$, and $\mathbf{D}$ are all negative, and the eigenvalues have negative real parts. As a result, the DFE is LAS. □

Theorem 2. 

If ${\mathbf{R}}_0$ is greater than 1, the EE point is LAS.

Proof. 

Consider the Jacobian of the considered problem (1) at ${\mathcal{F}}_1$ is,

$$\begin{array}{c}\hfill J_0=\left(\begin{array}{cccccc}-\frac{\phi {\mathcal{I}}_{*}}{\mathrm{N}}-\frac{\phi {\mathcal{H}}_{*}q}{\mathrm{N}}-{\eta }_0& 0& -\frac{\phi {\mathcal{S}}_{*}}{\mathrm{N}}& 0& -\frac{\phi q{\mathcal{S}}_{*}}{\mathrm{N}}& 0\\ \frac{\phi {\mathcal{I}}_{*}}{\mathrm{N}}+\frac{\phi {\mathcal{H}}_{*}q}{\mathrm{N}}& -(\chi +{\eta }_0)& \frac{\phi {\mathcal{S}}_{*}}{\mathrm{N}}& 0& \frac{\phi q{\mathcal{S}}_{*}}{\mathrm{N}}& 0\\ 0& \chi \xi & Q_3& 0& 0& 0\\ 0& \chi (1-\xi )& 0& -({\eta }_0+{\eta }_1)& 0& 0\\ 0& 0& {\vartheta }_a& 0& -({\vartheta }_{\vartheta }+{\eta }_0)& 0\\ 0& 0& {\vartheta }_1& 0& {\vartheta }_{\vartheta }& -{\eta }_0\end{array}\right).\end{array}$$

(6)

After performing a row operation and simplifying the resulting expressions, we obtain the following Jacobian matrix:

$$\begin{array}{c}\hfill \left(\begin{array}{cccccc}-\frac{\phi {\mathcal{I}}_{*}}{\mathrm{N}}-\frac{\phi {\mathcal{H}}_{*}q}{\mathrm{N}}-{\eta }_0& 0& -\frac{\phi {\mathcal{S}}_{*}}{\mathrm{N}}& 0& -\frac{\phi q{\mathcal{S}}_{*}}{\mathrm{N}}& 0\\ 0& -(\chi +{\eta }_0)(\frac{\phi {\mathcal{I}}_{*}}{\mathrm{N}}+\frac{\phi {\mathcal{H}}_{*}q{\mathcal{S}}_{*}}{\mathrm{N}}+{\eta }_0)& -\frac{\phi {\mathcal{S}}_{*}}{\mathrm{N}}& 0& -\frac{\phi q{\mathcal{S}}_{*}}{\mathrm{N}}& 0\\ 0& 0& {\mathbb{Z}}_1& {\eta }_2\xi & 0& 0\\ 0& 0& 0& {\mathbb{Z}}_2& {\mathbb{Z}}_3& {\mathbb{Z}}_4\\ 0& 0& 0& 0& {\mathbb{Z}}_5& {\eta }_0{\vartheta }_{\vartheta }\\ 0& 0& 0& 0& 0& {\mathbb{Z}}_6\end{array}\right),\end{array}$$

(7)

where

$$\begin{array}{ccc}\hfill & & {\mathbb{Z}}_1=-\xi Q_3,\hfill \\ \\ \hfill & & {\mathbb{Z}}_2=-{\eta }_2(\kappa +{\eta }_0)(\phi {\mathcal{I}}_{*}+\phi q{\mathcal{H}}_{*}+{\eta }_0){\vartheta }_1,\hfill \\ \\ \hfill & & {\mathbb{Z}}_3=-{\vartheta }_{\vartheta }\phi {\mathcal{S}}_{*}{\eta }_0\chi \xi ,\hfill \\ \\ \hfill & & {\mathbb{Z}}_4=-{\eta }_0\phi {\mathcal{S}}_{*}\chi (1-\xi ),\hfill \\ \\ \hfill & & {\mathbb{Z}}_5=-Q_2{\vartheta }_1-{\vartheta }_a-{\vartheta }_{\vartheta },\hfill \\ \\ \hfill & & {\mathbb{Z}}_6=-{\eta }_0{Q}_3({\eta }_2(\chi +{\eta }_0)(\phi {\mathcal{I}}_{*}\hfill \\ \\ \hfill & & +\phi q{\mathcal{H}}_{*}+{\eta }_0){\vartheta }_1({\mathbf{R}}_0-1)\times [(\chi +{\eta }_0)Q_3](\phi {\mathcal{S}}_{*}{\eta }_0\phi ({\eta }_2\xi {\vartheta }_I)({\vartheta }_{\vartheta }+{\eta }_0){\vartheta }_1+{\vartheta }_a\hfill \\ \\ \hfill & & +Q_3({\eta }_0\xi Q_3({\eta }_2(\chi +{\eta }_0))(\phi {\mathcal{I}}_{*}+\phi q{\mathcal{H}}_{*}+{\eta }_0){\vartheta }_1\hfill \\ \\ \hfill & & +({\vartheta }_{\vartheta }\phi {\mathcal{S}}_{*}{\eta }_0\chi \xi ){\eta }_2\xi {\vartheta }_1.\hfill \\ \end{array}.$$

The eigenvalues are given by

$$\begin{array}{ccc}\hfill & & {\zeta }_1=-\frac{\phi {\mathcal{I}}_{*}}{\mathrm{N}}-\frac{\phi {\mathcal{H}}_{*}q}{\mathrm{N}}-{\eta }_0<0,\hfill \\ \\ \hfill & & {\zeta }_2=-(\chi +{\eta }_0)(\frac{\phi {\mathcal{I}}_{*}}{\mathrm{N}}+\frac{\phi {\mathcal{H}}_{*}q{\mathcal{S}}_{*}}{\mathrm{N}}+{\eta }_0)<0,\hfill \\ \\ \hfill & & {\zeta }_3=-\xi ({\vartheta }_a+{\vartheta }_1+{\eta }_0+{\eta }_1)={\mathbb{Z}}_1<0,\hfill \\ \\ \hfill & & {\zeta }_4=-{\eta }_2(\kappa +{\eta }_0)(\phi {\mathcal{I}}_{*}+\phi q{\mathcal{H}}_{*}+{\eta }_0){\vartheta }_1={\mathbb{Z}}_2,<0,\hfill \\ \\ \hfill & & {\zeta }_5=-({\vartheta }_{\vartheta }+{\eta }_0){\vartheta }_1-{\vartheta }_a-{\vartheta }_{\vartheta }={\mathbb{Z}}_5<0,\hfill \\ \\ \hfill & & {\zeta }_6={\mathbb{Z}}_6<0.\hfill \end{array}.$$

(8)

As per the findings reported in [43], when ${\mathbf{R}}_0>1$, all of the eigenvalues have nonpositive real parts, which indicates that the EE point is LAS. □

4.2. Analysis of Global Stability

The next theorem presents that the said system is globally asymptotically stable (GAS) for the DFE and EE point.

Theorem 3. 

The DFE of the system is GAS for ${\mathbf{R}}_0<1$, otherwise unstable.

Proof. 

We define the Lyapunov function as follows:

$$U(t)=k_1(\mathcal{S}-{\mathcal{S}}_0)+k_2\mathcal{E}+k_3\mathcal{I}+k_4\mathcal{A}+k_5\mathcal{H}.$$

(9)

We differentiate Equation (9) and obtain:

$$\begin{array}{ccc}\hfill U^{\prime }(t)& =& k_1{\mathcal{S}}^{\prime }+k_2{\mathrm{E}}^{\prime }+k_3{\mathcal{I}}^{\prime }+k_4{\mathcal{A}}^{\prime }+k_5{\mathcal{H}}^{\prime }.\hfill \end{array}$$

(10)

Using model (1), we obtain

$$\begin{array}{ccc}\hfill U^{\prime }(t)& =& k_1[b\mathrm{N}-\frac{\phi \mathcal{I}\mathcal{S}}{\mathrm{N}}-\frac{\phi q\mathcal{H}\mathcal{S}}{\mathrm{N}}-{\eta }_{0}S]+k_2[\frac{\phi \mathcal{I}\mathcal{S}}{\mathrm{N}}+\frac{\phi q\mathcal{H}\mathcal{S}}{\mathrm{N}}-(\chi +{\eta }_0)\mathcal{E}]\hfill \\ & +& k_3[\chi \xi \mathcal{E}-({\vartheta }_a+{\vartheta }_1)\mathcal{I}-({\eta }_0+{\eta }_1)\mathcal{I}]+k_4[\chi (1-\xi )\mathcal{E}-({\eta }_0+{\eta }_2)\mathcal{A}]\hfill \\ & +& k_5[{\vartheta }_a\mathcal{I}-{\vartheta }_{\vartheta }\mathcal{H}-{\eta }_0\mathcal{H}].\hfill \end{array}$$

If we choose the positive parameter values $k_1=k_2=k_4=\xi$, $k_3=1$, and $k_5=q\phi \rho$, and simplify, we obtain:

$$\begin{array}{ccc}\hfill U^{\prime }(t)& =& -\xi {\eta }_0(\mathcal{S}-{\mathcal{S}}_0)-2\frac{\phi q\mathcal{H}\mathcal{S}}{\mathrm{N}}-Q_1\xi {\eta }_0\mathcal{E}-({\eta }_0+{\eta }_2)\mathcal{A}\hfill \\ & -& Q_3[1-{\mathbf{R}}_0]-Q_2{Q}_3(\kappa +{\eta }_0)-Q_2^2\mathcal{H}.\hfill \end{array}$$

where

$${\mathcal{S}}_0=\frac{b\mathrm{N}}{{\eta }_0},$$

Let $U^{\prime }(t)$ be a function of time t, and let S and ${\mathbf{R}}_0$ be constants. If $S>{\mathcal{S}}_0$ and ${\mathbf{R}}_0<1$, then $U^{\prime }(t)$ is negative. If $\mathcal{S}={\mathcal{S}}_0$, then $U^{\prime }(t)=0$. According to the LaSalle invariance principle [44,45], if $\mathcal{E}=\mathcal{I}=\mathcal{A}=\mathcal{H}=0$, then the set of initial conditions for which $U^{\prime }(t)$ approaches zero as t approaches infinity is an invariant set.

Therefore, the DFE ${\mathcal{F}}_0$ is GAS. □

Theorem 4. 

When ${\mathbf{R}}_0>1$, the EE point is GAS at ${\mathcal{F}}_1$, and unstable when ${\mathbf{R}}_0<1$.

Proof. 

For the GAS of the EE point, we define the Lyapunov function as:

$$U(t)=\frac{1}{2}{[p_1(\mathcal{S}-{\mathcal{S}}_{*})+p_2(\mathcal{E}-{\mathcal{E}}_{*})+p_3(\mathcal{I}-{\mathcal{I}}_{*})+p_4(\mathcal{A}-{\mathcal{A}}_{*})+p_5(\mathcal{H}-{\mathcal{H}}_{*})]}^2,$$

(11)

and we introduce the constants $p_1,p_2,p_3,p_4$, and $p_5$, which will be chosen later. Upon differentiating Equation (11), we obtain:

$$\begin{array}{ccc}\hfill U^{\prime }(t)& =& [p_1(\mathcal{S}-{\mathcal{S}}_{*})+p_2(\mathcal{E}-{\mathcal{E}}_{*})+p_3(\mathcal{I}-{\mathcal{I}}_{*})+p_4(\mathcal{A}-{\mathcal{A}}_{*})\hfill \\ & +& p_5(\mathcal{H}-{\mathcal{H}}_{*})][p_1(\frac{ds}{dt})+p_2(\frac{dE}{dt})+p_3(\frac{dI}{dt})+p_4(\frac{dA}{dt})+p_5(\frac{dH}{dt})],\hfill \end{array}$$

$$\begin{array}{ccc}\hfill U^{\prime }(t)& =& p_1(\mathcal{S}-{\mathcal{S}}_{*})+p_2(\mathcal{E}-{\mathcal{E}}_{*})+p_3(\mathcal{I}-{\mathcal{I}}_{*})+p_4(\mathcal{A}-{\mathcal{A}}_{*})+p_5(\mathcal{H}-{\mathcal{H}}_{*})(p_1(b\mathrm{N}-\frac{\phi \mathcal{I}\mathcal{S}}{\mathrm{N}}-\frac{\phi q\mathcal{H}\mathcal{S}}{\mathrm{N}}-{\eta }_0\mathcal{S})\hfill \\ & +& p_2(\frac{\phi \mathcal{I}\mathcal{S}}{\mathrm{N}}+\frac{\phi q\mathcal{H}\mathcal{S}}{\mathrm{N}}-(\chi +{\eta }_0)\mathcal{E}+p_3(\chi \xi \mathcal{E}-({\vartheta }_a+{\vartheta }_1)\mathcal{I}-({\eta }_0+{\eta }_1)\mathcal{I})+p_4(\chi (1-\xi )\mathcal{E}-({\eta }_0+{\eta }_2)\mathcal{A}\hfill \\ & +& p_5({\vartheta }_a\mathcal{I}-{\vartheta }_{\vartheta }\mathcal{H}-{\eta }_0\mathcal{H}).\hfill \end{array}$$

After some calculation, we obtain, and utilizing the values of $p_1,\dots p_5$, we obtain

$$U^{\prime }(t)=-\frac{(\phi {\mathcal{S}}_{*})}{\mathrm{N}}(\mathcal{S}-{\mathcal{S}}_{*})-({\mathbf{R}}_0-1)Q_1{Q}_2-\frac{\phi q{\mathcal{H}}_{*}{\mathcal{S}}_{*}}{\mathrm{N}}[\mathcal{E}+\mathcal{I}+\mathcal{A}]-({\vartheta }_a\mathcal{I}-{\vartheta }_{\vartheta })Q_2{Q}_3^2{\mathcal{H}}_{*}$$

For $\mathcal{S}={\mathcal{S}}_{*}$ and $({\vartheta }_a>{\vartheta }_{\vartheta })$ for ${\mathbf{R}}_0$ is greater than 1; thus, the proof is finished. □

5. Results and Discussion

In this context, we substantiate our analytical discoveries through the application of the fourth-order Runge–Kutta method [46]. We select certain parameters for illustrative purposes, while obtaining others from published data sources [38]. The parameters employed in the simulation are chosen with careful consideration of their biological plausibility. The ensuing set of parameters is utilized for the subsequent analysis. $\phi =0.007; q=0.003; \chi =0.005;$ $\xi =0.0001; {\eta }_0=0.0003; {\eta }_1=0.0001; \xi =0.002; {\vartheta }_1=0.001; {\vartheta }_a=0.000001; {\vartheta }_{\vartheta }=0.0007; \mathrm{and}$ $bN=0.00004.$ To validate the analytical findings of the proposed model concerning the DFE, we employed the aforementioned parameter values. Subsequently, we computed the DFE point’s coordinates and the threshold parameter ${\mathbf{R}}_0$ as $(7.98337196,0,0,0,0)$ and $(0.043732)$, respectively. The simulation outcomes utilizing the aforementioned parameters are depicted in Figure 2 and Figure 3, thereby substantiating the analytical conclusions outlined in the theorem. To corroborate this, we employed the linear stability analysis technique and introduced perturbations to the initial compartmental population values. Remarkably, these perturbed values consistently converged to the DFE, underscoring its robustness against varying initial conditions, $\mathcal{S}(0)=1000$, $\mathcal{E}(0)=800$, $\mathcal{I}(0)=600$, $\mathcal{A}(0)=500$, $\mathcal{H}(0)=400$, and $\mathcal{R}(0)=300$. Drawing from the theoretical interpretation of the data, a definitive conclusion can be drawn: when the value of ${\mathbf{R}}_0$ is below 1, the disease transmission will inevitably diminish over time. This is evidenced by the convergence of every solution curve to a stable position, as depicted in the corresponding plots.

$$\begin{array}{ccc}\hfill \frac{{\mathcal{S}}^{i+1}-{\mathcal{S}}^i}{l}& =& bN-\frac{\phi {\mathcal{I}}^i{\mathcal{S}}^{i+1}}{\mathrm{N}}-\frac{\phi q{\mathcal{H}}^i{\mathcal{S}}^{i+1}}{\mathrm{N}}-{\eta }_0{\mathcal{S}}^{i+1},\hfill \\ \hfill \frac{{\mathcal{E}}^{i+1}-{\mathcal{E}}^i}{l}& =& \frac{\phi {\mathcal{I}}^i{\mathcal{S}}^{i+1}}{\mathrm{N}}+\frac{\phi q{\mathcal{H}}^i{\mathcal{S}}^{i+1}}{\mathrm{N}}-(\chi +{\eta }_0){\mathcal{E}}^{i+1},\hfill \\ \hfill \frac{{\mathcal{I}}^{i+1}-{\mathcal{I}}^i}{l}& =& \chi \xi {\mathcal{E}}^{i+1}-({\vartheta }_a+{\vartheta }_1){\mathcal{I}}^{i+1}-({\eta }_0+{\eta }_1){\mathcal{I}}^{i+1},\hfill \\ \hfill \frac{{\mathcal{A}}^{i+1}-{\mathcal{A}}^i}{l}& =& \chi (1-\xi ){\mathcal{E}}^{i+1}-({\eta }_0+{\eta }_2){\mathcal{A}}^{i+1},\hfill \\ \hfill \frac{{\mathcal{H}}^{i+1}-{\mathcal{H}}^i}{l}& =& {\vartheta }_a{\mathcal{I}}^{i+1}-{\vartheta }_{\vartheta }{\mathcal{H}}^{i+1}-{\eta }_0{\mathcal{H}}^{i+1},\hfill \\ \hfill \frac{{\mathcal{R}}^{i+1}-{\mathcal{R}}^i}{l}& =& {\vartheta }_1{\mathcal{I}}^{i+1}+{\vartheta }_{\vartheta }{\mathcal{H}}^{i+1}-{\eta }_0{\mathcal{R}}^{i+1}.\hfill \end{array}$$

Used Algorithm

Step 1: $({\mathcal{S}}_0,{\mathcal{E}}_0,{\mathcal{I}}_0,{\mathcal{A}}_0,{\mathcal{H}}_0,{\mathcal{R}}_0=0).$

Step 2: Let $i=1,2\dots n-1.$

$$\begin{array}{ccc}\hfill {\mathcal{S}}^{i+1}& =& \frac{\mathrm{N}lb}{\phi {\mathcal{I}}^i{\mathcal{S}}^{i+1}+\phi q{\mathcal{H}}^i{\mathcal{S}}^{i+1}l+{\eta }_{0}l\mathrm{N}}+\frac{{\mathcal{S}}^{i+1}}{\phi {\mathcal{I}}^i{\mathcal{S}}^{i+1}+\phi q{\mathcal{H}}^i{\mathcal{S}}^{i+1}l+{\eta }_{0}l\mathrm{N}},\hfill \\ \hfill {\mathcal{E}}^{i+1}& =& \frac{l\phi {\mathcal{I}}^i{\mathcal{S}}^{i+1}}{N(1+l(\chi +{\eta }_0))}+\frac{l\phi q{\mathcal{H}}^i{\mathcal{S}}^{i+1}}{N(1+l(\chi +{\eta }_0))}+\frac{{\mathcal{E}}^{i+1}}{(1+l(\chi +{\eta }_0))},\hfill \\ \hfill {\mathcal{I}}^{i+1}& =& \frac{l\chi \xi {\mathcal{E}}^{i+1}}{(1+l({\vartheta }_a+{\vartheta }_1)+({\eta }_0+{\eta }_1)l)}+\frac{{\mathcal{I}}^{i+1}}{1+l({\vartheta }_a+{\vartheta }_1)+({\eta }_0+{\eta }_1)l},\hfill \\ \hfill {\mathcal{A}}^{i+1}& =& \frac{l\chi (1-\xi ){\mathcal{E}}^{i+1}}{1+({\eta }_0+{\eta }_2)l}+\frac{{\mathcal{A}}^{i+1}}{1+({\eta }_0+{\eta }_2)l},\hfill \\ \hfill {\mathcal{H}}^{i+1}& =& \frac{l{\vartheta }_a{\mathcal{I}}^{i+1}}{1+{\eta }_{0}l+{\vartheta }_{\vartheta }l}+\frac{{\mathcal{H}}^{i+1}}{1+{\vartheta }_{\vartheta }+{\eta }_{0}l},\hfill \\ \hfill {\mathcal{R}}^{i+1}& =& \frac{l{\vartheta }_1{\mathcal{I}}^{i+1}}{1+{\eta }_{0}l}+\frac{l{\vartheta }_{\vartheta }{\mathcal{H}}^{i+1}}{1+{\eta }_{0}l}+\frac{{\mathcal{R}}^{i+1}}{1+{\eta }_{0}l}.\hfill \end{array}$$

Step 3: Let $i=1,2,3,\dots ,n-1$, by letting “${\mathcal{S}}_{*}(t_i)={\mathcal{S}}_{*},$ ${\mathcal{E}}_{*}(t_i)={\mathcal{E}}_{*},$ ${\mathcal{I}}_{*}(t_i)={\mathcal{I}}_{*},$ ${\mathcal{A}}_{*}(t_i)={\mathcal{A}}_{*},$ ${\mathcal{H}}_{*}(t_i)={\mathcal{H}}_{*}$, ${\mathcal{R}}_{*}(t_i)={\mathcal{R}}_{*}$.”

Subsequently, we proceed to explore the system’s dynamics around the EE by assuming an alternate set of parameters: $\phi =0.17$, $q=0.03$, $\chi =0.05$, $\xi =0.01$, ${\eta }_0=0.03$, ${\eta }_1=0.031$, $\xi =0.052$, ${\vartheta }_1=0.041$, ${\vartheta }_a=0.000001$, ${\vartheta }_{\vartheta }=0.0007$, and $bN=0.004$. Using these parameter values acquired earlier, we calculate the endemic equilibrium points and the associated ${\mathbf{R}}_0$ for the model (1). When ${\mathbf{R}}_0>1$, the endemic equilibrium point is determined to be $(40.76549,110.908700,70.45321,$ $85.934214,85.7659321,120.7659321)$, with a calculated value of ${\mathbf{R}}_0=7.13587$. Assuming the same initial population sizes for the compartments as in the previous analysis, the graphical results indicate that the populations of susceptible, exposed, infected asymptomatic, hospitalized, and recovered individuals initially undergo fluctuations before eventually stabilizing at their respective equilibrium values. For the parameter values employed in this study, the equilibrium point is $(24.76549,99.908700,22.45321,90.934214,85.7659321,120.7659321)$.

6. Analysis of Optimal Control

Here, we aim to establish an effective control strategy to prevent the spread of MERS in the population. Optimal control theory is a powerful mathematical technique that can be applied to design control schemes for a variety of infectious diseases. To achieve this, we apply optimal control theory, as described in previous works [46,47,48,49], to establish an appropriate control strategy. Our objective in this study is to reduce the prevalence of MERS in the populace by increasing the number of persons who recover from the disease, denoted as $\mathcal{R}$, and decreasing the number of individuals who are infectious, denoted as $\mathcal{I}$, and hospitalized, denoted as $\mathcal{H}$, by implementing time-dependent control variables such as treatment $v_1(t)$ and care $v_2(t)$. In model $(1)$, we take six state variables $\mathcal{S}$, $\mathcal{E}$, $\mathcal{I}$, $\mathcal{A}$, $\mathcal{H}$, and $\mathcal{R}$. Now, for the control problem, we take the two control variables, that is treatment $v_1(t)$ and care $v_2(t)$. Hence, we have the successive optimal control problem to reduce the objective functional

$$J(v_1,v_2)={\int }_0^T[c_1\mathcal{I}(t)+c_2\mathcal{H}(t)+\frac{1}{2}(c_3{v}_1^2(t)+c_4{v}_2^2(t)]dt$$

(12)

subject to

$$\begin{array}{cc}& \dot{\mathcal{S}}(t)=b\mathrm{N}-\frac{\phi \mathcal{I}\mathcal{S}}{\mathrm{N}}-\frac{\phi q\mathcal{H}\mathcal{S}}{\mathrm{N}}-{\eta }_0\mathcal{S},\hfill \\ & \dot{\mathcal{E}}(t)=\frac{\phi \mathcal{I}\mathcal{S}}{\mathrm{N}}+\frac{\phi q\mathcal{H}\mathcal{S}}{\mathrm{N}}-(\chi +{\eta }_0)\mathcal{E},\hfill \\ & \dot{\mathcal{I}}(t)=\chi \xi \mathcal{E}-({\vartheta }_a+{\vartheta }_1)\mathcal{I}-({\eta }_0+{\eta }_1)-v_1\mathcal{I},\hfill \\ & \dot{\mathcal{A}}(t)=\chi (1-\xi )\mathcal{E}-({\eta }_0+{\eta }_2)\mathcal{A},\hfill \\ & \dot{\mathcal{H}}(t)={\vartheta }_a\mathcal{I}-{\vartheta }_{\vartheta }\mathcal{H}-{\eta }_0\mathcal{H}-v_2\mathcal{H},\hfill \\ & \dot{\mathcal{R}}(t)={\vartheta }_1\mathcal{I}+{\vartheta }_{\vartheta }\mathcal{H}-{\eta }_0\mathcal{R}+v_1\mathcal{I}+v_2\mathcal{H},\hfill \end{array}$$

(13)

with initial conditions

$$“\mathcal{S}(0)\ge 0,\mathcal{E}(0)\ge 0,\mathcal{I}(0)\ge 0,\mathcal{A}(0)\ge 0,\mathcal{H}(0)\ge 0,\mathcal{R}(0)\ge 0”.$$

Equation (12) includes weight constants $c_1$, $c_2$, $c_3$, and $c_4$ that correspond to the relative importance of infected people I and hospitalized individual H in the objective function. The parameters $\frac{1}{2}{c}_3{v}_1^2$ and $\frac{1}{2}{c}_4{v}_2^2$ represent the costs associated with self-care and treatment. The primary objective is to evaluate the control function to achieve a specific goal.

$$J(v_1^{*},v_2^{*})=min\{J(v_1,v_2),v_1,v_2\in U\}$$

(14)

dependent on control system (13), where U in Equation (14) is known as the control set and is presented as,

$$“U=\{(v_1,v_2)/v_i(t)\mathrm{is}\mathrm{Lebesgue}\mathrm{measurable}\mathrm{on}[0,1],0\le v_i(t)\le 1,i=1,2\}.”$$

(15)

Before proceeding, it is important to establish the existence of control variables. According to Kamien and Aldila’s study [47], a solution for a state system can be found when the controls are bounded and Lebesgue measurable, in addition to satisfying the initial conditions. Thus, we can suppose that the considered control model can be formulated in the manner presented below.

$$\frac{d\varphi }{dt}=\mathcal{A}\varphi +\mathcal{B}\varphi .$$

From the above equation $\varphi =(\mathcal{S},\mathcal{E},\mathcal{I},\mathcal{A},\mathcal{H},\mathcal{R})$, where $\mathcal{A}(\varphi )$ and $\mathcal{B}(\varphi )$ denote the linear and nonlinear bounded coefficient ∋

$$\begin{array}{c}\hfill J_0=\left(\begin{array}{cccccc}-{\eta }_0& 0& 0& 0& 0& 0\\ 0& -(\chi +{\eta }_0)& 0& 0& 0& 0\\ 0& \chi \xi & -Q_3& 0& 0& 0\\ 0& \chi (1-\xi )& 0& -{\eta }_2& 0& 0\\ 0& 0& {\vartheta }_a& 0& -Q_2& 0\\ 0& 0& {\vartheta }_1& 0& {\vartheta }_{\vartheta }& -{\eta }_0\end{array}\right).\end{array}$$

(16)

$$\begin{array}{c}\hfill B(\varphi )=\left(\begin{array}{c}b\mathrm{N}-\frac{\phi \mathcal{I}\mathcal{S}}{\mathrm{N}}-\frac{\phi q\mathcal{H}\mathcal{S}}{\mathrm{N}}\\ \frac{\phi \mathcal{I}\mathcal{S}}{\mathrm{N}}+\frac{\phi q\mathcal{H}\mathcal{S}}{\mathrm{N}}\\ 0\\ 0\\ 0\\ 0\end{array}\right).\end{array}$$

(17)

Letting $\mathrm{L}(\varphi )=\mathcal{A}\varphi +\mathrm{F}\varphi ,$

$$\begin{array}{ccc}\hfill |\mathrm{F}({\varphi }_1)-\mathrm{F}({\varphi }_2)|& \le & m_1|{\mathcal{S}}_1-{\mathcal{S}}_2|+m_2|{\mathcal{E}}_1-{\mathcal{E}}_2|+m_3|{\mathcal{I}}_1-{\mathcal{I}}_2|+m_4|{\mathcal{A}}_1-{\mathcal{A}}_2|+m_5|{\mathcal{H}}_1-{\mathcal{H}}_2|+m_6|{\mathcal{R}}_1-{\mathcal{R}}_2|\hfill \\ & \le & \mathrm{N}|{\mathcal{S}}_1-{\mathcal{S}}_2|+|{\mathcal{E}}_1-{\mathcal{E}}_2|+|{\mathcal{I}}_1-{\mathcal{I}}_2|+|{\mathcal{A}}_1-{\mathcal{A}}_2|+|{\mathcal{H}}_1-{\mathcal{H}}_2|+|{\mathcal{R}}_1-{\mathcal{R}}_2|.\hfill \end{array}$$

Here, $\mathrm{N}=max(m_1,m_2,m_3,m_4,m_5,m_6)$ is a constant that is independent of the state variables in the aforementioned system. We also express

$$|\mathrm{L}({\varphi }_1)-\mathrm{L}({\varphi }_2)|\le \mathrm{M}|({\varphi }_1)-({\varphi }_2)|.$$

The solution for (13) exists due to the nonnegativity of the model state variables $\mathcal{S}$, $\mathcal{E}$, $\mathcal{I}$, $\mathcal{A}$, $\mathcal{H}$, and $\mathcal{R}$. Furthermore, it has been shown that the function L is Lipschitz uniformly continuous, and where $\mathrm{M}=(\mathrm{N},\parallel \mathcal{K}\parallel )<\infty$. Based on the properties mentioned earlier, we present the following theorem to establish the existence of a solution for model (1), which we then proceed to prove.

Theorem 5. 

For the control problem in Equations (12) and (13) there exists an optimal control as $v^{*}=(v_1^{*},v_2^{*})\in U$.

Proof. 

It is evident that the control and state variables in system (1) are positive. Additionally, the control variables set $\mathrm{U}$ is a closed and convex set, as mentioned in the problem statement. Furthermore, the control system is bounded, implying the compactness of the system. The integral in the objective function of the optimization problem, given by $c_1\mathcal{I}+c_2\mathcal{H}+\frac{1}{2}(c_3{v}_1^2(t)+c_4{v}_2^2(t))$, is also convex w.r.t the control set $\mathrm{U}$. This convexity guarantees the existence results for optimal control for the optimal control variables $(v_1^{*},v_2^{*})$. □

6.1. Methods

Next to show the optimal solution to the control model (12) and (13), we can apply the Lagrangian and Hamiltonian methods, described in the equation below

$$\begin{array}{ccc}\hfill \mathrm{L}(\mathcal{I},\mathcal{H},v_1,v_2)& =& c_1\mathcal{I}+c_2\mathcal{H}+\frac{1}{2}(c_3{v}_1^2(t)+c_4{v}_2^2(t).\hfill \end{array}$$

To describe the Hamiltonian ($\mathcal{H}$), by utilizing the notation $\vartheta$ = $({\vartheta }_1,{\vartheta }_2,{\vartheta }_3,{\vartheta }_4,{\vartheta }_5,{\vartheta }_6)$ and $y=(y_1,y_2,y_3,y_4,y_5,y_6)$, thus

$$\mathcal{H}(x,v,\vartheta )=\mathrm{L}(x,v)+\vartheta \mathrm{Z}(x,v),$$

where

$$\begin{array}{ccc}\hfill & & {\mathrm{Z}}_1=b\mathrm{N}-\frac{\phi \mathcal{I}\mathcal{S}}{\mathrm{N}}-\frac{\phi q\mathcal{H}\mathcal{S}}{\mathrm{N}}-{\eta }_0\mathcal{S},\hfill \\ \\ \hfill & & {\mathrm{Z}}_2=\frac{\phi \mathcal{I}\mathcal{S}}{\mathrm{N}}+\frac{\phi q\mathcal{H}\mathcal{S}}{\mathrm{N}}-(\chi +{\eta }_0)\mathcal{E},\hfill \\ \\ \hfill & & {\mathrm{Z}}_3=\chi \xi \mathcal{E}-({\vartheta }_a+{\vartheta }_1)\mathcal{I}-({\eta }_0+{\eta }_1)-v_1(t)\mathcal{I}(t),\hfill \\ \\ \hfill & & {\mathrm{Z}}_4=\chi (1-\xi )\mathcal{E}-({\eta }_0+{\eta }_2)\mathcal{A},\hfill \\ \\ \hfill & & {\mathrm{Z}}_5={\vartheta }_a\mathcal{I}-{\vartheta }_{\vartheta }\mathcal{H}-{\eta }_0\mathcal{H}-v_2(t)\mathcal{H}(t),\hfill \\ \\ \hfill & & {\mathrm{Z}}_6={\vartheta }_1\mathcal{I}+{\vartheta }_{\vartheta }\mathcal{H}-{\eta }_0\mathcal{R}+v_1(t)\mathcal{I}(t)+v_2(t)\mathcal{H}(t),\hfill \end{array}$$

(18)

Thus we apply the Pontryagin Maximum Principle [50,51] to the Hamiltonian in order to determine the optimal solution. According to this principle, if $(x^{*},v^{*})$ is an optimal solution, then there must exist a function $\vartheta$∋:

$$\begin{array}{c}\hfill \frac{dx}{dt}=\frac{\partial \mathcal{H}}{\partial \vartheta },0=\frac{\partial \mathcal{H}}{\partial u},\\ \hfill \vartheta {(t)}^{{}^{\prime }}=-\frac{\partial \mathcal{H}}{\partial x}.\end{array}$$

$$\mathcal{H}(t,x^{*},v^{*},\vartheta )\partial x=ma{x}_{v_1,v_2,v_3,v_4\in [0,1]}\mathcal{H}(x^{*}(t),v_1,v_2\vartheta (t));$$

(19)

with

$$\vartheta (t_f)=0,$$

(20)

The principles outlined in Equation (19) are utilized to determine the adjoint system (adjoint variables) and optimal control variables. Based on these principles, the following result can be obtained.

Theorem 6. 

Suppose ${\mathcal{S}}_{*}$, ${\mathcal{E}}_{*}$, ${\mathcal{I}}_{*}$, ${\mathcal{A}}_{*}$, ${\mathcal{H}}_{*}$, and ${\mathcal{R}}_{*}$ represent the optimal state solutions for the system, obtained using the combined optimal control variables $(v_1^{*},v_2^{*})$ that were derived through the numerical solution of the optimality system. The optimal control problem is defined by the objective function (12) and the control system (13). Then ∃ adjoint variables ${\vartheta }_1(t)$, ${\vartheta }_2(t)$, ${\vartheta }_3(t)$, and ${\vartheta }_4(t)$, ${\vartheta }_5(t)$, ${\vartheta }_6(t)$ satisfy

$$\begin{array}{ccc}\hfill {\vartheta }_1^{{}^{\prime }}(t)& =& -{\mathcal{A}}_1+({\vartheta }_2-{\vartheta }_1)\phi {\mathcal{I}}_{*}+({\vartheta }_2-{\vartheta }_1)\phi q{\mathcal{H}}_{*}-{\eta }_0{\vartheta }_1,\hfill \\ \hfill {\vartheta }_2^{{}^{\prime }}(t)& =& -{\mathcal{A}}_2+({\vartheta }_4-{\vartheta }_2)\phi {\mathrm{N}}^{*}+(\chi +{\eta }_0){\vartheta }_2-{\vartheta }_1{v}_1^{*}-{\vartheta }_3\xi ,\hfill \\ \hfill {\vartheta }_3^{{}^{\prime }}(t)& =& -{\mathcal{A}}_3+({\vartheta }_2-{\vartheta }_1)\phi {\mathcal{S}}_{*}+({\vartheta }_5-{\vartheta }_3){\vartheta }_a+({\vartheta }_6-{\vartheta }_3){\vartheta }_1+({\vartheta }_6-{\vartheta }_3)v_1(t)-v_1{\vartheta }_3,\hfill \\ \hfill {\vartheta }_4^{{}^{\prime }}(t)& =& -{\mathcal{A}}_4+(v_2-u_0){\vartheta }_4,\hfill \\ \hfill {\vartheta }_5{}^{\prime }(t)& =& -{\mathcal{A}}_5-({\vartheta }_2-{\vartheta }_1)\phi q{\mathcal{S}}_{*}+({\vartheta }_6-{\vartheta }_5){\vartheta }_{\vartheta }-(u_0+v_2){\vartheta }_5,\hfill \\ \hfill {\vartheta }_6{}^{\prime }(t)& =& -{\mathcal{A}}_6+u_0{\vartheta }_6,\hfill \end{array}$$

(21)

with boundary conditions.

Additionally, the optimal control parameters $v_1^{}(t)$ and $v_2^{}(t)$ are obtained through numerical solutions of the optimality problem and are presented below.

$$\begin{array}{ccc}\hfill v_1^{*}(t)& =& max\{min\{\frac{({\vartheta }_6-{\vartheta }_3){\mathcal{I}}_{*}}{B_1},1\},0\},\hfill \end{array}$$

(22)

$$\begin{array}{ccc}\hfill v_2^{*}(t)& =& max\{min\{\frac{({\vartheta }_6-{\vartheta }_5){\mathcal{H}}_{*}}{B_2},1\},0\}.\hfill \end{array}$$

(23)

Proof. 

The adjoint problem described by Equation (21) is obtained through the uses of the Pontryagin Maximum Principle given by Equation (19), while the transversal conditions arise from $\vartheta (T)=0$. The set of optimal functions $v_1^{*},v_2^{{}^{\prime }}$ is obtained using $\frac{\partial H}{\partial u}$. In the following section, we present numerical solutions to the optimality system in order to provide a clearer understanding for the reader, as opposed to relying solely on analytical results. The optimality problem is expressed by several components, including the control problem (13), the adjoint model (21), the boundary (terminal) conditions, and the optimal control functions. By solving these components numerically, we can gain valuable insights into the behavior of the system and assess its performance. □

6.2. Results and Discussion for Optimal Control

We utilize the Runge–Kutta method of order four to solve the optimal control system (13), in order to investigate the effects of self-care and treatment. To find the solution of the state system (12) with initial conditions in the time interval [0, 50], we employ the forward Runge–Kutta procedure. Similarly, the backward Runge–Kutta technique is used to solve the adjoint system (21) in the same interval with the assistance of the transversality condition. Below are the parameters that we used for the simulation: $bN=0.0071;$ $\phi =0.00041;$ $q=0.0000123;$ $\chi =0.0000123;$ $\xi =0.0000123;$ ${\vartheta }_1=0.003907997;$ ${\vartheta }_a=0.98;$ ${\vartheta }_{\vartheta }=0.0000404720925;$ $q=0.017816;$ $\rho =0.00007;and$ ${\eta }_0=0.00997.$ The weight constants $c_1,c_2,c_3,$ and $c_4$ were chosen based on biological feasibility. Specifically, we set $c_1=0.6610000,$ $c_2=0.54450,$ $c_3=0.0090030,$ $c_4=0.44440.$ The results obtained from the simulations are presented in Figure 4 and Figure 5.

Figure 4 and Figure 5 show the variations in the number of all compartments with and without control measures implemented.

7. Conclusions

The objective of this study is to developed a more realistic mathematical model that captures the transmission dynamics of the MERS-CoV. This is accomplished by introducing new parameters for the birth and death rates in the host populace. The threshold number ${\mathbf{R}}_0$ is a measure used to estimate the potential spread of a disease within a populace, and it can be calculated from a model to quantify the transmissibility of MERS-CoV. The model is analyzed using stability theory to identify conditions for local and global stability, and the most sensitive parameter is determined through a sensitivity analysis of ${\mathbf{R}}_0$. An optimal control problem is formulated with the goal of minimizing the number of infected persons and maximizing the number of recoveries in the population. The effectiveness of the approach is verified through numerical simulations, which demonstrate the stability of the results.