Modeling the epidemic trend of middle eastern respiratory syndrome coronavirus with optimal control

1.   Introduction

The coronavirus family comprises a diverse range of viruses that can be found in different animal species, including cats, camels, bats, and cattle. In rare cases, animal coronaviruses can infect humans and spread among them, as has been seen with SARS, MERS, and the current COVID-19 pandemic ^[1,2]. Since it first emerged in Saudi Arabia in 2012, the Middle East respiratory syndrome (MERS-CoV) has claimed the lives of 1791 individuals, according to various sources ^[3,4]. Meanwhile, the 2003 outbreak of severe acute respiratory syndrome (SARS) resulted in the deaths of 774 individuals ^[5].

Camels have been identified as the primary carriers of the MERS-CoV virus according to scientific research. Human-to-human transmission is the leading cause of MERS-CoV cases, responsible for 75 to 88 percent of all cases, while the remaining cases are caused by transmission from camels to humans. It is important to note that the virus can spread through respiratory discharge from infected individuals, such as coughing. In addition, close contact, including caring for or living with an infected person, can also result in the transmission of the virus. Since the discovery of MERS-CoV in April 2012, there have been a total of 536 reported cases, with 145 resulting in death. This gives a case fatality rate of 27 percent. The majority of cases have been recorded in the Middle East, specifically in countries such as Saudi Arabia, Jordan, and Qatar, as noted in ^[6]. It is crucial to maintain awareness of this disease and take appropriate precautions, particularly in areas where the virus has been reported, in order to prevent its spread.

The transmission of MERS-CoV between camels and humans is influenced by various environmental factors. The Hajj and Umrah pilgrimages are significant contributors to the spread of the virus, as these events attract more than 10 million individuals from different parts of the world to Saudi Arabia. Mathematical modeling has proven to be a valuable tool for understanding the outbreak, developing effective control strategies, and exploring the immune response to MERS-CoV. Several modeling studies have been conducted to investigate the MERS-CoV outbreak, as highlighted in ^[7,8,9]. One of the most extensive MERS-CoV epidemics was documented by Assire et al. in ^[10], who provided evidence suggesting that the virus can be transmitted from person to person. The Kingdom of Saudi Arabia (KSA) has reported the highest number of cases, with most of the cases being recorded there. It is imperative to take appropriate measures to mitigate the spread of the virus, especially in areas that have reported cases, in order to prevent further outbreaks. The consumption of unpasteurized camel milk, which is a common practice in KSA, is a potential cause of camel-to-human transmission of MERS-CoV, as suggested by ^[11]. Furthermore, Poletto et al. have proposed that the movement and mingling of individuals during the Hajj and Umrah events may play a significant role in the spread of MERS-CoV, as noted in ^[12]. Other activities, such as camel racing and the opening and closing of camel markets, have also been identified as potential contributors to the transmission of MERS-CoV. Several researchers have constructed mathematical models to study various diseases and real-world problems, including MERS-CoV, as mentioned in ^{[13,14,15,16]}. These models have proven to be useful in predicting the spread of the virus, designing effective control strategies, and exploring the immune response to MERS-CoV. It is important to continue this research in order to better understand the disease and limit its impact on public health.

This study utilizes the next-generation matrix (NGM) approach to model the transmission and spread of MERS-CoV between humans and camels. The researchers calculate the fundamental reproductive number and determine the local stability of the model using the Routh-Hurwitz (RH) criterion. Furthermore, the global stability of the model is assessed through the Castillo-Chavez and Lyapunov type function methods, and the stability conditions are determined in terms of $\mathcal{R}_0$. The parameters that impact the transmission of the disease are analyzed through sensitivity analysis of the fundamental reproductive numbers. On the other hand, there have been numerous effective studies have been conducted related to the modelling of infectious diseases ^[17,18,19], their stability analyses ^[20,21,22], bifurcation and chaos properties ^{[23,24,25,26,27,28,29,30]}.

In addition, the study employs an optimal control analysis to minimize the number of infected individuals and increase the number of cured individuals in the community. By identifying the factors that contribute to the spread of MERS-CoV, this research can inform effective control strategies and minimize the impact of the disease on public health.

2.   Formulation of the model

In this section, we introduce a transmission model for MERS-CoV that accounts for transmission between people-camel and human-human. The model is formulated using a set of differential equations that describe the dynamics of six distinct population groups. These groups include the susceptible population (S(t)), the exposed population (E(t)), the symptomatic and infectious population (I(t)), the asymptomatic but infectious population (A(t)), the hospitalized population (H(t)), and the recovered population (R(t)). To simplify the modeling process, we make the following four assumptions.

a. All the parameters and variables are non-negative.

b. Four transmission routes are considered for the disease transmission, which is from individuals symptomatic to asymptomatic, from which to hospitalize, and then reservoir, which are camels for MERS-CoV.

c. The rate of death because of MERS-CoV is considered in the compartment that contains the infection.

d. We suppose two types of recoveries, the first one is natural and the second one is with treatment.

Utilizing the above-considered assumptions, we obtain the non-linear system of ODEs as,

with the initial conditions

3.   Model well-posedness

Consider, $\mathrm{N}_p(\mathrm{t})$ which represents the total population of humans in such a manner that $\mathrm{N}_p(\mathrm{t}) = \mathrm{E}(\mathrm{t})+\mathrm{A}(\mathrm{t})+\mathrm{S}(\mathrm{t})+\mathrm{I}(\mathrm{t})+\mathrm{R}(\mathrm{t})+\mathrm{H}(\mathrm{t})$, then $\mathrm{N}_p(\mathrm{t})$ is bounded with lower bound to be 0 and the upper-bound $\frac{\phi}{\varpi_0}$, i.e., $0\leq \mathrm{N}_p(\mathrm{t})\leq \frac{\phi}{\varpi_0}$.

Using this fact, we present the following theorem:

Theorem 1. If $\mathrm{N}_p(\mathrm{t})$ represents the number of human and $0\leq \mathrm{N}_p(\mathrm{t})\leq \frac{\phi}{\varpi_0}$ and $\mathrm{N}_p(\mathrm{t})\leq\frac{\phi}{\varpi_0},$ then suggested model (2.1) is well defined in the region as follows:

Let us adopt $\mathbb{B}$ as Banach space, and positive $\mathrm{u} = t_+$, so

where the norm on the space $\mathbb{B}$ is supposed to be as $\|\Pi\| = \sum_{i = 1}^{7}\|\Pi_j\| = (\Pi_1, \Pi_2, \Pi_3, \Pi_4, \Pi_5, \Pi_6, \Pi_7)\in \mathbb{B}.$

Further, $\mathbb{B}_+$ represents cone($+ive$) of $\varpi^1(0, \mathrm{u}),$ so from Eq (3.1), $\mathbb{B}_+$ is given as

Hence the state space of system (2.1) yields:

We suppose an operator which is linear as ${\bf L}$ and vector ${\rm{ \mathsf{ ψ}}} = (\mathrm{S}, \mathrm{I}, \mathrm{A}, \mathrm{E}, \mathrm{H}, \mathrm{C}, \mathrm{R}),$ implies that ${\bf L}{\rm{ \mathsf{ ψ}}} = ({\bf L}_i)^T$, here $i = 1, 2, \ldots, 7$ where

and domain $\mathcal{D}({\bf L})$ is

$\mathcal{D}({\bf L}) = \bigg\{{\rm{ \mathsf{ ϕ}}}\in \mathbb{B}: {\rm{ \mathsf{ ψ}}} \in {\bf L}\mathcal{C}[0, \mathrm{u})$, ${\rm{ \mathsf{ ϕ}}}(0) = Ics\bigg\}$. Here, ${\bf L}\mathcal{C}[0, \mathrm{u})$ represent the set containing continuous functions which is defined on the $[0, \mathrm{u})$. Consider $\mathscr{O}$ is the nonlinear operator, that is $\mathscr{O}:B\rightarrow B$ defined as,

Suppose $V(\mathrm{t}) = (\mathrm{S}(\mathrm{t}), \mathrm{I}(\mathrm{t}), \mathrm{H}(\mathrm{t}), \mathrm{E}(\mathrm{t}), \mathrm{A}(\mathrm{t}), \mathrm{R}(\mathrm{t}), \mathrm{C}(\mathrm{t}))$ then the suggested system can be written as

where $V_{0} = (Ics)^T$. Utilizing the results in ^[31,32], we present the existence of the system's (3.2) solution, so we define following theorem:

Theorem 2. For each $V_{0}\in B_+,$ there arises an interval (maximal) $[0, t_0)$, and unique continuous solution $V(t, V_0)$, in such a way that,

Theorem 3. The suggested system (2.1) is invariant (positively) subjected to the non-negative $R^7+$.

Proof. Consider ${\rm{ \mathsf{ ψ}}}$ and $h_1 = (\xi+\varpi_0)$, $h_2 = ({\rm{ \mathsf{ σ}}}_1+{\rm{ \mathsf{ σ}}}_2+\varpi_0+\varpi_1)$, $h_3 = (\nu+\varpi_0)$, $h_4 = ({\rm{ \mathsf{ σ}}}_1+\nu)$, $h_5 = (\theta\varpi_0),$

It could be noted from Eq (3.3), matrix ${\bf D}$ is positive, while the off-diagonal of ${\bf L}$ are non-negative, so the properties of Metzler type matrix holds. Thus the suggested system is invariant in $R^7.$ □

Theorem 4. We assume a positive initial population value for the problem specified in Eq (2.2) and, if the solutions to the model in Eq (2.1) exist, they will be positive for all $\mathrm{u}$.

Proof. Let us consider the first equation

By constant formula of alternation, we obtain the solution (3.4),

□

$\mathrm{S}(\mathrm{t}) > 0$, in the same pattern one can present that, the remaining equations in (2.1) are positive.

4.   Stability analysis

The main focus of our study is to examine the mathematical and biological plausibility of the system described in Eq (2.1). To achieve this, we carry out a qualitative analysis of the system dynamics. Initially, we compute the threshold parameter $\mathcal{R}_0$, which is commonly referred to as the basic reproduction number. This metric allows us to evaluate the inherent capacity of the disease to spread, and determine whether or not an epidemic will persist or eventually fade out. Additionally, we investigate the equilibria of the system and discuss the factors that lead to system stability.

4.1.   Equilibria and basic reproductive number

In this study, we have performed a qualitative analysis of the suggested system in order to identify the conditions under which it remains stable. To achieve this, we have calculated the equilibria of the mathematical model described in Eq (2.1).

One of the most important equilibrium points for this system is the disease-free equilibrium (DFE), which represents the state of the system when no disease is present. In order to determine the DFE point for the system, we equate the right-hand side of the equations to zero, with the exception of the susceptible class $\mathrm{S}$, which we set to its initial value $\mathrm{S}_0$. By doing so, we are able to obtain the DFE point, which we denote as $\mathcal{D}_0$. This point represents an important baseline for the system, against which we can compare the behavior of the system in the presence of disease.

Overall, our qualitative analysis of the suggested system has allowed us to gain a deeper understanding of its behavior under various conditions, including the DFE point which is a key reference point for the system.

We utilize linear stability to study the dynamics of the DFE point and calculate the condition if the equilibrium point turns towards stability and the model becomes under control.

The endemic equilibrium (EE) point is expressed by $\mathcal{D}_1 = (\mathrm{S}^{*}, \mathrm{E}^{*}, \mathrm{I}^{*}, \mathrm{A}^{*}, \mathrm{H}^{*}, \mathrm{R}^{*}, \mathrm{C}^{*}),$ and it occurs in the presence of disease

The above equations present that, EE of the model (2.1) exists only, if $\mathcal{R}_0$ is greater than one. Thus we state the following theorem.

Theorem 5. The EE point $\mathcal{D}_1 = (\mathrm{S}^{*}, \mathrm{E}^{*}, \mathrm{I}^{*}, \mathrm{A}^{*}, \mathrm{H}^{*}, \mathrm{R}^{*}, \mathrm{C}^{*})$ exists only in case $\mathcal{R}_0$ is greater than one.

The definition of $(\mathcal{R}_0)$ can be described as the number of individuals who become infected after being in contact with an infected individual in a population that is initially fully susceptible and without any prior infections. If $\mathcal{R}_0 > 1$, it means that an epidemic is likely to occur, while if $\mathcal{R}_0 < 1$, an outbreak is unlikely. The value of $\mathcal{R}_0$ is crucial in determining the strength of control measures that need to be implemented to contain the epidemic. In order to calculate $\mathcal{R}_0$ for the suggested model (2.1), we use the method described in ^[33], we have

where $b_1 = (\xi+\varpi_0)$, $b_2 = ({\rm{ \mathsf{ σ}}}_1+{\rm{ \mathsf{ σ}}}_2+\varpi_0+\varpi_1)$, $b_3 = -\xi(1-{\rm{ \mathsf{ ρ}}})$, $b_4 = (\nu+\varpi_0)$ $b_5 = ({\rm{ \mathsf{ σ}}}_3+\varpi_0).$ $\mathcal{R}_0$ represents spectral-radius of NGM $\bar{H} = \mathcal{F}\mathcal{V}^{-1}$.

So $\mathcal{R}_0$ for model (2.1) is

where

The $\mathcal{R}_0$ of this model is composed of four components: transmission from individuals who are symptomatic to those who are asymptomatic, transmission from asymptomatic individuals to those who require hospitalization, transmission from hospitalization to the reservoir (camels for MERS-CoV), and transmission from the reservoir to susceptible individuals. These four modes of transmission collectively determine the risk of disease spread during this epidemic.

We study the dynamics of the proposed system (2.1) at DFE with aid of Theorem 6 as follows:

Theorem 6. The DFE point $\mathcal{D}_0 = (\mathrm{S}_0, \mathrm{I}_{0}, \mathrm{E}_{0}, \mathrm{H}_{0}, \mathrm{A}_{0}, \mathrm{C}_{0}, \mathrm{R}_{0}),$ is asymptotically stable (locally) if $\mathcal{R}_0\leq1$.

Proof. The Jacobian-matrix of the model at DFE point $(D_0, 0, 0, 0, 0, 0),$ is:

Characteristic equation of Jacobian matrix (4.2) is:

where

$a_1a_3a_2 > a_2^2a_4+a_2^2$ if $\mathcal{R}_0 < 1$. By RH criteria the real parts of all the roots for characteristic polynomial $P({\rm{ \mathsf{ ζ}}})$ are negative, which shows that $D_0$ is asymptotically local stable ^[34,35]. □

4.2.   Global stability analysis

The upcoming proof presents the global stability at DFE point $D_0.$ To analyze the global stability analysis at $F_0$ we introduce the Lyapunov function as follows.

Theorem 7. When the reproductive number $\mathcal{R}_0$ is less than 1, the disease-free equilibrium of the system is globally and asymptotically stable.

Proof. Consider the Lyapunov function as

Here $d_i$ for $i = 1, 2, 3, 4, 5$ are arbitrary constants, to be considered after differentiating Eq $(4.4)$, and using (2.1), so we obtain

By considering the $+ive$ parameter $d_1 = d_2 = d_3 = QQ_1$, $d_4 = \frac{1}{Q_2}$, $d_5 = \varpi_0$ and after the interpretation we obtain

where

$\mathfrak{U}'(\mathrm{t})$ is negative when $\mathrm{S} > \mathrm{S}^0$ and $\mathcal{R}_0 < 1$ and $\mathfrak{U}'(\mathrm{t}) = 0$ in case if $\mathrm{S} = \mathrm{S}^0$ by the LaSalle's invariance principle ^[36,37], and $\mathrm{E} = \mathrm{A} = \mathrm{I} = \mathrm{H} = \mathrm{C} = 0.$ Thus the DFE is globally asymptotic stable in $F_0.$ □

Theorem 8. If the threshold value is greater than 1, then the model (2.1) around EE point $\mathcal{D}_1$ is locally as well as globally asymptotically stable.

Proof. The linearization of model (2.1) around EE point $\mathcal{D}_1$ is,

where

Using row transformation, we obtain:

where $\mathfrak{B} = (\xi+\varpi_0)(\eta_1\mathrm{I}^*+\eta_2{\rm{ \mathsf{ ϕ}}} \mathrm{A}^*+\eta_3q\mathrm{H}^*+\eta_4\mathrm{S}^*),$ $\mathfrak{B}_1 = ({\rm{ \mathsf{ σ}}}_1+{\rm{ \mathsf{ σ}}}_2+\varpi_0+\varpi_1)(\xi+\varpi_0)\mathrm{A}^*, \mathfrak{B}_2(\nu+\varpi_0)(\xi+\varpi_0)(\mathcal{R}_0-1)\mathrm{A}^*,$ $\mathfrak{B}_3({\rm{ \mathsf{ σ}}}_1+{\rm{ \mathsf{ σ}}}_2+\varpi_0+\varpi_1)({\rm{ \mathsf{ σ}}}_3+\varpi_0),$ $B_4 = \varpi_2\xi{\rm{ \mathsf{ ρ}}}({\rm{ \mathsf{ σ}}}_1+{\rm{ \mathsf{ σ}}}_2+\varpi_0+\varpi_1).$

When $\mathcal{R}_0 > 1$ the real parts of eigenvalues are negative, hence the model (2.1) is asymptotically locally stable at $D_1$ ^[38]. □

Theorem 9. If $\mathcal{R}_0$ is greater than 1, then EE point $D_1$ is globally asymptotically stable and is not stable if less than 1.

Proof. In order to show the asymptotic global stability of the considered model (2.1) at EE point $D_1$, we utilize the Castillo-Chavez technique ^[39,40]. Now let us take the sub-system of (2.1),

Consider $\mathcal{P}$ and $\mathcal{P}^{\mid2\mid}$ be the linearized matrix and second-additive of the model which contains the first three equations of system (2.1), which becomes

Let $\mathbb{Q}(\chi) = { \mathbb{Q}}(\mathrm{S}(t), \mathrm{E}(t), \mathrm{I}(t)) = diag\left\{\frac{\mathrm{S}}{\mathrm{E}}, \frac{\mathrm{S}}{\mathrm{E}}, \frac{\mathrm{S}}{\mathrm{E}}\right\},$ then ${ \mathbb{Q}}^{-1}(\chi) = diag\left\{\frac{\mathrm{E}}{\mathrm{S}}, \frac{\mathrm{E}}{\mathrm{S}}, \frac{\mathrm{E}}{\mathrm{S}}\right\},$ the derivative of ${ \mathbb{Q}}_f(\chi)$ w.r.t time, implies that

Now ${ \mathbb{Q}}_f{ \mathbb{Q}}^{-1} = diag\{K_{1}, K_{1}, K_{1}\}$ and $\mathbb{Q}\mathcal{P}^{\mid2\mid}_2 \mathbb{Q}^{-1} = \mathcal{P}^{\mid2\mid}_2$, where $K_{1} = \dot{\mathrm{S}}{\mathrm{S}}-\frac{\dot{\mathrm{E}}}{\mathrm{E}}$ ${\bf A} = \mathbb{Q}_f \mathbb{Q}^{-1}+ \mathbb{Q}\mathcal{P}^{\mid2\mid}_2 \mathbb{Q}^{-1}$, and

Let $({n}_1, {n}_2, {n}_3)$ be vector $in$ $R^{3}$ and $\|.\|$ of $({n}_1, {n}_2, {n}_3)$ presented by,

Here we consider the Lozinski-measure introduced in ^[41], $\delta({\bf A})\leq sup\{{\varrho}_1, {\varrho}_2\} = sup\{\delta({\bf A}_{11})+\|{\bf A}_{12}\|, \delta({\bf A}_{22})+\|{\bf A}_{21}\|\},$ where $h_i = \delta({\bf A}_{ii})+\|{\bf A}_{ij}\|$ for $i = 1, 2$ and $i\neq j,$ $\Rightarrow$

where $\delta({\bf A}_{11}) = \frac{\dot{\mathrm{S}}}{\mathrm{S}}-\frac{\dot{\mathrm{E}}}{\mathrm{E}}-\eta_1\mathrm{I}-\eta_2{\rm{ \mathsf{ ϕ}}} \mathit{A}-\eta_3q\mathrm{H}-(\xi+\varpi_0)$, $\delta({\bf A}_{22}) = max \left\{\frac{\dot{\mathrm{S}}}{\mathrm{S}}-\frac{\dot{\mathrm{E}}}{\mathrm{E}}-\eta_1\mathrm{I}-\eta_2{\rm{ \mathsf{ ϕ}}} \mathit{A}-\eta_3 q\mathrm{H}-\varpi_0, \eta_1 \mathrm{I}-\eta_2{\rm{ \mathsf{ ϕ}}} \mathit{A}-\eta_3 q \mathrm{H}-\eta_4 \mathrm{C} \right\}\\ = \frac{\dot{\mathrm{S}}}{\mathrm{S}}-\frac{\dot{\mathrm{E}}}{\mathrm{E}}-\eta_1\mathrm{I}-\eta_2{\rm{ \mathsf{ ϕ}}} \mathit{A}-\eta_3 q\mathrm{H}-\eta_4 \mathrm{C}\}$, $\|{A} _{12}\| = \eta_1\mathrm{S}$ and $\|{\bf A}_{21}\| = max\{\xi{\rm{ \mathsf{ ρ}}}, 0\} = \xi{\rm{ \mathsf{ ρ}}}.$ Therefore ${\varrho}_1$ and ${\varrho}_2$ become, i.e, ${\varrho}_1\leq\frac{\dot{\mathrm{S}}}{\mathrm{S}}-2\varpi_0-\xi{\rm{ \mathsf{ ρ}}}$ and ${\varrho}_2\leq\frac{\dot{\mathrm{S}}}{\mathrm{S}}-2\varpi_0-\xi-min\{{\rm{ \mathsf{ σ}}}_1, \nu{\rm{ \mathsf{ ρ}}}\}$, which presents $\delta({\bf A})\leq\left\{\frac{\dot{\mathrm{S}}}{\mathrm{S}}+{\rm{ \mathsf{ σ}}}_1-min\{\xi, {\rm{ \mathsf{ σ}}}_1\}-2\varpi_0\right\}.$ Hence $\delta({B})\leq\frac{\dot{\mathrm{S}}}{\mathrm{S}}-2\varpi_0$. Integrating $\delta({\bf A})$ in $[0, t]$ and also considering $\lim_{t\rightarrow \infty},$ we have

So that the system of the first three compartments of model (2.1) is globally asymptotically stable. □

4.3.   Sensitivity analysis

To find the relation of parameters to ${\mathcal{R}_0}$ in the disease transmission we use the formula $\Delta^{\mathcal{R}_0}_h = \frac{\partial \mathcal{R}_0}{\partial k}\frac{h}{\mathcal{R}_0}$ where $h$ is the parameter, introduced by ^[34,35]. This makes it easy to identify the variables that have a substantial impact on reproduction number, using the above formula we have

These demonstrate the relevance of many factors in disease transmission. It also measures the change in ${\mathcal{R}_0}$ as a function of a parameter modification. The sensitivity indices show that there is indeed a direct relationship between ${\mathcal{R}_0}$ and a set of variables $S_1 = [{\eta_1, \eta_2, \eta_3, \eta_4, \phi, \psi_1, \psi_2}]$, while has an inverse relation with $\mathrm{S}_2 = [{\varpi_0, \varpi_2, {\rm{ \mathsf{ σ}}}_1, {\rm{ \mathsf{ σ}}}_2, \nu}]$. This demonstrates that higher the value of parameters $\mathrm{S}_1$ increases the value of threshold quantity greatly, but increasing the value for parameters $\mathrm{S}_ 2$ decreases the value of threshold value.

5.   Numerical simulation

In this part, we validate our analytical conclusion. We employ the Runge-Kutta technique of fourth order ^[42]. Some factors are chosen for demonstration purposes, while others are derived from publicly available data. The parameters are chosen in a way that is more biologically realistic. For the simulation, we use the following parameters. ${\rm{ \mathsf{ ϕ}}} = 0.00004; \eta_1 = 0.007; \varpi_0 = 0.0003; \eta_2 = 0.003;$ $\psi_2 = 0.00008; \varpi_1 = 0.0001; \eta_3 = 0.005; \xi = 0.002;$ $\eta_4 = 0.0001; {\rm{ \mathsf{ σ}}}_2 = 0.000001; \phi = 0.016; q = 0.00007; \varpi_2 = 0.00003; {\rm{ \mathsf{ σ}}}_1 = 0.001; {\rm{ \mathsf{ σ}}}_3 = 0.0007;$ $\psi_1 = 0.0006; \nu = 0.000002.$ Figures 3 and 4 depict the performance of the proposed model based on the aforementioned parameters, which validate the theorem's analytical discovery (4.2). According to the theoretical understanding of these findings, whenever $\mathrm{R}_0 < 1,$ each curve of solution of the sensitive population takes 150–300 days to achieve equilibrium. Likewise, the exposed community takes 250 to 150 days, the infected populace takes 200 to 100 days, and the asymptomatic populace, hospitalized, and recovered requires 100 to 50 days. Camel dynamics initially grow and then achieve an equilibrium state, as seen in illustration 0.

Next, we consider the parameters $\eta_1 = 0.007;$ ${\rm{ \mathsf{ σ}}}_1 = 0.001;$ $\varpi_1 = 0.0001;$ $\psi_1 = 0.0006;$ $\eta_2 = 0.003;$ $\varpi_2 = 0.00003;$ ${\rm{ \mathsf{ ϕ}}} = 0.00004;$ ${\rm{ \mathsf{ σ}}}_2 = 0.000001;$ $\eta_3 = 0.005;$ $\varpi_0 = 0.0003;$ $\eta_4 = 0.0001;$ $\xi = 0.002;$ $\phi = 0.016;$ ${\rm{ \mathsf{ σ}}}_3 = 0.0007;$ $q = 0.00007;$ $\nu = 0.000002.$ $\psi_2 = 0.00008;$ and find $\mathcal{R}_0 = 9.89887$ which is greater than $1$. We investigate the dynamics of the given model in the vicinity of the EE point. The numerical simulations based on the aforementioned settings are displayed in Figures (5) and (6), which validate the finding presented in Theorem (8).

6.   Optimal control (OC) strategy for the reduction of infected with MERS-CoV

We develop control techniques based on sensitivity as well as model dynamics (2.1). The maximal sensitivity indices parameter is ($\eta_1, \eta_2, \eta_3, \eta_4$), and increasing this value by 10 percent, raises the threshold value. To limit the progress of the illness, we must minimize these parameters by using the control variables $\mathbb{E}_1(\mathrm{t}), \mathbb{E}_2(\mathrm{t}), \mathbb{E}_3(\mathrm{t}), \mathbb{E}_4(\mathrm{t})$ to represent (awareness about the mask, isolation people (infected), oxygen therapies, ventilation and self-care from the camels.)

Our main goals are the reduction of MERS-CoV in the populace with increasing $\mathrm{R}(\mathrm{t})$ and decreasing $\mathrm{I}(\mathrm{t})$, $\mathit{A}(\mathrm{t})$ and $\mathrm{H}(\mathrm{t})$, reservoir $\mathrm{C}(\mathrm{t})$ with applying control parameters (time-dependent) $\mathbb{E}_1(\mathrm{t}), \mathbb{E}_2(\mathrm{t}), \mathbb{E}_3(\mathrm{t}), \mathbb{E}_4(\mathrm{t})$.

$i.$ $\mathbb{E}_1(\mathrm{t})$ represents the control parameter (time-dependent) represents the awareness concerning surgical masks and hand washing.

$ii.$ $\mathbb{E}_2(\mathrm{t})$ represents the control parameter (time-dependent) represents quarantining of infected persons.

$iii.$ $\mathbb{E}_3(\mathrm{t})$ represents the control parameter (time-dependent) represents mechanical ventilation (oxygen therapy).

$iv.$ $\mathbb{E}_4(\mathrm{t})$ represents the control parameter (time-dependent) self-care that is keeping distance from camels, avoiding raw milk, or eating improperly cooked meat.

By the use of these control parameters in our suggested optimal control problem which we obtain by modifying model (2.1):

with the initial conditions

The purpose here is to demonstrate that it is feasible to apply time-dependent control mechanisms while reducing the expense of doing so. We assume that the expenses of control schemes are nonlinear and take a quadratic shape, ^[43], which are cost variables that balance the size and importance of the sections of the optimization problem. As a result, we select the observable (cost) function as,

In Eq (6.2) $\zeta_1$, $\zeta_2$, $\zeta_3$, $\zeta_4$, $\zeta_5$, $\zeta_6$, $\zeta_7$, $\zeta_8$, stand for weight constants. $\zeta_1$, $\zeta_2$, $\zeta_3$, $\zeta_4$ express relative costs of infected ($\mathrm{I}$), asymptomatic ($\mathrm{A}$), hospitalized ($\mathrm{H}$) and reservoir ($\mathrm{C}$), while $\zeta_5$, $\zeta_6$, $\zeta_7$, $\zeta_8$ show associated-cost of control parameters. $\frac{1}{2}\zeta_5 \mathbb{E}_1^2$, $\frac{1}{2}\zeta_8 \mathbb{E}_4^2$, $\frac{1}{2}\zeta_6 \mathbb{E}_2^2$, $\frac{1}{2}\zeta_7 \mathbb{E}_3^2$, describe self care, treatment and isolation.

Our objective is to obtain OC pair $\mathbb{E}_1^{*}$, $\mathbb{E}_2^{*}$, $\mathbb{E}_3^{*}$, $\mathbb{E}_4^{*}$, i.e.,

dependent on model (6.1), we consider, the control-set of parameters as:

We obey the result ^[44], stating that the solution of the system exists in the case when control parameters are bounded as well as Lebesgue measurable. So, we consider that the suggested control system can be presented as:

In above system $\Omega = (\mathrm{S}, \mathrm{I}, \mathrm{E}, \mathrm{H}, \mathrm{A}, \mathrm{C})$, $\mathrm{A}(\Omega)$ and $\mathscr{B}(\Omega)$ represent linear and nonlinear bounded coefficient, respectively, so that

where $y_1 = (\nu+\varpi_0+ \mathbb{E}_1)$, $y_2 = ({\rm{ \mathsf{ σ}}}_1+{\rm{ \mathsf{ σ}}}_2+\varpi_0+\varpi_1+ \mathbb{E}_2)$, $y_3 = (\epsilon+\varpi_0+\varpi_2+ \mathbb{E}_2).$

Considering $L(\Omega) = F\Omega + A\Omega,$

Here $P = max(p_1, p_2, p_3, p_4, p_5, p_6, p_7, p_8)$ does not depend on the suggested model state-classes. We can also express

where $W = (P, \|A\|)$ is less than $\infty,$ $L$ is continuous in the Lipschitz sense, and from the description the system classes are non-negative, it obviously shows that the solution of model (6.1) exists. For the existence of the solution let us consider, and prove the following theorem:

Theorem 10. There exist an OC $\mathbb{E}^{*} = (\mathbb{E}_1^{*}, \mathbb{E}_2^{*}, \mathbb{E}_3^{*}, \mathbb{E}_4^{*})\in \mathbb{E}$, to control-system presented in Eqs (6.1) and (6.2).

Proof. As it is obvious that the control and system variables are not negative. It is also worth noting that U (set of variables) is closed and convex by expression. Furthermore, the control problem is bounded, indicating the problem's compactness. The expression $\zeta_1\mathrm{I}+\zeta_2\mathrm{A}+\zeta_3\mathrm{H}+\zeta_4\mathrm{C}+\frac{1}{2}(\zeta_5 \mathbb{E}_1^2(\mathrm{t})+\zeta_6 \mathbb{E}_2^2(\mathrm{t})+\zeta_7 \mathbb{E}_3^2(\mathrm{t})+\zeta_8 \mathbb{E}_4^2(\mathrm{t})$ is convex as well, w.r.t the set $U$. It guarantees the existence of OC for OC variables $(\mathbb{E}^{*}_1, \mathbb{E}^{*}_2, \mathbb{E}^{*}_3, \mathbb{E}^{*}_4)$. □

6.1.   Methods

Here, we determine the best solution to control problems (6.1) and (6.2). For this, we employ the Lagrangian, and Hamiltonian equations, as shown below:

To define the Hamiltonian (H) associated with the model, we use the notion $\Theta$ = $(\Theta_1, \Theta_2, \Theta_3, \Theta_4, \Theta_5, \Theta_6, \Theta_7)$ and $\Upsilon = (\Upsilon_1, \Upsilon_2, \Upsilon_3, \Upsilon_4, \Upsilon_5, \Upsilon_6, \Upsilon_7)$ then,

where

and $\mathfrak{Z}(x, u) = \Upsilon_1(x, u), \Upsilon_2(x, u), \Upsilon_3(x, u), \Upsilon_4(x, u), \Upsilon_5(x, u), \Upsilon_6(x, u), \Upsilon_7(x, u)$.

Here we utilize, the principle ^[45,46] to Hamiltonian, in order to obtain an optimality solution, which is stated that if the solution expressed with $(x^*, u^*)$ is optimal, then $\exists$ a function $\Theta$, such that

and the condition of transversality

Thus to obtain the adjoint variables and OC variables, we use the principles Eq (6.7). So we get

Theorem 11. Suppose that optimal and control-parameters are expressed by $\mathrm{S}^{*}$, $\mathrm{E}^{*}$, $\mathrm{A}^{*}$, $\mathrm{I}^{*}$, $\mathrm{H}^{*}$, $\mathrm{C}^{*}$, $\mathrm{R}^{*}$ be the optimal-state $(\mathbb{E}_1^{*}, \mathbb{E}_2^{*}, \mathbb{E}_3^{*}, \mathbb{E}_4^{*})$ for system (6.1)-(6.2). Then $\Theta\mathrm{A}(\mathrm{t})$ (adjoint variables set) satisfies:

having terminal condition

The OC variables $\mathbb{E}_1^{*}(\mathrm{t})$, $\mathbb{E}_2^{*}(\mathrm{t})$, $\mathbb{E}_3^{*}(\mathrm{t})$, $\mathbb{E}_4^{*}(\mathrm{t})$ are

Proof: The adjoint-model (6.9) is obtained by applying the principle (6.7) and the transversality conditions from the outcomes of $\Theta\mathrm{A}(\mathrm{t}) = 0$. For optimal functions set $\mathbb{E}_1^*, \mathbb{E}_2^*, \mathbb{E}_3^*$ and $\mathbb{E}_4^*,$ we utilized $\frac{\partial \mathrm{H}}{\partial u}$. In the next part, we evaluate the optimality problem numerically. Since it will be easier for such readers to understand as compared to analytical data. The optimization problem system is defined by its control-system (6.1), adjoint-system (6.9), boundary conditions, and OC functions.

6.2.   Numerical simulation of optimal control results

Using the RK technique of order four, we calculate the optimal control model (6.1) to observe the influence of masks, treatment, isolations, and self-care from camels. We employ the forward RK technique to get the solution of system (2.1) with starting conditions in the time interval [0, 50]. To obtain a solution to the adjoint-system (6.9), we apply the backward RK technique in the same domain with the assistance of the transversality constraint. For simulation purposes, we consider the parameters as: $\nu = 0.0071;$ $\varpi_1 = 0.014567125;$ $\eta_1 = 0.00041;$ $\eta_3 = 0.0000123;$ $\eta_4 = 0.0000123;$ $\theta = 0.98;$ ${\rm{ \mathsf{ σ}}}_1 = 0.0000404720925;$ ${\rm{ \mathsf{ σ}}}_3 = 0.00135;$ $q = 0.017816;$ $\psi_1 = 0.05;$ $\varpi_0 = 0.00997;$ $\eta_2 = 0.0000123;$ $\phi = 0.003907997;$ ${\rm{ \mathsf{ σ}}}_2 = 0.000431;$ ${\rm{ \mathsf{ ρ}}} = 0.00007;$ $\varpi_2 = 0.014567125$ $\psi_2 = 0.06;$. These parameters are considered in such a manner that more feasible biologically. Weight constants, here are taken as $\zeta_1 = 0.6610000;$ $\zeta_2 = 0.54450;$ $\zeta_3 = 0.0090030;$ $\zeta_4 = 0.44440;$ $\zeta_6 = 0.3550;$ $\zeta_7 = 0.67676;$ $\zeta_8 = 0.999.$ So we get the upcoming behaviours shown in Figures 7a, 7b, 7c, 7d, 8a, 8b, 8c.

These figures reflect the dynamics of susceptibles, exposed, infected people, asymptomatic persons, hospitalized, recovered, and MERS reservoirs, i.e., camels with and also without control. Our major goal in using the OC tool is to reduce the number of persons who are infected while increasing the number of those who are not infected, as demonstrated by numerical findings.

7.   Conclusions

In this study, we developed a mathematical model to analyze the transmission of MERS-CoV between people and its reservoir (camels), with the goal of assessing the transmission risk of MERS-CoV. We calculated the model's fundamental reproductive number, $\mathcal{R}_0$, and employed stability theory to examine the local and global behavior of the model and determine the conditions that lead to stability. We also evaluated the sensitivity of $\mathrm{R}_0$ to understand the impact of each epidemiological parameter on disease transmission. To minimize the number of infected individuals and intervention costs, we incorporated optimal control into the model, which included time-dependent control variables such as supportive care, surgical masks, treatment, and public awareness campaigns about the use of masks. Furthermore, our biological interpretation of the results indicates that if the basic reproduction number is less than one, the susceptible population decreases for up to 60 days, and then becomes stable, indicating that the population will remain stable. Our numerical simulations validated the effectiveness of our control strategies in reducing the number of infected individuals, asymptomatic cases, hospitalizations, and MERS-CoV reservoir, while increasing the susceptible and recovered populations. These simulations support our analytical work.

Acknowledgement

This study is supported via funding from Prince Sattam bin Abdulaziz University project number PSAU/2023/R/1444.