Mathematical Analysis of Fractal-Fractional Mathematical Model of COVID-19

Abstract

In this work, we modified a dynamical system that addresses COVID-19 infection under a fractal-fractional-order derivative. The model investigates the psychological effects of the disease on humans. We establish global and local stability results for the model under the aforementioned derivative. Additionally, we compute the fundamental reproduction number, which helps predict the transmission of the disease in the community. Using the Carlos Castillo-Chavez method, we derive some adequate results about the bifurcation analysis of the proposed model. We also investigate sensitivity analysis to the given model using the criteria of Chitnis and his co-authors. Furthermore, we formulate the characterization of optimal control strategies by utilizing Pontryagin’s maximum principle. We simulate the model for different fractal-fractional orders subject to various parameter values using Adam Bashforth’s numerical method. All numerical findings are presented graphically.

1. Introduction

Recently, the new coronavirus has become a severe global issue that is wreaking havoc on the global economy and the health of people. Human society, endowed with modern technology, has been currently infected with the COVID-19 disease, which is a dreadful infection that has mostly harmed underdeveloped countries more than developed countries. According to reports, the world’s first case of the new coronavirus was reported at the end of 2019. The study of infectious diseases is crucial to investigate the spread of various diseases in communities, and researchers have investigated the effect of these infectious diseases on the behavior of living organisms. It is worth noting that infectious diseases have affected humans, as well as other animals and plants, very badly for many centuries. For instance, see [1]. The pandemic has not only affected the world’s economy but also caused the loss of human lives. The main source of COVID-19 infection is droplets from the afflicted person’s mouth or nose when speaking or sneezing, and people who are in close proximity to the afflicted person are at risk of contracting the sickness. Almost every government in the world has embraced the lockdown policy as a preventive measure to safeguard the safety of its population, and in such cases, doctors and paramedics have pledged to provide health services to those who have been affected. Since the beginning of the COVID-19 pandemic, many researchers in their respective fields have studied coronavirus infection and attempted to identify possible ways to immunize humans against this virus in order to follow the World Health Organization’s (WHO) rules and regulations and immunize people to the maximum extent.

Since the onset of the COVID-19 pandemic, numerous scholars across various disciplines have conducted research on coronavirus infection. Determining the apex of infection within a particular country can inform decisions on how to mitigate these occurrences and the probability of future outbreaks with fewer infected individuals. The World Health Organization (WHO) has established guidelines and protocols to optimize the immunization of individuals. This investigation endeavors to identify individual perspectives on adhering to stay-at-home orders during the COVID-19 pandemic through the use of metaphors, and explores the correlation between these perspectives, stress, depression, and anxiety. It is important to acknowledge that measures taken to contain the spread of the pandemic may also result in diverse psychological complications.

A novel mathematical model has been created and analyzed to evaluate the interplay between COVID-19, Zika, chikungunya, and dengue co-dynamics. The aim is to determine the influence of COVID-19 on the dynamics of Zika, dengue, and chikungunya, as well as the reverse effect [2]. A proposed SEQAIHR model with saturated treatment is presented, and the biological feasibility of the model solutions is assessed, along with the calculation of the basic reproduction number (R0). The model displays transcritical, backward bifurcation, and forward bifurcation with hysteresis under certain restrictions with respect to different parameters. To validate the model, it is fitted with actual COVID-19 infection data from Hong Kong between December 19th, 2021 and April 3rd, 2022, and the model parameters are estimated. Sensitivity analysis is conducted to identify the most sensitive parameters that affect R0. R0 is estimated using actual initial growth data of COVID-19, and the effective reproduction number is calculated for the same period. Finally, an optimal control problem is proposed, taking into account effective vaccination and saturated treatment for the hospitalized class, to reduce the density of the infected class and minimize the implemented cost [3].

Quantitative data were analyzed using the Kruskal–Wallis test, while qualitative data was subjected to content analysis. The findings indicated that a significant proportion of the participants perceived staying at home as a state of confinement, ennui/dejection, and powerlessness. Conversely, some participants regarded staying at home as a conscientious act, a chance, and a prerequisite for ensuring safety. Participants who viewed staying at home as a form of confinement or a cause of ennui/dejection exhibited more psychological challenges, whereas those who perceived it as a responsibility or an opportunity experienced fewer difficulties.

Since its inception as a subject over a century ago, mathematical modeling and investigation of infectious diseases have been central to the study of infectious disease transmission. With the advancement of better processing, data software, and the availability of potential for transmitting and preserving information over the internet, electronic surveillance of infectious diseases has become common. Anderson et al. [4] found that mathematical models have been used as powerful tools for investigating epidemiological features of various diseases. Researchers have subsequently proposed appropriate control and preventative techniques to better understand infectious illnesses. Infectious illnesses and their transmission have been heavily influenced by human behavior and societal responses. Control efforts have been implemented to better understand the effects of human behavior and societal responses on the spread of infectious illnesses. The general population reacted to avian influenza in behavioral, psychological, and societal ways. Individuals’ psychological reactions throughout the epidemic affected the pace with which it spread, as well as the potential for emotional and societal difficulties.

In this paper, we develop a COVID-19 epidemic model that includes a nonlinear incidence rate for the spread of the COVID-19 virus from infective humans to susceptible humans, which accounts for the population saturation effect. Additionally, we consider a non-monotone incidence rate for the spread of the COVID-19 virus from infective humans to susceptible humans, which captures the psychological effect within humans.

Fractal-fractional calculus has gained much popularity in recent times. In this branch, the word ‘fractal derivative,’ also called the Hausdorff derivative, is a non-Newtonian generalization of the traditional order derivative. This derivative deals with the measurement of fractals, which are concepts defined in fractal geometry. The concept was introduced when the classical derivative failed to explain the anomalous behavior of the diffusion process in the fractal nature of the media. This derivative is local in nature compared to the ordinary fractional order derivative, which is non-local. It is worth noting that certain phenomena related to porous media, turbulence, aquifers, and other media usually display fractal properties. The use of the fractal-fractional derivative has become increasingly popular in recent years due to its ability to model complex systems with high accuracy. This type of derivative takes into account the fractal nature of many natural phenomena, such as the branching patterns of trees or the irregular shapes of coastlines. By incorporating fractal geometry into the traditional calculus framework, the fractal-fractional derivative provides a more comprehensive understanding of the behavior of complex systems. This has important applications in fields such as physics, engineering, and finance, where accurate modeling of complex systems is crucial. Additionally, the use of the fractal-fractional derivative has the potential to lead to new insights and discoveries in these fields, making it an exciting area of research.

The mentioned form of the derivative has been used effectively in recent years to investigate various real-world problems. For some recent work, we refer to [5,6,7,8,9,10,11]. Therefore, keeping this in mind, we also investigate our proposed model in terms of the Caputo fractal-fractional order derivative. We establish a detailed analysis, including the construction of basic reproductive numbers, and global and local stability analysis using various approaches. We also use mathematical arguments to develop sufficient conditions for bifurcation and sensitivity analysis for the proposed model. Finally, we present various numerical simulations using the Adam Bashforth method, and all the results are displayed graphically.

The manuscript is structured as follows: In Section 2, we present fundamental definitions and theorems that underpin this study. In Section 3, we determine the disease-free and endemic equilibrium points and calculate the basic reproduction number at the disease-free equilibrium point. Additionally, we conduct a stability analysis based on the equilibrium point. In Section 4, we compute the saddle node resulting from the zero eigenvalue and derive several outcomes concerning bifurcation analysis. In Section 5, we perform sensitivity analysis based on the basic reproduction number ${\mathcal{R}}_0$. In Section 6, we apply vaccination and treatment control variables to prevent the population from contracting the COVID-19 virus again. In Section 7, we establish the numerical scheme for the fractal-fractional problem from [12] for the numerical simulation of the COVID-19 model. The numerical values and initial conditions to the compartments/classes with descriptions are given in Table 1 and Table 2. Finally, we present the outcomes and discussion and summarize the numerical simulations and the paper in Section 8 and Section 9, respectively.

2. Preliminaries

In this section, we recall fractional-fractional operators such that:

Theorem 1

([9]). Assume that the function $s\left(\mathfrak{t}\right)$ is differentiable in the open interval $(a,b)$. Then, the fractal-fractional derivative of s of order λ in the Caputo sense with power law is given as:

(1)$\begin{array}{c}\hfill \begin{array}{cc}\hfill {}_a^{\mathbf{FFP}}{\mathcal{D}}_{\mathfrak{t}}^{\kappa ,\lambda }s\left(\mathfrak{t}\right)=& \frac{1}{\mathsf{\Gamma }[m-\lambda ]}{\int }_a^{\mathfrak{t}}\frac{ds\left(\xi \right)}{d{\xi }^{\kappa }}{(\mathfrak{t}-\xi )}^{m-\lambda -1}d\xi ,m-1<\lambda \le m,0<m-1<\kappa \le m\hfill \\ \hfill \frac{ds\left(\xi \right)}{d{\xi }^{\kappa }}=& \underset{\mathfrak{t}-\xi }{lim}\frac{s\left(\mathfrak{t}\right)-s\left(\xi \right)}{{\mathfrak{t}}^{\kappa }-{\xi }^{\kappa }}.\hfill \end{array}\end{array}$

 if s is fractal differentiable on $(a,b)$ with order κ.

Theorem 2

([9]). The fractal-fractional integral of function s with fractal order $\kappa >0$ and fractional order $\alpha >0$ is defined as:

(2)$\begin{array}{c}\hfill {}_0^{\mathbf{FFP}}{\mathcal{I}}_{\mathfrak{t}}^{\kappa ,\lambda }s\left(\mathfrak{t}\right)=\frac{\kappa }{\mathsf{\Gamma }\left(\lambda \right)}{\int }_0^{\mathfrak{t}}{\xi }^{\lambda -1}s\left(\xi \right){(\mathfrak{t}-\xi )}^{\lambda -1}d\xi .\end{array}$

 if $s\left(\mathfrak{t}\right)$ is continuous in an opened interval $[a,b]$.

Definition 1.

Suppose we consider the fractal-fractional non-linear ODE such that

(3) $\begin{array}{cc}\hfill {}_0^{\mathbf{FFP}}{\mathcal{D}}^{\kappa ,\lambda }u\left(t\right)& =v(t,u\left(t\right)),withu\left(0\right)=u_0.\hfill \end{array}$

From Zhang et al. [12] the resultant numerical scheme for the Equation (3) can be written as:

(4) $\begin{array}{c}\hfill \begin{array}{cc}\hfill u(k+1)=& u\left(0\right)+\frac{\kappa {\hslash }^{\lambda }}{\mathsf{\Gamma }(\lambda +2)}\sum _{p=0}^m\left[t_{\left(p\right)}^{\kappa -1}v\left(u_{\left(p\right)},t_p\right)\right.\hfill \\ \hfill & \times \left({(m+1-p)}^{\kappa }(m-p+2+\kappa )-{(m-p)}^{\kappa }(m-p+2+2\kappa )\right)\hfill \\ \hfill & \left.-t_{(p-1)}^{\kappa -1}v\left(u_{(p-1)},t_{(p-1)}\right)\left({(m+1-p)}^{\kappa }+1-{(m-p)}^{\kappa }(m-p+1+\kappa )\right)\right].\hfill \end{array}\end{array}$

Formulation of the Proposed Model of COVID-19

The model is modified from the existence SIQR model of [16]:

(5)$\begin{array}{c}\hfill \left.\begin{array}{cc}\hfill \frac{dS}{dt}& =\mathsf{\Lambda }-\mu S-\frac{\beta SI}{N},\hfill \\ \hfill \frac{dI}{dt}& =\frac{\beta SI}{N}-(\mu +\gamma +\delta +\alpha )I,\hfill \\ \hfill \frac{dQ}{dt}& =\delta I-(\mu +ϵ+\alpha )Q,\hfill \\ \hfill \frac{dR}{dt}& =\gamma I+ϵQ-\mu R.\hfill \end{array}\right\}\end{array}$

The modified version of the model (5) is as follows in the fractal-fractional Caputo sense with power law kernel:

(6)$\begin{array}{c}\hfill \left.\begin{array}{cc}\hfill {}_0^{\mathbf{FFP}}{\mathcal{D}}^{\kappa ,\lambda }S\left(t\right)& =\upsilon -\tau S\left(t\right)-\frac{\vartheta S\left(t\right)I\left(t\right)}{1+\omega I\left(t\right)},\hfill \\ \hfill {}_0^{\mathbf{FFP}}{\mathcal{D}}^{\kappa ,\lambda }I\left(t\right)& =\frac{\vartheta S\left(t\right)I\left(t\right)}{1+\omega I\left(t\right)}-(\tau +\chi +\nu +{\varrho }_1)I\left(t\right),\hfill \\ \hfill {}_0^{\mathbf{FFP}}{\mathcal{D}}^{\kappa ,\lambda }Q\left(t\right)& =\nu I\left(t\right)-(\tau +\varpi +{\varrho }_2)Q\left(t\right),\hfill \\ \hfill {}_0^{\mathbf{FFP}}{\mathcal{D}}^{\kappa ,\lambda }\mathcal{P}\left(t\right)& =\chi I\left(t\right)+\varpi Q\left(t\right)-\tau \mathcal{P}\left(t\right),\hfill \\ \hfill {}_0^{\mathbf{FFP}}{\mathcal{D}}^{\kappa ,\lambda }{D}_{IQ}& ={\varrho }_{1}I\left(t\right)+{\varrho }_{2}Q\left(t\right).\hfill \end{array}\right\}\end{array}$

where the parameters are defined as: $\upsilon$ = recruitment rate, $\tau$ natural death rate, $\vartheta$ transmission rate, $\omega$ psychological effect on humans [14,17], ${\varrho }_1,{\varrho }_2$ death rates of infected and quarantine population due to disease, respectively, $\chi ,\varpi$ are recovery rates of infected and quarantine population, respectively, and $\nu$ is the rate of infection. The functions or classes are defined as: $S\left(t\right)$ is the susceptible human population, $I\left(t\right)$ is the infected human population, $Q\left(t\right)$ is the quarantine population, $\mathcal{P}\left(t\right)$ is the protected population, and $D_{IQ}$ represents deaths in infected and quarantine populations.

3. Equilibrium Points and Stability

Here we derive our main results in the following sub-sections.

3.1. Equilibrium Points

We have two types of equilibrium points such as disease-free and endemic equilibrium points. In the disease-free equilibrium point, we have no infection in the population. Thus in the mathematical model the infected compartments become zero while in endemic equilibrium the model counts as a whole. We denote the disease-free equilibrium point by ${\mathcal{E}}^0$ while the disease-endemic equilibrium point is denoted by ${\mathcal{E}}^{*}$. The equilibrium points are as follows:

(7)$\begin{array}{c}\hfill \begin{array}{c}\hfill {\mathcal{E}}^0=\left(\frac{\nu }{\tau },0,0,0,0\right).\end{array}\end{array}$

while the disease-endemic equilibrium is computed in terms of one class, such that

(8)$\begin{array}{c}\hfill \begin{array}{c}\hfill {\mathcal{E}}^{*}=\left({\mathcal{S}}^{*},{\mathcal{I}}^{*},{\mathcal{Q}}^{*},{\mathcal{P}}^{*},{\mathcal{D}}^{*}\right),\end{array}\end{array}$

where

(9)$\begin{array}{c}\hfill \left.\begin{array}{cc}\hfill {\mathcal{S}}^{*}& =\frac{\nu (1+\omega {\mathcal{I}}^{*})}{(1+\omega {\mathcal{I}}^{*})+\vartheta {\mathcal{I}}^{*}},\hfill \\ \hfill {\mathcal{Q}}^{*}& =-\frac{{\varrho }_1}{{\varrho }_2}{\mathcal{I}}^{*},\hfill \\ \hfill {\mathcal{P}}^{*}& =\frac{1}{\tau }\left(\chi \tau -\frac{\varpi }{{\varrho }_2}{\varrho }_1{\mathcal{I}}^{*}\right).\hfill \end{array}\right\}\end{array}$

3.2. Basic Reproduction Number

Here we are going to compute the fundamental reproduction number as

(10)${\left.\frac{dZ}{dt}\right|}_{{\mathcal{E}}^0}=\mathbf{f}-\mathbf{v}.$

where the non-linear and linear terms from the infected classes in matrix $\mathbf{f}$ and $\mathbf{v}$, respectively are given as

(11)$\begin{array}{c}\hfill \mathbf{f}=\left(\begin{array}{c}\frac{\vartheta S\left(t\right)I\left(t\right)}{1+\omega I\left(t\right)}\\ 0\end{array}\right),\mathbf{v}=\left(\begin{array}{c}-(\tau +\chi +\nu +{\varrho }_1)I\left(t\right)\\ (\tau +\varpi +{\varrho }_2)Q\left(t\right)-\nu I\left(t\right)\end{array}\right).\end{array}$

Now, the Jacobian matrix of $\mathbf{f}$ and $\mathbf{v}$ is given by

(12)$\begin{array}{c}\hfill F=\left(\begin{array}{cc}\frac{\vartheta S^0\left(t\right)I^0\left(t\right)}{1+\omega I^0\left(t\right)}-\frac{\vartheta S^0\left(t\right)I^0\left(t\right)\omega \vartheta }{{(1+\omega I^0\left(t\right))}^2}& 0\\ 0& 0\end{array}\right),\mathrm{and}V=\left(\begin{array}{cc}\tau +\chi +\nu +{\varrho }_1& 0\\ -\nu & \tau +\varpi +{\varrho }_2\end{array}\right).\end{array}$

We now calculate the inverse of matrix V and the next generation matrix $G\left({\mathcal{E}}^0\right)$, such that

(13)$\begin{array}{c}\hfill V^{-1}=\left(\begin{array}{cc}\frac{1}{\tau +\chi +\nu +{\varrho }_1}& 0\\ \frac{\nu }{({\varrho }_1+\tau +\varpi )(\tau +\chi +\nu +{\varrho }_1)}& \frac{1}{{\varrho }_2+\tau +\varpi }\end{array}\right),\mathrm{and}G\left({\mathcal{E}}^0\right)=\left(\begin{array}{cc}\frac{\vartheta S^0}{(\tau +\chi +\nu +{\varrho }_1)}& 0\\ 0& 0\end{array}\right).\end{array}$

Thus, the non-zero and largest eigenvalue is the basic reproduction number ${\mathcal{R}}_0$ is

(14)${\mathcal{R}}_0=\frac{\vartheta S^0}{(\tau +\chi +\nu +{\varrho }_1)},$

where $S^0=\frac{\upsilon }{\tau }.$

3.3. Stability Analysis

Let a fractal-fractional order linear homogenous system of the following form be given:

(15)$\begin{array}{c}\hfill \begin{array}{cc}\hfill {}_0^{\mathbf{FFP}}{\mathcal{D}}_{\mathfrak{t}}^{\kappa ,\lambda }s\left(t\right)& =\mathcal{A}s\left(\mathfrak{t}\right),\hfill \\ \hfill s\left(0\right)& =s_0.\hfill \end{array}\end{array}$

where $\mathcal{A}\in N_{n\times n}\left(\mathcal{R}\right)$ & $0<\kappa ,\lambda \le 1$. The following theorems are on the stability of a linear homogenous system (15).

Theorem 3.

The COVID-19 model at the disease-free equilibrium point ${\mathcal{E}}_0$ is locally asymptotically stable if ${\mathcal{R}}_0<1$, otherwise unstable.

Proof. 

The Jacobian matrix of the system (6) at disease-free equilibrium point ${\mathcal{E}}^0$ is given by

(16)$\begin{array}{c}\hfill \mathcal{J}\left({\mathcal{E}}^0\right)=\left(\begin{array}{ccccc}-\tau & -\vartheta S^0& 0& 0& 0\\ 0& \vartheta S^0-({\varrho }_1+\nu +\tau +\chi )& 0& 0& 0\\ 0& \upsilon & -({\varrho }_2+\tau +\varpi )& 0& 0\\ 0& \chi & \varpi & -\tau & 0\\ 0& {\varrho }_1& {\varrho }_2& 0& 0\end{array}\right).\end{array}$

The two eigenvalues are ${\zeta }_1={\zeta }_2=-\tau$. The reduced matrix becomes

(17)$\begin{array}{c}\hfill \mathcal{J}\left({\mathcal{E}}^0\right)=\left(\begin{array}{ccc}\vartheta S^0-({\varrho }_1+\nu +\tau +\chi )& 0& 0\\ \upsilon & -({\varrho }_2+\tau +\varpi )& 0\\ {\varrho }_1& {\varrho }_2& 0\end{array}\right).\end{array}$

The characteristic equation becomes

(18)$\begin{array}{c}\hfill \zeta ({\zeta }^2+{\alpha }_1\zeta +{\alpha }_2)=0.\end{array}$

Thus, the roots of the characteristic Equation (18) are given by

(19)$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\zeta }_1& =-({\varrho }_1+\nu +\tau +\chi )(1-{\mathcal{R}}_0),\hfill \\ \hfill {\zeta }_2& =-({\varrho }_2+\tau +\varpi ),\hfill \\ \hfill {\zeta }_3& =0.\hfill \end{array}\end{array}$

We have a single zero eigenvalue while according to the eigenvalues the system is locally asymptotically stable. For the zero eigenvalue, we proceed with the saddle-node bifurcation in the next section. For this, we take basic reproduction number ${\mathcal{R}}_0$ equal to 1 such that

(20)${\mathcal{R}}_0=\frac{\vartheta S^0}{(\tau +\chi +\nu +{\varrho }_1)}=1.$

From Equation (20), we have the parameter ${\vartheta }^0=\vartheta$, that is

(21)${\vartheta }^0=\frac{\tau (\tau +\chi +\nu +{\varrho }_1)}{\nu }.$

Thus, Equation (21) is the bifurcation parameter. □

4. Bifurcation Analysis

In this section, we demonstrate bifurcation analysis of the system (6). We follow theorem A of Carlos Castillo-Chavez and Song [18] and Buonomo and Lacitignola [19], which is discussed above to obtain the required result. Recall the system (6):

(22)$\begin{array}{c}\hfill \left.\begin{array}{cc}\hfill {}_0^{\mathbf{FFP}}{\mathcal{D}}^{\kappa ,\lambda }\mathcal{S}\left(t\right)& =\upsilon -\tau S\left(t\right)-\frac{\vartheta S\left(t\right)I\left(t\right)}{1+\omega I\left(t\right)},\hfill \\ \hfill {}_0^{\mathbf{FFP}}{\mathcal{D}}^{\kappa ,\lambda }\mathcal{I}\left(t\right)& =\frac{\vartheta S\left(t\right)I\left(t\right)}{1+\omega I\left(t\right)}-(\tau +\chi +\nu +{\varrho }_1)I\left(t\right),\hfill \\ \hfill {}_0^{\mathbf{FFP}}{\mathcal{D}}^{\kappa ,\lambda }\mathcal{Q}\left(t\right)& =\nu I\left(t\right)-(\tau +\varpi +{\varrho }_2)Q\left(t\right),\hfill \\ \hfill {}_0^{\mathbf{FFP}}{\mathcal{D}}^{\kappa ,\lambda }\mathcal{P}\left(t\right)& =\chi I\left(t\right)+\varpi Q\left(t\right)-\tau \mathcal{P}\left(t\right),\hfill \\ \hfill {}_0^{\mathbf{FFP}}{\mathcal{D}}^{\kappa ,\lambda }{\mathcal{D}}_{IQ}& ={\varrho }_{1}I\left(t\right)+{\varrho }_{2}Q\left(t\right).\hfill \end{array}\right\}\end{array}$

Now, put $\mathcal{S}\left(t\right)=u_1$, $\mathcal{I}\left(t\right)=u_2$, $\mathcal{Q}\left(t\right)=u_3$, $\mathcal{P}\left(t\right)=u_4$, & ${\mathcal{D}}_{IQ}\left(t\right)=u_5$, while the right hand side equals $g_k$ for $k=1,2,3,4,5$. Thus the system (22) becomes

(23)$\begin{array}{c}\hfill \left.\begin{array}{cc}\hfill g_1& =\upsilon -\tau u_1-\frac{\vartheta u_1{u}_2}{1+\omega u_2},\hfill \\ \hfill g_2& =\frac{\vartheta u_1{u}_2}{1+\omega u_2}-(\tau +\chi +\nu +{\varrho }_1)u_2,\hfill \\ \hfill g_3& =\nu u_2-(\tau +\varpi +{\varrho }_2)u_3,\hfill \\ \hfill g_4& =\chi u_2+\varpi u_3-\tau u_4,\hfill \\ \hfill g_5& ={\varrho }_1{u}_2+{\varrho }_2{u}_3.\hfill \end{array}\right\}\end{array}$

Theorem 4.

By applying Theorem A we show that system (6) can exhibit a forward bifurcation, when $\vartheta ={\vartheta }^0$.

Proof. 

According to Equation (20), the disease-free equilibrium $E^0$ exhibits a single zero eigenvalue, while the remaining eigenvalues are either real or possess negative real components. Consequently, when $\vartheta ={\vartheta }^0$ (or $R_0=1$), $E^0$ represents a non-hyperbolic equilibrium, thereby satisfying assumption (A1) of Theorem A. Let $w={\left(w_1,w_2,w_3,w_4,w_5\right)}^T$ denote a right eigenvector linked to the zero eigenvalue ${\zeta }_3=0$. It follows that

(24)$\begin{array}{c}\hfill \left(\begin{array}{c}-\tau w_1-\frac{\vartheta w_2\nu }{\tau }\\ [\vartheta \frac{\nu }{\tau }-({\varrho }_1+\nu +\tau +\chi )]w_2\\ \nu w_2-({\varrho }_2+\tau +\varpi )w_2\\ \chi w_2+\varpi w_3-\tau w_4\\ {\varrho }_1{w}_2+{\varrho }_2{w}_3\end{array}\right)=\left(\begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\end{array}\right).\end{array}$

Thus, the non-zero ws are given as

(25)$\begin{array}{c}\hfill w_1=\frac{\vartheta w_3{S}^0}{\tau {\varrho }_1},w_2=-\frac{{\varrho }_2{w}_3}{{\varrho }_1},\mathrm{and}{w}_4=\left(\varpi -\frac{\chi {\varrho }_2}{\tau }\right)w_3.\end{array}$

Furthermore, the left eigenvector $v=\left(v_1,v_2,v_3,v_4,v_5\right)$ satisfying $v·w=1$ is given by:

(26)$\begin{array}{c}\hfill \left(\begin{array}{c}-\tau v_1\\ [\vartheta S^0-({\varrho }_1+\nu +\tau +\chi )]v_2-\vartheta S^0{v}_1+\nu v_3+\chi v_4+{\varrho }_1{v}_5\\ \varpi v_4-({\varrho }_2+\tau +\varpi )v_3+{\varrho }_2{v}_5\\ -\tau v_4\\ 0\end{array}\right)=\left(\begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\end{array}\right).\end{array}$

Thus, the non-zero vs are given as

(27)$\begin{array}{c}\hfill v_2=\frac{v_3[\nu +\frac{{\varrho }_2}{{\varrho }_1}({\varrho }_2+\tau +\varpi )]}{({\varrho }_1+\nu +\tau +\chi )-\vartheta S^0},v_5=v_3\left[\frac{{\varrho }_2+\tau +\varpi }{{\varrho }_2}\right].\end{array}$

The coefficients a and b defined in Theorem A are given by

(28)$\begin{array}{c}\hfill \left.\begin{array}{cc}\hfill a& =\sum _{k,i,j=1}^5{v}_k{w}_i{w}_j\frac{{\partial }^2{g}_k}{\partial u_i\partial u_j}\left(E_0,{\vartheta }^0\right),\hfill \\ \hfill b& =\sum _{k,i=1}^5{v}_k{w}_i\frac{{\partial }^2{g}_k}{\partial u_i\partial \vartheta }\left(E_0,{\vartheta }^0\right).\hfill \end{array}\right\}\end{array}$

Now, from Equations (23), (25), (27) and (28) we have

(29)$\begin{array}{c}\hfill \begin{array}{cc}\hfill a& =v_2{w}_1^2\frac{{\partial }^2{g}_2}{\partial u_1^2}\left(E^0,{\vartheta }^0\right)+v_2{w}_1{w}_2\frac{{\partial }^2{g}_2}{\partial u_1\partial u_2}\left(E^0,{\vartheta }^0\right)+v_2{w}_1{w}_3\frac{{\partial }^2{g}_2}{\partial u_1\partial u_3}\left(E^0,{\vartheta }^0\right)\hfill \\ \hfill & +v_2{w}_1{w}_4\frac{{\partial }^2{g}_2}{\partial u_1\partial u_4}\left(E^0,{\vartheta }^0\right)+v_3{w}_2{w}_1\frac{{\partial }^2{g}_3}{\partial u_2\partial u_1}\left(E^0,{\vartheta }^0\right)+v_3{w}_2^2\frac{{\partial }^2{g}_3}{\partial u_2^2}\left(E^0,{\vartheta }^0\right)\hfill \\ \hfill & +v_3{w}_2{w}_3\frac{{\partial }^2{g}_3}{\partial u_2\partial u_3}\left(E^0,{\vartheta }^0\right)+v_3{w}_2{w}_4\frac{{\partial }^2{g}_3}{\partial u_2\partial u_4}\left(E^0,{\vartheta }^0\right)+v_5{w}_3{w}_1\frac{{\partial }^2{g}_5}{\partial u_5\partial u_1}\left(E^0,{\vartheta }^0\right)\hfill \\ \hfill & +v_5{w}_3{w}_2\frac{{\partial }^2{g}_5}{\partial u_5\partial u_2}\left(E^0,{\vartheta }^0\right)+v_5{w}_3{w}_3\frac{{\partial }^2{g}_5}{\partial u_5\partial u_3}\left(E^0,{\vartheta }^0\right)+v_5{w}_3{w}_4\frac{{\partial }^2{g}_5}{\partial u_5\partial u_4}\left(E^0,{\vartheta }^0\right),\hfill \end{array}\end{array}$

which yields that

(30)$\begin{array}{c}\hfill \begin{array}{c}\hfill a=v_3\tau {\left[\frac{{\varrho }_2{w}_3(\tau +\chi +\nu +{\varrho }_1)}{{\varrho }_1}\right]}^2\left[\frac{\nu +\frac{{\varrho }_2}{{\varrho }_1}({\varrho }_2+\tau +\varpi )}{({\varrho }_1+\nu +\tau +\chi )-\frac{\vartheta \nu }{\tau }}\right],\end{array}\end{array}$

and b is given as

(31)$\begin{array}{c}\hfill \begin{array}{cc}\hfill b& =v_2{w}_1\frac{{\partial }^2{g}_2}{\partial u_1\partial \vartheta }\left(E^0,{\vartheta }^0\right)+v_2{w}_2\frac{{\partial }^2{g}_2}{\partial u_2\partial \vartheta }\left(E^0,{\vartheta }^0\right)+v_2{w}_3\frac{{\partial }^2{g}_2}{\partial u_3\partial \vartheta }\left(E^0,{\vartheta }^0\right)\hfill \\ \hfill & +v_2{w}_4\frac{{\partial }^2{g}_2}{\partial u_4\partial \vartheta }\left(E^0,{\vartheta }^0\right)+v_3{w}_1\frac{{\partial }^2{g}_3}{\partial u_1\partial \vartheta }\left(E^0,{\vartheta }^0\right)+v_3{w}_2\frac{{\partial }^2{g}_3}{\partial u_2\partial \vartheta }\left(E^0,{\vartheta }^0\right)\hfill \\ \hfill & +v_3{w}_3\frac{{\partial }^2{g}_3}{\partial u_3\partial \vartheta }\left(E^0,{\vartheta }^0\right)+v_3{w}_4\frac{{\partial }^2{g}_3}{\partial u_4\partial \vartheta }\left(E^0,{\vartheta }^0\right)+v_5{w}_1\frac{{\partial }^2{g}_5}{\partial u_1\partial \vartheta }\left(E^0,{\vartheta }^0\right)\hfill \\ \hfill & +v_5{w}_2\frac{{\partial }^2{g}_5}{\partial u_2\partial \vartheta }\left(E^0,{\vartheta }^0\right)+v_5{w}_2\frac{{\partial }^2{g}_5}{\partial u_2\partial \vartheta }\left(E^0,{\vartheta }^0\right)+v_5{w}_3\frac{{\partial }^2{g}_5}{\partial u_3\partial \vartheta }\left(E^0,{\vartheta }^0\right)\hfill \\ \hfill & +v_5{w}_4\frac{{\partial }^2{g}_5}{\partial u_4\partial \vartheta }\left(E^0,{\vartheta }^0\right).\hfill \end{array}\end{array}$

After evaluation of b, such that

(32)$\begin{array}{c}\hfill \begin{array}{c}\hfill b=\frac{{\varrho }_2{v}_3{w}_3\left[\nu +\frac{{\varrho }_2}{{\varrho }_1}({\varrho }_2+\tau +\varpi )\right]}{{\varrho }_1\left[({\varrho }_1+\nu +\tau +\chi )-\frac{\vartheta \nu }{\tau }\right]}\left[\omega -\frac{\nu }{\tau }\right],\end{array}\end{array}$

eventually, a and b are both positive. Hence, it is backward bifurcation.  □

5. Sensitivity Analysis

Here, we use the criterion proposed by Chitnis et al. [20] to determine how sensitive the parameters in the proposed model are, in order to investigate the transmission of infections or disorders. We compute the partial derivative with regard to each component in the reproduction number in order to ascertain the sensitivity of ${\mathcal{R}}_0$ in terms of each of its parameters. In other words, if p is an arbitrary parameter, then $R_0$ versus p sensitivity index is determined as follows

${\mathbb{S}}_p^{{\mathcal{R}}_0}=\frac{p}{{\mathcal{R}}_0}\left[\frac{\partial {\mathcal{R}}_0}{\partial p}\right].$

Now, according to the above relation, we have

(33)$\begin{array}{c}\hfill \left.\begin{array}{cc}\hfill {\mathbb{S}}_{\nu }^{{\mathcal{R}}_0}& =-\frac{\nu }{{\mathcal{R}}_0}\left[\frac{\vartheta }{\tau {\left(\chi +{\varrho }_1+\nu +\tau \right)}^2}\right]<0,\hfill \\ \hfill {\mathbb{S}}_{\tau }^{{\mathcal{R}}_0}& =-\frac{\tau }{{\mathcal{R}}_0}\left[\frac{\upsilon \vartheta }{\tau {\left(\chi +{\varrho }_1+\nu +\tau \right)}^2}+\frac{\upsilon \vartheta }{{\tau }^2\left(\chi +{\varrho }_1+\nu +\tau \right)}\right]<0,\hfill \\ \hfill {\mathbb{S}}_{\vartheta }^{{\mathcal{R}}_0}& =\frac{\vartheta }{{\mathcal{R}}_0}\left[\frac{\upsilon }{\tau \left(\chi +{\varrho }_1+\nu +\tau \right)}\right]>0,\hfill \\ \hfill {\mathbb{S}}_{\chi }^{{\mathcal{R}}_0}& =-\frac{\chi }{{\mathcal{R}}_0}\left[\frac{\upsilon \vartheta }{\tau {\left(\chi +{\varrho }_1+\nu +\tau \right)}^2}\right]<0,\hfill \\ \hfill {\mathbb{S}}_{{\varrho }_1}^{{\mathcal{R}}_0}& =-\frac{{\varrho }_1}{{\mathcal{R}}_0}\left[\frac{\upsilon \vartheta }{\tau {\left(\chi +{\varrho }_1+\nu +\tau \right)}^2}\right]<0,\hfill \\ \hfill {\mathbb{S}}_{\upsilon }^{{\mathcal{R}}_0}& =\frac{\upsilon }{{\mathcal{R}}_0}\left[\frac{\vartheta }{\tau {\left(\chi +{\varrho }_1+\nu +\tau \right)}^2}\right]>0.\hfill \end{array}\right\}\end{array}$

Now, we present a Table 3 for the numerical values of the sensitivity analysis such that

The graphical representation is the following:

In Figure 1, the sensitivity indices are graphically presented as a bar chart for the recognition of dominant parameters in which, by increase or decrease, there is a large effect in the model (6) associated to those parameters. In the presented bar chart (Figure 1), the most sensitive parameters are $\vartheta$ and $\upsilon$. Therefore, by increasing the values of these sensitive parameters, the basic reproduction number also increases.

6. Optimal Control Strategies

The COVID-19 model, as presented in Equation (6), is subjected to the optimal control strategy. To reduce coronavirus infection and further propagation in the population, we explore four potential controls. These controls can be categorized as follows: the vaccine is specified as the control variable ${\mathbf{u}}_1$, and the treatment of the infected human population is defined as the control variable ${\mathbf{u}}_2$.

Consider the model (6) with optimal control variables ${\mathbf{u}}_1$ and ${\mathbf{u}}_2$.

(34)$\begin{array}{c}\hfill \left.\begin{array}{cc}\hfill {}_0^{\mathbf{FFP}}{\mathcal{D}}^{\kappa ,\lambda }S\left(t\right)& =\upsilon -({\mathbf{u}}_1+\tau )S\left(t\right)-\frac{\vartheta S\left(t\right)I\left(t\right)}{1+\omega I\left(t\right)}-\theta P\left(t\right),\hfill \\ \hfill {}_0^{\mathbf{FFP}}{\mathcal{D}}^{\kappa ,\lambda }I\left(t\right)& =\frac{\vartheta S\left(t\right)I\left(t\right)}{1+\omega I\left(t\right)}-(\tau +{\mathbf{u}}_2+\nu +{\varrho }_1)I\left(t\right),\hfill \\ \hfill {}_0^{\mathbf{FFP}}{\mathcal{D}}^{\kappa ,\lambda }Q\left(t\right)& =\nu I\left(t\right)-(\tau +\varpi +{\varrho }_2)Q\left(t\right),\hfill \\ \hfill {}_0^{\mathbf{FFP}}{\mathcal{D}}^{\kappa ,\lambda }\mathcal{P}\left(t\right)& ={\mathbf{u}}_{1}S\left(t\right)+{\mathbf{u}}_{2}I\left(t\right)+\varpi Q\left(t\right)-\tau \mathcal{P}\left(t\right)-\theta P\left(t\right),\hfill \\ \hfill {}_0^{\mathbf{FFP}}{\mathcal{D}}^{\kappa ,\lambda }{D}_{IQ}& ={\varrho }_{1}I\left(t\right)+{\varrho }_{2}Q\left(t\right).\hfill \end{array}\right\}\end{array}$

The primary aim is to reduce the objective function, which is defined as a function of the state variables subject to non-negativity constraints and non-negative initial conditions.

(35)$J\left(u_1,u_2\right)={\int }_0^{T_f}\left(z_{1}I+z_{2}Q+\frac{1}{2}\left[B_4{u}_1^2+B_5{u}_2^2\right]\right)dt,$

subject to the control system (34). The parameters governing weight and balance are denoted by the constants $B_j$, where j ranges from 1 to 5, and the terminal time is denoted by $T_f$. Due to the nonlinearity of the control cost associated with intervention, quadratic control functions are employed, indicating that the relationship between the cost of intervention for infected individuals and the impact of intervention is not linear. The aforementioned controls are Lebesgue-integrable functions that are bounded. Our objective is to obtain optimal controls $u_i^{*}$ for i equal to 1 or 2.

(36)$J\left(u_i^{*}\right)=\underset{U}{min}\left(J\left(u_i\right)\right),$

where U is denoted to be the control set given by

(37)$U=\left\{\left(u_i\right):\left[0,T_f\right]\to [0,1],\left(u_i\right),\mathrm{is}\mathrm{Lebesgue}\mathrm{measurable}\right\}.$

Characterization of an Optimal Control

The theory of Pontryagin’s maximum principle [21] provides the necessary requirements that an optimal control must meet. The systems (34) and (35) are transformed into a problem of minimizing a Hamiltonian H pointwise with regard to the controls ${\mathbf{u}}_i,i=1,2$ using this principle. To begin, we must write the Hamiltonian H, which is given by

(38)$\begin{array}{cc}\hfill H=& z_{1}I+z_{2}Q\frac{1}{2}+\left[{\mathsf{\rho }}_1{\mathbf{u}}_1^2+{\rho }_2{\mathbf{u}}_2^2\right]\hfill \\ \hfill & +{\mathsf{\Phi }}_1\left[\upsilon -({\mathbf{u}}_1+\tau )S\left(t\right)-\frac{\vartheta S\left(t\right)I\left(t\right)}{1+\omega I\left(t\right)}-\theta P\left(t\right)\right]\hfill \\ \hfill & +{\mathsf{\Phi }}_2\left[\frac{\vartheta S\left(t\right)I\left(t\right)}{1+\omega I\left(t\right)}-(\tau +{\mathbf{u}}_2+\nu +{\varrho }_1)I\left(t\right)\right]\hfill \\ \hfill & +{\mathsf{\Phi }}_3\left[\nu I\left(t\right)-(\tau +\varpi +{\varrho }_2)Q\left(t\right)\right]\hfill \\ \hfill & +{\mathsf{\Phi }}_4\left[{\mathbf{u}}_{1}S\left(t\right)+{\mathbf{u}}_{2}I\left(t\right)+\varpi Q\left(t\right)-\tau \mathcal{P}\left(t\right)-\theta P\left(t\right)\right]\hfill \\ \hfill & +{\mathsf{\Phi }}_5\left[{\varrho }_{1}I\left(t\right)+{\varrho }_{2}Q\left(t\right)\right],\hfill \end{array}$

where ${\mathsf{\Phi }}_i$ for $i=1,\dots ,5$ denote the adjoint variables associated to the state variables $S,I,Q,P$, and R. For each state variable, the adjoint system can be obtained by taking partial derivatives of the Hamiltonian (38) with respect to the state variables.

Theorem 5.

If the control system (38) is minimized over U by the given optimal controls $u_i^{}$ and solutions $S^{},I^{},Q^{},P^{},D^{}$, then there exist adjoint variables ${\lambda }_i$ that satisfy the equation $\frac{\partial {\mathsf{\Phi }}_i}{dt}=-\frac{\partial H}{\partial i}$, with the transversal conditions given by

(39) ${\mathsf{\Phi }}_i\left(T_f\right)=0,i=S,I,Q,P,D.$

Now, consider the characteristic equations of the control variables, such that

(40) $\left.\begin{array}{cc}\hfill {\mathbf{u}}_1& :={\rho }_1{\mathbf{u}}_1+{\mathsf{\Phi }}_4{S}^{*}-{\mathsf{\Phi }}_1{S}^{*}=0,\hfill \\ \hfill {\mathbf{u}}_2& :={\rho }_2{\mathbf{u}}_2-I^{*}{\mathsf{\Phi }}_2+I^{*}{\mathsf{\Phi }}_4=0.\hfill \end{array}\right\}$

and the optimality condition is given by $\frac{\partial H}{\partial u_i}=0,i=1,\dots ,4.$ Furthermore, we have the controls $u_i^{*}$

(41) $\left.\begin{array}{cc}\hfill & {\mathbf{u}}_1^{*}=min\left\{1,max\left[0,\frac{S^{*}\left({\mathsf{\Phi }}_1-{\mathsf{\Phi }}_4\right)}{{\rho }_1}\right]\right\},\hfill \\ \hfill & {\mathbf{u}}_2^{*}=min\left\{1,max\left[0,\frac{I^{*}\left({\mathsf{\Phi }}_2-{\mathsf{\Phi }}_4\right)}{{\rho }_2}\right]\right\}.\hfill \end{array}\right\}$

Proof. 

The findings of Fleming and Rishel [22] and Sinan et al. [23] provide assurance regarding the presence of an optimal control problem and the corresponding adjoint system, which are derived through the partial derivatives of the Hamiltonian H in Equation (35) and the assessment of the optimal control variables ${\mathbf{u}}_1$ and ${\mathbf{u}}_2$, respectively.

(42)$\left.\begin{array}{cc}\hfill {}_0^{\mathbf{FFP}}{\mathcal{D}}^{\kappa ,\lambda }{\mathsf{\Phi }}_1\left(t\right)& ={\mathsf{\Phi }}_4{\mathbf{u}}_1-{\mathsf{\Phi }}_1\left(\tau +{\mathbf{u}}_1+\frac{I^{*}\vartheta }{I^{*}\omega +1}\right)+\frac{I^{*}{\mathsf{\Phi }}_2\vartheta }{I^{*}\omega +1},\hfill \\ \hfill {}_0^{\mathbf{FFP}}{\mathcal{D}}^{\kappa ,\lambda }{\mathsf{\Phi }}_2\left(t\right)& =z_1+{\mathsf{\Phi }}_5{\varrho }_1+{\mathsf{\Phi }}_3\nu +{\mathsf{\Phi }}_4{\mathbf{u}}_2-{\mathsf{\Phi }}_2\left({\varrho }_1+\nu +\tau +{\mathbf{u}}_2-\frac{S^{*}\vartheta }{I^{*}\omega +1}+\frac{I^{*}{S}^{*}\omega \vartheta }{{\left(I^{*}\omega +1\right)}^2}\right)\hfill \\ \hfill & -\frac{{\mathsf{\Phi }}_1{S}^{*}\vartheta }{{\left(I^{*}\omega +1\right)}^2},\hfill \\ \hfill {}_0^{\mathbf{FFP}}{\mathcal{D}}^{\kappa ,\lambda }{\mathsf{\Phi }}_3\left(t\right)& =z_2+{\mathsf{\Phi }}_5{\varrho }_2+{\mathsf{\Phi }}_4\varpi -{\mathsf{\Phi }}_3\left({\varrho }_2+\tau +\varpi \right),\hfill \\ \hfill {}_0^{\mathbf{FFP}}{\mathcal{D}}^{\kappa ,\lambda }{\mathsf{\Phi }}_4\left(t\right)& ={\mathsf{\Phi }}_1\theta -{\mathsf{\Phi }}_4\left(\tau +\theta \right),\hfill \\ \hfill {}_0^{\mathbf{FFP}}{\mathcal{D}}^{\kappa ,\lambda }{\mathsf{\Phi }}_5\left(t\right)& =0.\hfill \end{array}\right\}$

This completes the proof. For the graphical representation of the optimal control we provide Figure 2 and Figure 3 in which the vaccination and treatment rate show best results in protected population $P\left(t\right)$. □

7. Numerical Scheme for the Fractal-Fractional Model

Consider the fractal-fractional model (6) and apply the scheme such that

(43)$\left.\begin{array}{cc}\hfill {}_0^{\mathbf{FFP}}{\mathcal{D}}^{\kappa ,\lambda }S\left(t\right)& =\upsilon -\tau S\left(t\right)-\frac{\vartheta S\left(t\right)I\left(t\right)}{1+\omega I\left(t\right)},\hfill \\ \hfill {}_0^{\mathbf{FFP}}{\mathcal{D}}^{\kappa ,\lambda }I\left(t\right)& =\frac{\vartheta S\left(t\right)I\left(t\right)}{1+\omega I\left(t\right)}-(\tau +\chi +\nu +{\varrho }_1)I\left(t\right),\hfill \\ \hfill {}_0^{\mathbf{FFP}}{\mathcal{D}}^{\kappa ,\lambda }Q\left(t\right)& =\nu I\left(t\right)-(\tau +\varpi +{\varrho }_2)Q\left(t\right),\hfill \\ \hfill {}_0^{\mathbf{FFP}}{\mathcal{D}}^{\kappa ,\lambda }\mathcal{P}\left(t\right)& =\chi I\left(t\right)+\varpi Q\left(t\right)-\tau \mathcal{P}\left(t\right),\hfill \\ \hfill {}_0^{\mathbf{FFP}}{\mathcal{D}}^{\kappa ,\lambda }{D}_{IQ}& ={\varrho }_{1}I\left(t\right)+{\varrho }_{2}Q\left(t\right).\hfill \end{array}\right\}$

Now, the system of Equation (43) implies that

(44)$\left.\begin{array}{c}\hfill \begin{array}{cc}\hfill & S\left(t\right)=S\left(0\right)+\frac{\kappa }{\mathsf{\Gamma }\left(\lambda \right)}{\int }_0^t{r}^{\kappa -1}{(t-x)}^{\lambda -1}{v}_1\left(S,x\right)dx,\hfill \\ \hfill & I\left(t\right)=I\left(0\right)+\frac{\kappa }{\mathsf{\Gamma }\left(\lambda \right)}{\int }_0^t{r}^{\kappa -1}{(t-x)}^{\lambda -1}{v}_2\left(I,x\right)dx,\hfill \\ \hfill & Q\left(t\right)=Q\left(0\right)+\frac{\kappa }{\mathsf{\Gamma }\left(\lambda \right)}{\int }_0^t{r}^{\kappa -1}{(t-x)}^{\lambda -1}{v}_3\left(Q,x\right)dx,\hfill \\ \hfill & P\left(t\right)=P\left(0\right)+\frac{\kappa }{\mathsf{\Gamma }\left(\lambda \right)}{\int }_0^t{r}^{\kappa -1}{(t-x)}^{\lambda -1}{v}_4\left(P,x\right)dx,\hfill \\ \hfill & D_{IQ}\left(t\right)=D_{IQ}\left(0\right)+\frac{\kappa }{\mathsf{\Gamma }\left(\lambda \right)}{\int }_0^t{r}^{\kappa -1}{(t-x)}^{\lambda -1}{v}_5\left(D_{IQ},x\right)dx.\hfill \end{array}\end{array}\right\}$

The resultant scheme for the system (44) is given as

(45)$\begin{array}{cc}\hfill S(m+1)=& S\left(0\right)+\frac{\kappa {\hslash }^{\lambda }}{\mathsf{\Gamma }(\lambda +2)}\sum _{p=0}^m\left[t_{\left(p\right)}^{\kappa -1}v\left(S_{\left(p\right)},t_p\right)\right.\hfill \\ \hfill & \times \left({(m+1-p)}^{\kappa }(m-p+2+\kappa )-{(m-p)}^{\kappa }(m-p+2+2\kappa )\right)\hfill \\ \hfill & \left.t_{(p-1)}^{\kappa -1}v\left(S_{(p-1)},t_{(p-1)}\right)\left({(m+1-p)}^{\kappa }+1-{(m-p)}^{\kappa }(m-p+1+\kappa )\right)\right],\hfill \end{array}$

(46)$\begin{array}{cc}\hfill I(m+1)=& I\left(0\right)+\frac{\kappa {\hslash }^{\lambda }}{\mathsf{\Gamma }(\lambda +2)}\sum _{p=0}^m\left[t_{\left(p\right)}^{\kappa -1}v\left(I_{\left(p\right)},t_p\right)\right.\hfill \\ \hfill & \times \left({(m+1-p)}^{\kappa }(m-p+2+\kappa )-{(m-p)}^{\kappa }(m-p+2+2\kappa )\right)\hfill \\ \hfill & \left.-t_{(p-1)}^{\kappa -1}v\left(I_{(p-1)},t_{(p-1)}\right)\left({(m+1-p)}^{\kappa }+1-{(m-p)}^{\kappa }(m-p+1+\kappa )\right)\right],\hfill \end{array}$

(47)$\begin{array}{cc}\hfill Q(m+1)=& Q\left(0\right)+\frac{\kappa {\hslash }^{\lambda }}{\mathsf{\Gamma }(\lambda +2)}\sum _{p=0}^m\left[t_{\left(p\right)}^{\kappa -1}v\left(Q_{\left(p\right)},t_p\right)\right.\hfill \\ \hfill & \times \left({(m+1-p)}^{\kappa }(m-p+2+\kappa )-{(m-p)}^{\kappa }(m-p+2+2\kappa )\right)\hfill \\ \hfill & \left.-t_{(p-1)}^{\kappa -1}v\left(Q_{(p-1)},t_{(p-1)}\right)\left({(m+1-p)}^{\kappa }+1-{(m-p)}^{\kappa }(m-p+1+\kappa )\right)\right],\hfill \end{array}$

(48)$\begin{array}{cc}\hfill P(m+1)=& P\left(0\right)+\frac{\kappa {\hslash }^{\lambda }}{\mathsf{\Gamma }(\lambda +2)}\sum _{p=0}^m\left[t_{\left(p\right)}^{\kappa -1}v\left(P_{\left(p\right)},t_p\right)\right.\hfill \\ \hfill & \times \left({(m+1-p)}^{\kappa }(m-p+2+\kappa )-{(m-p)}^{\kappa }(m-p+2+2\kappa )\right)\hfill \\ \hfill & \left.-t_{(p-1)}^{\kappa -1}v\left(P_{(p-1)},t_{(p-1)}\right)\left({(m+1-p)}^{\kappa }+1-{(m-p)}^{\kappa }(m-p+1+\kappa )\right)\right],\hfill \end{array}$

(49)$\begin{array}{cc}\hfill D_{IQ}(m+1)=& D_{IQ}\left(0\right)+\frac{\kappa {\hslash }^{\lambda }}{\mathsf{\Gamma }(\lambda +2)}\sum _{p=0}^m\left[t_{\left(p\right)}^{\kappa -1}v\left(D_{IQ\left(p\right)},t_p\right)\right.\hfill \\ \hfill & \times \left({(m+1-p)}^{\kappa }(m-p+2+\kappa )-{(m-p)}^{\kappa }(m-p+2+2\kappa )\right)\hfill \\ \hfill & \left.-t_{(p-1)}^{\kappa -1}v\left(D_{IQ(p-1)},t_{(p-1)}\right)\left({(m+1-p)}^{\kappa }+1-{(m-p)}^{\kappa }(m-p+1+\kappa )\right)\right].\hfill \end{array}$

8. Results and Discussion

We have simulated our model for various values of fractal dimensions and different fractional orders, and the results are presented in Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13, Figure 14, Figure 15, Figure 16, Figure 17, Figure 18, Figure 19, Figure 20, Figure 21, Figure 22, Figure 23, Figure 24, Figure 25, Figure 26, Figure 27, Figure 28, Figure 29, Figure 30 and Figure 31. To protect the human population from reinfection and treat infected humans in hospitals, we implemented an optimal control strategy using Pontryagin’s maximum principle, as used by [24,25,26]. Two optimal control variables, ${\mathbf{u}}_1$ and ${\mathbf{u}}_2$, were employed to represent vaccination of the susceptible class and treatment of the infected class, respectively. The resulting figure, shown in Figure 2, demonstrates an increase in the protected human population for fractal-fractional orders with $\kappa =0.5$ and $\lambda =0.5$. The control variables are presented in Figure 3 and Figure 32. A sensitivity analysis was conducted on each parameter’s basic reproduction number, and the findings were displayed in a bar chart in Figure 1. Additionally, Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8 present three-dimensional plots of each parameter. The results of the sensitivity analysis indicate that the parameters $\tau$, $\vartheta$, $\chi$, and $\upsilon$ are more sensitive than the other parameters. The behaviour of surfaces in Figure 33, Figure 34, Figure 35, Figure 36 and Figure 37 show the effect in $R_0$ due to the variation in parameters $\upsilon$ and $\vartheta$, $\tau$ and $\vartheta$, $\chi$ and $\vartheta$, $\nu$ and $\vartheta$ while ${\varrho }_1$ and $\vartheta$. By examining the sensitivity of certain parameters, a meaningful connection can be established among the uncertain parameters in the model. From Figure 19, Figure 20, Figure 21 and Figure 22, it reflects that the COVID-19 pandemic has had a significant impact on the mental health and well-being of people worldwide. The uncertainty, fear, and stress caused by the pandemic have led to a range of psychological effects, including: Anxiety and Depression: The pandemic has caused a lot of anxiety and depression among people due to the fear of getting infected, losing loved ones, and the uncertainty of the future. Social Isolation: The pandemic has led to social distancing measures, which have resulted in social isolation and loneliness, leading to depression and anxiety. Post-Traumatic Stress Disorder (PTSD): People who have been infected with COVID-19 or have lost loved ones to the virus may develop PTSD, which can cause flashbacks, nightmares, and anxiety. Substance Abuse: The pandemic has led to an increase in substance abuse, including alcohol and drug abuse, as people try to cope with the stress and anxiety caused by the pandemic. Burnout: Healthcare workers and essential workers have been working tirelessly during the pandemic, leading to burnout and exhaustion. Overall, the COVID-19 pandemic has had a significant impact on the mental health and well-being of people worldwide. So, the factors described above can cause deaths.

9. Conclusions

In this manuscript, we have established a comprehensive analysis for a fractal-fractional-type model addressing the COVID-19 disease. We have computed the fundamental threshold number, which allows us to predict the transmission of the disease in the community. Additionally, we have carried out a sensitivity analysis of some parameters in the considered model and developed a bifurcation analysis for the proposed model. Bifurcation analysis is a powerful procedure to investigate the steady-state nonlinear dynamics of systems, and we have used the methodology adopted by Castillo-Chavez and Song. We have also established an optimal control strategy using Pontryagin’s maximum principle in the proposed model. Finally, we have simulated the proposed model by considering various values of fractal and fractional order corresponding to some real values. The sensitivity and bifurcation analysis of the dynamical system enables us to investigate the abrupt qualitative changes in the considered system, and the fractal-fractional order derivative allows for a more precise explanation of the dynamics. In the future, the analysis established in this paper can be extended to a more complex biological model. The coronavirus pandemic has had a significant impact on the mental health and well-being of people worldwide. The fear and uncertainty surrounding the virus, coupled with the social isolation and economic stress caused by lockdowns and restrictions, have led to increased levels of anxiety, depression, and other mental health issues. People have also experienced grief and loss due to the death of loved ones, as well as the loss of jobs, businesses, and social connections. The pandemic has also highlighted existing inequalities and disparities in access to healthcare and resources, which can further exacerbate mental health issues. However, it is important to note that people have also shown resilience and adaptability in the face of the pandemic. Many have found ways to cope with the challenges and have even experienced personal growth and positive changes in their lives. Seeking support from mental health professionals, staying connected with loved ones, and practicing self-care can all help mitigate the psychological effects of the pandemic.

Author Contributions

Conceptualization, N.H.A.; Writing—original draft, M.S. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

All the data are within the manuscript.

Acknowledgments

This research was supported by the Deanship of Scientific Research, Imam Mohammad Ibn Saud Islamic University (IMSIO), Saudi Arabia, Grant No. (21-13-18-043).

Conflicts of Interest

The authors have no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Plot of Sensitivity Analysis based on [Forumla omitted. See PDF.] with each parameter associated with [Forumla omitted. See PDF.].

Figure 2. The implementation of fractal-fractional orders with [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] can lead to optimal control of COVID-19. This is due to the influence of vaccination and treatment factors, represented by [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], respectively.

Figure 3. The dynamics of vaccine control variable [Forumla omitted. See PDF.] for [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Figure 4. The kinetics of the human population that is susceptible to a given condition is being studied for a fractal order of [Forumla omitted. See PDF.] and different fractional orders denoted by [Forumla omitted. See PDF.].

Figure 5. Dynamics of Infected Human Population for Fractal Order [Forumla omitted. See PDF.] and various Fractional Orders, [Forumla omitted. See PDF.].

Figure 6. Dynamics of Quarantine Human Population for Fractal Order [Forumla omitted. See PDF.] and various Fractional Orders, [Forumla omitted. See PDF.].

Figure 7. Dynamics of Protected Human Population for Fractal Order [Forumla omitted. See PDF.] and various Fractional Orders, [Forumla omitted. See PDF.].

Figure 8. Dynamics of Deaths in Quarantine and Infected Human Populations for Fractal Order [Forumla omitted. See PDF.] and various Fractional Orders, [Forumla omitted. See PDF.].

Figure 9. Dynamics of Susceptible Human Populations for Fractional Order [Forumla omitted. See PDF.] and various Fractal Orders, [Forumla omitted. See PDF.].

Figure 10. Dynamics of Infected Human Populations for Fractional Order [Forumla omitted. See PDF.] and various Fractal Orders, [Forumla omitted. See PDF.].

Figure 11. Dynamics of Quarantine Human Populations for Fractional Order [Forumla omitted. See PDF.] and various Fractal Orders, [Forumla omitted. See PDF.].

Figure 12. Dynamics of Protected Human Populations for Fractional Order [Forumla omitted. See PDF.] and various Fractal Orders, [Forumla omitted. See PDF.].

Figure 13. Dynamics of Deaths in Quarantine and Infected Human Populations for Fractional Order [Forumla omitted. See PDF.] and various Fractal Orders, [Forumla omitted. See PDF.].

Figure 14. Dynamics of Susceptible Human Population for various Fractal-Fractional Orders.

Figure 15. Dynamics of Infected Human Population for various Fractal-Fractional Orders.

Figure 16. Dynamics of Quarantine Human Population for various Fractal-Fractional Orders.

Figure 17. The study investigates the behavior of safeguarded human communities under different fractal-fractional orders.

Figure 18. Dynamics of Deaths in Quarantine and Infected Human Populations for various Fractal-Fractional Orders.

Figure 19. The study investigates the impact of the psychological effect parameter, [Forumla omitted. See PDF.], on the dynamics of the susceptible human population.

Figure 20. This research examines the influence of the psychological factor parameter, [Forumla omitted. See PDF.], on the kinetics of the infected human population.

Figure 21. The study investigates the behavior of protected population under different values of the psychological impact parameter, denoted as [Forumla omitted. See PDF.].

Figure 22. Dynamics of Deaths in Quarantine and Infected Human Populations for various values of Parameter of Psychological Effect, [Forumla omitted. See PDF.].

Figure 23. Phase Portraits of Susceptible Human vs. Infected Human Populations with Fractal–Fractional orders [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.] and with Initial condition larger than [Forumla omitted. See PDF.].

Figure 24. Phase Portraits of Susceptible Human vs. Quarantine Human Populations with Fractal-Fractional orders [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.] and with Initial condition larger than [Forumla omitted. See PDF.].

Figure 25. Phase Portraits of Susceptible Human vs. Protected Human Populations with Fractal-Fractional orders [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.] and with Initial condition larger than [Forumla omitted. See PDF.].

Figure 26. Phase Portraits of Susceptible Human vs. [Forumla omitted. See PDF.] Human Populations with Fractal-Fractional orders [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.] and with Initial condition larger than [Forumla omitted. See PDF.].

Figure 27. In this figure, we check for the stability of the solution of Susceptible Human Populations with Fractal-Fractional orders [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.] in which we experience the fluctuation with different initial conditions starting from [Forumla omitted. See PDF.] and above. Light blue color nearly, green at [Forumla omitted. See PDF.], indigo color at [Forumla omitted. See PDF.], yellow at [Forumla omitted. See PDF.], red at [Forumla omitted. See PDF.], and dark blue at [Forumla omitted. See PDF.].

Figure 28. Stability Curves of Infected Human Populations with Fractal-Fractional orders [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.] and with with various Initial conditions in which the Initial conditions larger than [Forumla omitted. See PDF.] starts fluctuation.

Figure 29. Stability Curves of Quarantine Human Populations with Fractal-Fractional orders [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.] and with with various Initial conditions.

Figure 30. Stability Curves of Protected Human Populations with Fractal-Fractional orders [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.] and with with various Initial conditions.

Figure 31. Stability Curves of [Forumla omitted. See PDF.] Human Populations with Fractal-Fractional orders [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.] and with with various Initial conditions.

Figure 32. The dynamics of treatment control variable [Forumla omitted. See PDF.] for [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Figure 33. The 3D dynamics of [Forumla omitted. See PDF.] with parameters [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] associated with its contour plot.

Figure 34. The 3D dynamics of [Forumla omitted. See PDF.] with parameters [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] associated with its contour plot.

Figure 35. The 3D dynamics of [Forumla omitted. See PDF.] with parameters [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] associated with its contour plot.

Figure 36. The 3D dynamics of [Forumla omitted. See PDF.] with parameters [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] associated with its contour plot.

Figure 37. The 3D dynamics of [Forumla omitted. See PDF.] with parameters [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] associated with its contour plot.

Table 1

A tabular representation of the demographic subgroups with their respective Initial coditions.

Symbol	Description of Compartment/Class	Initial Conditions
$S\left(t\right)$	Susceptible Human Population	0.85
$I\left(t\right)$	Exposed Human Population	3
$Q\left(t\right)$	Infected Human Population	1.9
$P\left(t\right)$	Recovered Human Population	300
$D_{IQ}\left(t\right)$	Susceptible Vector Population	100

Table 2

A tabular representation of the value of Parameters and description.

Symbol	Description of Parameter	Value
$\upsilon$	Recruitment Rate [13].	$\tau \times N\left(0\right)$
$\tau$	Natural Death Rate [13].	$\frac{1}{67.7\times 365}$
$\vartheta$	Transmission rate.	0.2784
$\omega$	Psychological Effect on Humans [14].	[0,1]
$\chi$	Recovery Rate of Infected Population [15].	0.1
$\varpi$	Recovery Rate of Quarantine Population.	0.020
$\nu$	Incubation Period.	0.010
${\varrho }_1$	The mortality rate of the afflicted populace as a result of a pathological condition.	0.015
${\varrho }_2$	Mortality rate among individuals subjected to quarantine as a result of disease.	0.015

Table 3

Sensitivity of the $R_0$ versus proposed parameters.

Parameter	Sensitivity Index	Value	Parameter	Sensitivity Index	Value
$\nu$	${\mathbf{s}}_{\nu }^{{\mathcal{R}}_0}$	−0.0799	$\tau$	${\mathbf{s}}_{\tau }^{{\mathcal{R}}_0}$	−1.0003
$\vartheta$	${\mathbf{s}}_{\vartheta }^{{\mathcal{R}}_0}$	1	$\chi$	${\mathbf{s}}_{\chi }^{{\mathcal{R}}_0}$	−0.7997
${\varrho }_1$	${\mathbf{s}}_{{\varrho }_1}^{{\mathcal{R}}_0}$	−0.1199	$\upsilon$	${\mathbf{s}}_{\upsilon }^{{\mathcal{R}}_0}$	1

References

1. Anderson, P.K.; Cunningham, A.A.; Patel, N.G.; Morales, F.J.; Epstein, P.R.; Daszak, P. Emerging infectious diseases of plants: Pathogen pollution, climate change and agrotechnology drivers. Trends Ecol. Evol.; 2004; 19, pp. 535-544. [DOI: https://dx.doi.org/10.1016/j.tree.2004.07.021] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/16701319]

2. Omame, A.; Isah, M.E.; Abbas, M. An optimal control model for COVID-19, zika, dengue, and chikungunya co-dynamics with reinfection. Optim. Control Appl. Methods; 2023; 44, pp. 170-204. [DOI: https://dx.doi.org/10.1002/oca.2936] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/36248678]

3. Saha, P.; Biswas, S.K.; Biswas, M.H.A.; Ghosh, U. An SEQAIHR model to study COVID-19 transmission and optimal control strategies in Hong Kong, 2022. Nonlinear Dyn.; 2023; 111, pp. 6873-6893. [DOI: https://dx.doi.org/10.1007/s11071-022-08181-0] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/36644569]

4. Anderson, R.M.; May, R.M. Infectious Diseases of Humans: Dynamics and Control; Oxford University Press: Oxford, UK, 1992.

5. Khan, H.; Ahmad, F.; Tunç, O.; Idrees, M. On fractal-fractional COVID-19 mathematical model. Chaos Solitons Fractals; 2022; 157, 111937. [DOI: https://dx.doi.org/10.1016/j.chaos.2022.111937]

6. Hamou, A.A.; Azroul, E.; Hammouch, Z.; Alaoui, A.L. A fractional multi-order model to predict the COVID-19 outbreak in Morocco. Appl. Comput. Math.; 2020; 20, pp. 177-203.

7. Lu, Z.; Yu, Y.; Chen, Y.; Ren, G.; Xu, C.; Wang, S.; Yin, Z. A fractional-order SEIHDR model for COVID-19 with inter-city networked coupling effects. Nonlinear Dyn.; 2020; 101, pp. 1717-1730. [DOI: https://dx.doi.org/10.1007/s11071-020-05848-4]

8. Atangana, A.; Qureshi, S. Modeling attractors of chaotic dynamical systems with fractal–fractional operators. Chaos Solitons Fractals; 2019; 123, pp. 320-337. [DOI: https://dx.doi.org/10.1016/j.chaos.2019.04.020]

9. Atangana, A. Fractal-fractional differentiation and integration: Connecting fractal calculus and fractional calculus to predict complex system. Chaos Solitons Fractals; 2017; 102, pp. 396-406. [DOI: https://dx.doi.org/10.1016/j.chaos.2017.04.027]

10. Qureshi, S.; Atangana, A. Mathematical analysis of dengue fever outbreak by novel fractional operators with field data. Phys. Stat. Mech. Its Appl.; 2019; 526, 121127. [DOI: https://dx.doi.org/10.1016/j.physa.2019.121127]

11. Li, Z.; Liu, Z.; Khan, M.A. Fractional investigation of bank data with fractal-fractional Caputo derivative. Chaos Solitons Fractals; 2020; 131, 109528. [DOI: https://dx.doi.org/10.1016/j.chaos.2019.109528]

12. Zhang, L.; ur Rahman, M.; Haidong, Q.; Arfan, M.; Adnan,. Fractal-fractional Anthroponotic Cutaneous Leishmania model study in sense of Caputo derivative. Alex. Eng. J.; 2021; 61, pp. 4423-4433. [DOI: https://dx.doi.org/10.1016/j.aej.2021.10.001]

13. Shen, Z.H.; Chu, Y.M.; Khan, M.A.; Muhammad, S.; Al-Hartomy, O.A.; Higazy, M. Mathematical modeling and optimal control of the COVID-19 dynamics. Results Phys.; 2021; 31, 105028. [DOI: https://dx.doi.org/10.1016/j.rinp.2021.105028]

14. Liu, S.; Pang, L.; Ruan, S.; Zhang, X. Global dynamics of avian influenza epidemic models with psychological effect. Comput. Math. Methods Med.; 2015; 2015, 913726. [DOI: https://dx.doi.org/10.1155/2015/913726]

15. Gu, Y.; Khan, M.A.; Hamed, Y.; Felemban, B.F. A Comprehensive Mathematical Model for SARS-CoV-2 in Caputo Derivative. Fractal Fract.; 2021; 5, 271. [DOI: https://dx.doi.org/10.3390/fractalfract5040271]

16. Lu, H.; Ding, Y.; Gong, S.; Wang, S. Mathematical modeling and dynamic analysis of SIQR model with delay for pandemic COVID-19. Math. Biosci. Eng.; 2021; 18, pp. 3197-3214. [DOI: https://dx.doi.org/10.3934/mbe.2021159]

17. Bozdağ, F. The psychological effects of staying home due to the COVID-19 pandemic. J. Gen. Psychol.; 2021; 148, pp. 226-248. [DOI: https://dx.doi.org/10.1080/00221309.2020.1867494]

18. Castillo-Chavez, C.; Song, B. Dynamical models of tuberculosis and their applications. Math. Biosci. Eng.; 2004; 1, 361. [DOI: https://dx.doi.org/10.3934/mbe.2004.1.361]

19. Buonomo, B.; Lacitignola, D. On the backward bifurcation of a vaccination model with nonlinear incidence. Nonlinear Anal. Model. Control; 2011; 16, pp. 30-46. [DOI: https://dx.doi.org/10.15388/NA.16.1.14113]

20. Chitnis, N.; Hyman, J.M.; Cushing, J.M. Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model. Bull. Math. Biol.; 2008; 70, pp. 1272-1296. [DOI: https://dx.doi.org/10.1007/s11538-008-9299-0]

21. Pontryagin, L.S. Mathematical Theory of Optimal Processes; CRC Press: Boca Raton, FL, USA, 1987.

22. Fleming, W.H.; Rishel, R.W. Deterministic and Stochastic Optimal Control; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2012; Volume 1.

23. Sinan, M.; Shah, K.; Kumam, P.; Mahariq, I.; Ansari, K.J.; Ahmad, Z.; Shah, Z. Fractional order mathematical modeling of typhoid fever disease. Results Phys.; 2022; 32, 105044. [DOI: https://dx.doi.org/10.1016/j.rinp.2021.105044]

24. Mahmudov, E.N.; Mardanov, M.J. On duality in optimal control problems with second-order differential inclusions and initial-point constraints. Proc. Inst. Math. Mech. Nat. Acad. Sci. Azerb.; 2020; 46, pp. 115-128.

25. Alekseev, V.M.; Tikhomirov, V.M.; Fomin, S.V. Optimal Control; Springer: New York, NY, USA, 2013.

26. Shalsh, O.K.; Al-Nassir, S. Dynamics and optimal Harvesting strategy for biological models with Beverton-Holt growth. Iraqi. J. Sci.; 2020; pp. 223-232. [DOI: https://dx.doi.org/10.24996/ijs.2020.SI.1.30]