# Optimal control for a time delay multi-strain tuberculosis fractional model: a numerical approach

Optimal control for a time delay multi-strain tuberculosis fractional model: a numerical approach Abstract In this article, optimal control for a novel fractional multi-strain Tuberculosis model with time delay is presented. The proposed model is governed by a system of fractional delay differential equations, where the fractional derivative is defined in the Caputo sense. Modified parameters are introduced to account for the fractional order. Two delays in the proposed model representing the time delay on the diagnosis and commencement of treatment of individuals with active Tuberculosis infection in two and three strains. Necessary and sufficient conditions that guarantee the existence and the uniqueness of the solution of the resulting systems are given. Two control variables are proposed to minimize the cost of interventions. Two simple numerical methods are used to study the nonlinear fractional delay optimal control problem. The methods are the iterative optimal control method (IOCM) and the generalized Euler method (GEM). Comparative studies are implemented, it is found that the IOCM is better than the GEM. 1. Introduction Mathematical modelling with delay differential equations (DDEs) is widely used for analysis and predictions in epidemiology, immunology and physiology (Bocharov & Rihan, 2000; Rihan, 2000; Fowler & Mackey, 2002; Nelson & Perelson, 2002; Smith, 2011). The time delays in these models can be considered as dependence of the present state of the modelled system on its history. The general theory of differential equations with delays is discussed in the literature (Cooke & Bellman, 1963; Driver, 1977; Hale, 1977; Van den Driessche & Watmough, 2002; Bellen & Zennaro, 2003; Pimenov & Tashirova, 2013). In reality, some issues are rarely so instantaneous; but there is usually a propagation delay before the effects are felt. Recently, many optimal control models pertaining to epidemic diseases have appeared in the literature. They include, but not limited to, delayed SIRS epidemic model (Laarabi et al., 2015, where the authors proposed an optimal control model governed by a delayed SIRS model with a single time delay representing the incubation period), delayed SIR model (Abta et al., 2014, where the authors formulated an optimal control model to optimize the costs of vaccination and treatment on a delayed SIR model with a single delay representing the incubation period and a saturated incidence rate) and tuberculosis model (Silva et al., 2017), where the authors proposed optimal control of a Tuberculosis model with state and control delays using data from Angola). Sweilam & AL-Mekhlafi (2016b), studied optimal control of a general nonlinear multi-strain tuberculosis (TB) model. Four controls are used for preventing the failure of treatment in active TB. In Sweilam & AL-Mekhlafi (2017), we have also performed numerical studies of optimal control of HIV: HIV model (Hattaf & Yousfi, 2012, where the authors used optimal controls to represent the efficiency of drug treatment in inhibiting viral production and preventing new infections), swine flu (Boklund et al., 2009, where the authors studied epidemiological and economic consequences of some control strategies in a classical swine flu epidemic under Danish conditions with respect to herd demographics and geography and investigated the effect of extra biosecurity measures on farms and dengue fever (Aldila et al., 2013, where the authors designed an optimal control problem using four different control parameters and discussed some results about epidemic prevention). Tuberculosis (TB) can be considered as one of the most important infectious diseases, it is the second largest cause of mortality by infectious diseases and is a challenging disease to control (World Health Organization, 2014). Delays to the diagnosis of active TB present a major obstacle to the control of a TB epidemic (Uys et al., 2007); it may worsen the disease, increase the risk of death and enhance tuberculosis transmission to the community (Toman, 1979; Sreeramareddy et al., 2009). Both patient and the health system may be responsible for the diagnosis delay (Sreeramareddy et al., 2009). On the other hand, mathematical models are quite important and efficient tool to describe and investigate TB epidemics, for more details see (Styblo, 1978; Castillo-ćhavez & Feng, 1997; Cohen & Murray, 2004; Small & Fujiwara, 2001; Denysiuk et al., 2014). In this work, we consider a general model of multi-strain TB diseases with time-delay memory. We introduce a discrete time-delay in the state variable of active TB infection of two and three strain, that represents the delay in the diagnosis and commencement of treatment. The multi-strain TB model which incorporates three strains: drug-sensitive, emerging multi-drug resistant (MDR) and extensively drug-resistant (XDR) is developed by Arino & Soliman (2015). This model included several factors of spreading TB such as the fast infection, the exogenous reinfection and secondary infection along with the resistance factor. Sweilam and AL$$-$$Mekhlafi introduced some numerical studies of this model in (Sweilam & AL-Mekhlafi, 2016,a,b,c; Sweilam et al., 2017). Fractional differential equations have been the focus of many studies due to their frequent appearance in various sciences (Sweilam et al., 2011, 2012a,b; Nagy & Sweilam, 2014; Salahshour et al., 2015; Tariboon et al., 2015; Agarwal & Choi, 2016; Rihan et al., 2016; Sweilam & AL-Mekhlafi, 2016,a,b, 2017; Zhang et al., 2016; Sweilam et al., 2017; Zhou et al., 2017). Delayed fractional differential equations are used to describe dynamical systems (Allen, 2007; Rihan et al., 2016). Recently, DFDEs begin to raises the attention of many researchers (Bhalekar & Daftardar-Gejji, 2011; Wang et al., 2013; Rihan et al., 2015). Simulating these equations are an important technique in the research, accordingly, finding effective numerical methods for the DFDEs are necessary process. In this article, we introduce the fractional order multi-strain TB with memory time-delay and modified parameters, this modified parameters are introduced to account for the fractional order (Arenas et al., 2016). The aim of this article is to study optimal control of the proposed model for minimization of the number of active TB infectious of two and three strains. Two simple numerical methods are used to study the nonlinear fractional delay optimal control problem (FDOCP). The methods are the iterative optimal control method (IOCM) and the generalized Euler method (GEM). Comparative studies are implemented. This article is organized as follows: In Section 2, mathematical preliminaries of fractional calculus which are required for establishing the results are given. In Section 3, the proposed system is described by fractional order derivatives with memory time-delay is presented and the stability of equilibrium points is studied. Formulation of the optimal control problem and the necessary optimality conditions for the fractional order multi-strain TB model are derived in Section 4. Numerical methods are introduced for solving the proposed control model in Section 5. Numerical simulations are given in Section 6. In Section 6, the conclusions are given. 2. Definitions and preliminaries Let us consider the signal fractional differential equation (see Podlubny, 1999; Nagy & Sweilam, 2014), $$D_{t}^{\alpha}z(t)=g(t,z(t)),\,\,\,T \geq t \geq 0,\,\, and,\,z(t_{0})=0. \label{ss}$$ (1) Definition 2.1 Podlubny (1999) The Grünwald–Letnikovs fractional derivative is defined as: \begin{align} D^{\alpha}_{t}g(t)&=\lim_{h\rightarrow0}\frac{1}{h^{\alpha}}\sum_{i=0}^{[\frac{t}{h}]}q^{\alpha}_{i}f(t-ih),\notag \end{align} where $$[\frac{t}{h}]$$ denotes the integer part of $$\frac{t}{h}$$ and \begin{align} q^{\alpha}_{i}&=(-1)^{i}(^{\alpha}_{i})=\frac{{\it{\Gamma}}(\alpha+1)}{i!{\it{\Gamma}}(\alpha-i+1)},\notag \end{align} with $$(^{\alpha}_{i})$$ being the fractional binomial coefficients. It is clear that the coefficients $$q^{\alpha}_{i}$$ can be evaluated recursively as follows : $$q^{\alpha}_{0}=1,$$$$q^{\alpha}_{i}=(1-\frac{\alpha+1}{i})q^{\alpha}_{i-1},$$$$i\geq 1.$$ Definition 2.2 Podlubny (1999) The left Caputo fractional derivative is defined as follows: \begin{align} ^{c}_{a}D^{\alpha}_{t}g(t)&=\frac{1}{{\it{\Gamma}}(n-\alpha)}\int^{t}_{a}(t-\tau)^{n-\alpha-1}\frac{d^{n}}{d\tau^{n}}z(\tau)d\tau,\notag \end{align} the right Caputo fractional derivative is defined as follows : \begin{align} ^{c}_{a}D^{\alpha}_{t}g(t)&=\frac{(-1)^{n}}{{\it{\Gamma}}(n-\alpha)}\int^{b}_{t}(t-\tau)^{n-\alpha-1}\frac{d^{n}}{d\tau^{n}}z(\tau)d\tau.\notag \end{align} Where $${\it{\Gamma}}$$ is the Euler gamma function. 3. Mathematical model In the following, the control problem of the fractional model of multi-strain TB with time-delay memory and modified parameters are presented. This model incorporates three strains: drug-sensitive, MDR and XDR. The population is divided into eight compartments depending on their epidemiological stages as follows: susceptible $$(S)$$; latently infected with drug sensitive TB $$(L_{s})$$; latently infected with MDR-TB $$(L_{m})$$; latently infected with XDR-TB $$(L_{x})$$; sensitive drug TB infectious $$(I_{s})$$; MDR-TB infectious $$(I_{m})$$; XDR-TB infectious $$(I_{x})$$; recovered $$R$$. The new parameters of the model are described in Table 1, this modified parameters are introduced to account for the fractional order (Arenas et al., 2016). We introduce a discrete time-delay in the state variables $$I_{m}$$ and $$I_{x},$$ denoted by $$\tau_{1}$$ and $$\tau_{2}$$, these represents the delay to diagnosis and commencement of treatment of active TB infection of two and three strain. One of the main assumptions of this model is that, the total population, $$N;$$$$N = S(t)+L_{s}(t)+L_{m}(t)+L_{x}(t)+I_{s}(t)+I_{m}(t)+I_{x}(t)+R(t),$$ is a constant in time. In other words, we assume that the birth and death rates are equal and there are no TB-related deaths i.e., $$\delta_{s}^{\alpha}=\delta_{m}^{\alpha}=\delta_{x}^{\alpha}$$=0. Table 1 All adapted parameters and their interpretation of the system (2)–(9) Parameters Interpretation $$d^{\alpha}$$ Death and birth rate $$\lambda_{j}^{\alpha}$$ Rate of infected individuals move to $$L_{j}$$ with strain $$j \in\{s ,m,x\}\$$ $$1-\lambda_{j}^{\alpha}$$ Rate of newly infected individuals progressing to active TB with strain $$j \in\{s,m,x \}$$ $$\beta_{j}^{\alpha}$$ Transmission coefficient with strain $$j \in\{s,m,x \}$$ $$\varepsilon_{j}^{\alpha}$$ Rate of endogenous reactivation of $$L_{j}$$ $$\gamma_{j}^{\alpha}$$ Rate of natural recovery to the latent stage $$L_{j}$$ $$\delta_{j}^{\alpha}$$ Rate of death due to $$TB$$ of strain $$j$$ $$\alpha_{j1}^{\alpha},\alpha_{j2}^{\alpha}$$ Rate of exogenous reinfection of $$L_{j1}$$ due to contact with $$I_{j2}$$ $$1-\sigma_{j}^{\alpha}$$ Efficiency of treatment in preventing infection with strain $$j$$ $$P_{1}^{\alpha}$$ Probability of treatment success for $$L_{s}$$ $$1-P_{1}^{\alpha}$$ Proportion of treated $$L_{s}$$ moved to $$L_{m}$$ due to incomplete treatment or lack of strict compliance in the use of drugs $$P_{2}^{\alpha}$$ Probability of treatment success for $$I_{s}$$ $$1-P_{2}^{\alpha}$$ Proportion of treated $$I_{s}$$ moved to $$L_{m}$$ due to incomplete treatment or lack of strict compliance in the use of drugs $$P_{3}^{\alpha}$$ Probability of treatment success for $$I_{m}$$ $$1-P_{3}^{\alpha}$$ proportion of treated $$I_{m}$$ moved to $$L_{x}$$due to incomplete treatment or lack of strict compliance in the use of drugs $$t_{1s}^{\alpha}$$ Rate of treatment for $$L_{s}$$ $$t_{2j}^{\alpha}$$ Rate of treatment for $$I_{j}$$. Note that $$t_{2x}$$ is the rate of successful treatment of $$I_{x}, j\in\{x,m,s\}$$ Parameters Interpretation $$d^{\alpha}$$ Death and birth rate $$\lambda_{j}^{\alpha}$$ Rate of infected individuals move to $$L_{j}$$ with strain $$j \in\{s ,m,x\}\$$ $$1-\lambda_{j}^{\alpha}$$ Rate of newly infected individuals progressing to active TB with strain $$j \in\{s,m,x \}$$ $$\beta_{j}^{\alpha}$$ Transmission coefficient with strain $$j \in\{s,m,x \}$$ $$\varepsilon_{j}^{\alpha}$$ Rate of endogenous reactivation of $$L_{j}$$ $$\gamma_{j}^{\alpha}$$ Rate of natural recovery to the latent stage $$L_{j}$$ $$\delta_{j}^{\alpha}$$ Rate of death due to $$TB$$ of strain $$j$$ $$\alpha_{j1}^{\alpha},\alpha_{j2}^{\alpha}$$ Rate of exogenous reinfection of $$L_{j1}$$ due to contact with $$I_{j2}$$ $$1-\sigma_{j}^{\alpha}$$ Efficiency of treatment in preventing infection with strain $$j$$ $$P_{1}^{\alpha}$$ Probability of treatment success for $$L_{s}$$ $$1-P_{1}^{\alpha}$$ Proportion of treated $$L_{s}$$ moved to $$L_{m}$$ due to incomplete treatment or lack of strict compliance in the use of drugs $$P_{2}^{\alpha}$$ Probability of treatment success for $$I_{s}$$ $$1-P_{2}^{\alpha}$$ Proportion of treated $$I_{s}$$ moved to $$L_{m}$$ due to incomplete treatment or lack of strict compliance in the use of drugs $$P_{3}^{\alpha}$$ Probability of treatment success for $$I_{m}$$ $$1-P_{3}^{\alpha}$$ proportion of treated $$I_{m}$$ moved to $$L_{x}$$due to incomplete treatment or lack of strict compliance in the use of drugs $$t_{1s}^{\alpha}$$ Rate of treatment for $$L_{s}$$ $$t_{2j}^{\alpha}$$ Rate of treatment for $$I_{j}$$. Note that $$t_{2x}$$ is the rate of successful treatment of $$I_{x}, j\in\{x,m,s\}$$ Table 1 All adapted parameters and their interpretation of the system (2)–(9) Parameters Interpretation $$d^{\alpha}$$ Death and birth rate $$\lambda_{j}^{\alpha}$$ Rate of infected individuals move to $$L_{j}$$ with strain $$j \in\{s ,m,x\}\$$ $$1-\lambda_{j}^{\alpha}$$ Rate of newly infected individuals progressing to active TB with strain $$j \in\{s,m,x \}$$ $$\beta_{j}^{\alpha}$$ Transmission coefficient with strain $$j \in\{s,m,x \}$$ $$\varepsilon_{j}^{\alpha}$$ Rate of endogenous reactivation of $$L_{j}$$ $$\gamma_{j}^{\alpha}$$ Rate of natural recovery to the latent stage $$L_{j}$$ $$\delta_{j}^{\alpha}$$ Rate of death due to $$TB$$ of strain $$j$$ $$\alpha_{j1}^{\alpha},\alpha_{j2}^{\alpha}$$ Rate of exogenous reinfection of $$L_{j1}$$ due to contact with $$I_{j2}$$ $$1-\sigma_{j}^{\alpha}$$ Efficiency of treatment in preventing infection with strain $$j$$ $$P_{1}^{\alpha}$$ Probability of treatment success for $$L_{s}$$ $$1-P_{1}^{\alpha}$$ Proportion of treated $$L_{s}$$ moved to $$L_{m}$$ due to incomplete treatment or lack of strict compliance in the use of drugs $$P_{2}^{\alpha}$$ Probability of treatment success for $$I_{s}$$ $$1-P_{2}^{\alpha}$$ Proportion of treated $$I_{s}$$ moved to $$L_{m}$$ due to incomplete treatment or lack of strict compliance in the use of drugs $$P_{3}^{\alpha}$$ Probability of treatment success for $$I_{m}$$ $$1-P_{3}^{\alpha}$$ proportion of treated $$I_{m}$$ moved to $$L_{x}$$due to incomplete treatment or lack of strict compliance in the use of drugs $$t_{1s}^{\alpha}$$ Rate of treatment for $$L_{s}$$ $$t_{2j}^{\alpha}$$ Rate of treatment for $$I_{j}$$. Note that $$t_{2x}$$ is the rate of successful treatment of $$I_{x}, j\in\{x,m,s\}$$ Parameters Interpretation $$d^{\alpha}$$ Death and birth rate $$\lambda_{j}^{\alpha}$$ Rate of infected individuals move to $$L_{j}$$ with strain $$j \in\{s ,m,x\}\$$ $$1-\lambda_{j}^{\alpha}$$ Rate of newly infected individuals progressing to active TB with strain $$j \in\{s,m,x \}$$ $$\beta_{j}^{\alpha}$$ Transmission coefficient with strain $$j \in\{s,m,x \}$$ $$\varepsilon_{j}^{\alpha}$$ Rate of endogenous reactivation of $$L_{j}$$ $$\gamma_{j}^{\alpha}$$ Rate of natural recovery to the latent stage $$L_{j}$$ $$\delta_{j}^{\alpha}$$ Rate of death due to $$TB$$ of strain $$j$$ $$\alpha_{j1}^{\alpha},\alpha_{j2}^{\alpha}$$ Rate of exogenous reinfection of $$L_{j1}$$ due to contact with $$I_{j2}$$ $$1-\sigma_{j}^{\alpha}$$ Efficiency of treatment in preventing infection with strain $$j$$ $$P_{1}^{\alpha}$$ Probability of treatment success for $$L_{s}$$ $$1-P_{1}^{\alpha}$$ Proportion of treated $$L_{s}$$ moved to $$L_{m}$$ due to incomplete treatment or lack of strict compliance in the use of drugs $$P_{2}^{\alpha}$$ Probability of treatment success for $$I_{s}$$ $$1-P_{2}^{\alpha}$$ Proportion of treated $$I_{s}$$ moved to $$L_{m}$$ due to incomplete treatment or lack of strict compliance in the use of drugs $$P_{3}^{\alpha}$$ Probability of treatment success for $$I_{m}$$ $$1-P_{3}^{\alpha}$$ proportion of treated $$I_{m}$$ moved to $$L_{x}$$due to incomplete treatment or lack of strict compliance in the use of drugs $$t_{1s}^{\alpha}$$ Rate of treatment for $$L_{s}$$ $$t_{2j}^{\alpha}$$ Rate of treatment for $$I_{j}$$. Note that $$t_{2x}$$ is the rate of successful treatment of $$I_{x}, j\in\{x,m,s\}$$ Two control functions $$u_{1}(\cdot),$$$$u_{2}(\cdot),$$ and two real positive constants $$\epsilon_{1},$$$$\epsilon_{2},$$ will be added to the model. The controls $$u_{1},$$$$u_{2}$$ represent the effort in preventing the failure of treatment in active TB infectious individuals $$I_{m},$$$$I_{x},$$ respectively, e.g., supervising the patients, helping them to take the TB medications regularly and to complete the TB treatment. The parameters $$\epsilon_{i}$$$$\in$$$$(0, 1)$$, $$i = 1, 2,$$ measure the effectiveness of the controls $$u_{k},$$$$k = 1, 2,$$ respectively, i.e., these parameters measure the efficacy of treatment interventions for active TB infectious individuals $$I_{m},$$$$I_{x}$$. The resulting model is given by a system of nonlinear fractional order differential equations with time-delay memory as follows: \begin{align} ^{c}_{a}D^{\alpha}_{t}S=&\,d^{\alpha}N-d^{\alpha}S-\beta_{s}^{\alpha}\frac{SI_{s}}{N}-\beta_{m}^{\alpha}\frac{SI_{m}}{N}-\beta_{x}^{\alpha}\frac{SI_{x}}{N}, \label{1} \\ \end{align} (2) \begin{align} ^{c}_{a}D^{\alpha}_{t}L_{s}=&\,\lambda_{s}^{\alpha}\beta_{s}^{\alpha}\frac{SI_{s}}{N}+\sigma_{s}^{\alpha}\lambda_{s}^{\alpha}\beta_{s}^{\alpha}\frac{RI_{s}}{N}-\alpha_{ss}^{\alpha}\beta_{s}^{\alpha}\frac{L_{s}I_{s}}{N}-\alpha_{sm}^{\alpha}\beta_{m}^{\alpha}\frac{L_{s}I_{m}}{N}-\alpha_{sx}^{\alpha}\beta_{x}^{\alpha}\frac{L_{s}I_{x}}{N}+\gamma_{s}^{\alpha}I_{s}\notag\\ &-(d^{\alpha}+\epsilon_{s}^{\alpha}+t_{1s}^{\alpha})L_{s},\\ \end{align} (3) \begin{align} ^{c}_{a}D^{\alpha}_{t}L_{m}=&\,\lambda_{m}^{\alpha}\beta_{m}^{\alpha}\frac{SI_{m}}{N}+\sigma_{m}^{\alpha}\lambda_{m}^{\alpha}\beta_{m}^{\alpha}\frac{RI_{m}}{N}+\alpha_{sm}^{\alpha}\beta_{m}^{\alpha}\lambda_{m}^{\alpha}\frac{L_{s}I_{m}}{N}-\alpha_{mm}^{\alpha}\beta_{m}^{\alpha}\frac{L_{m}I_{m}}{N}-\alpha_{mx}^{\alpha}\beta_{x}^{\alpha}\frac{L_{m}I_{x}}{N}\notag\\ &+\gamma_{m}^{\alpha}I_{m}-(d^{\alpha}+\epsilon_{m}^{\alpha})L_{m}+t_{1s}^{\alpha}L_{s}-P_{1}^{\alpha}t_{1s}^{\alpha}L_{s}+t_{2s}^{\alpha}I_{s}%\notag\\& -P_{2}^{\alpha}t_{2s}^{\alpha}I_{s},\\ \end{align} (4) \begin{align} ^{c}_{a}D^{\alpha}_{t}L_{x}=&\,\lambda_{x}^{\alpha}\beta_{x}^{\alpha}\frac{SI_{x}}{N}+\sigma_{x}^{\alpha}\lambda_{x}^{\alpha}\beta_{x}^{\alpha}\frac{RI_{x}}{N}+\alpha_{sx}^{\alpha}\beta_{x}^{\alpha}\lambda_{x}^{\alpha}\frac{L_{s}I_{x}}{N} +\alpha_{mx}^{\alpha}\beta_{x}^{\alpha}\lambda_{x}^{\alpha}\frac{L_{m}I_{x}}{N}-\alpha_{xx}^{\alpha}\beta_{x}^{\alpha}\frac{L_{x}I_{x}}{N}\notag\\ &-(d^{\alpha}+\epsilon_{x}^{\alpha})L_{x}+\gamma_{x}^{\alpha}I_{x}+t_{2m}^{\alpha}I_{m}(t-\tau_{1})+\epsilon_{1}u_{1}(t)I_{m}-P_{3}^{\alpha}t_{2m}^{\alpha}I_{m},\\ \end{align} (5) \begin{align} ^{c}_{a}D^{\alpha}_{t}I_{s}=&\,\alpha_{ss}^{\alpha}\beta_{s}^{\alpha}\frac{L_{s}I_{s}}{N}+(1-\lambda_{s}^{\alpha})\beta_{s}^{\alpha} \left(\frac{SI_{s}}{N}+\sigma_{s}^{\alpha}\frac{RI_{s}}{N}\right)+\epsilon_{s}^{\alpha}L_{s}-(d^{\alpha}+\delta_{s}^{\alpha}+t_{2s}^{\alpha}+\gamma_{s}^{\alpha}),\\ \end{align} (6) \begin{align} ^{c}_{a}D^{\alpha}_{t}I_{m}=&\,\alpha_{mm}^{\alpha}\beta_{m}^{\alpha}\frac{L_{m}I_{m}}{N}+(1-\lambda_{m}^{\alpha})\beta_{m}^{\alpha} \left(\frac{SI_{m}}{N}+\sigma_{m}^{\alpha}\frac{RI_{m}}{N}+\alpha_{sm}^{\alpha}\frac{L_{s}I_{m}}{N}\right)+\epsilon_{m}^{\alpha}L_{m} \notag\\& -(d^{\alpha}+\delta_{m}^{\alpha}+\epsilon_{1}u_{1}(t)+\gamma_{m}^{\alpha})I_{m}-t_{2m}^{\alpha}I_{m}(t-\tau_{1}),\\ \end{align} (7) \begin{align} ^{c}_{a}D^{\alpha}_{t}I_{x}=&\,\alpha_{xx}^{\alpha}\beta_{x}^{\alpha}\frac{L_{x}I_{x}}{N}+(1-\lambda_{x}^{\alpha})\beta_{x}^{\alpha} \left(\frac{SI_{x}}{N}+\sigma_{x}^{\alpha}\frac{RI_{x}}{N}+\alpha_{sx}^{\alpha}\frac{L_{s}I_{x}}{N}+\alpha_{mx}^{\alpha}\frac{L_{m}I_{x}}{N}\right)+\epsilon_{x}^{\alpha}L_{x} \notag\\&-(d^{\alpha}+\delta_{x}^{\alpha}+\gamma_{x}^{\alpha}+\epsilon_{2}u_{2}(t))I_{x}-t_{2x}^{\alpha}I_{x}(t-\tau_{2}),\\ \end{align} (8) \begin{align} ^{c}_{a}D^{\alpha}_{t}R=&\,P_{1}^{\alpha}t_{1s}^{\alpha}L_{s}+P_{2}^{\alpha}t_{2s}^{\alpha}I_{s}+P_{3}^{\alpha}t_{2m}^{\alpha}I_{m}+t_{2x}^{\alpha}I_{x}(t-\tau_{2})+\epsilon_{2}u_{2}(t)I_{x}-\sigma_{s}^{\alpha}\beta_{s}^{\alpha}\frac{RI_{s}}{N} \notag\\&-\sigma_{m}^{\alpha}\beta_{m}^{\alpha}\frac{RI_{m}}{N}-\sigma_{x}^{\alpha}\beta_{x}^{\alpha}\frac{RI_{x}}{N} -d^{\alpha}R \label{8}. \end{align} (9) Where $$^{c}_{a}D^{\alpha}_{t}$$ is the Caputo fractional derivative (Sweilam & AL-Mekhlafi, 2016b). Because the model (2)–(9) monitors the dynamics of human populations, all the parameters are assumed to be nonnegative. The initial conditions for system (2)–(9) are: $$S(\vartheta)=\theta_{1}(\vartheta),$$$$L_{s}(\vartheta)=\theta_{2}(\vartheta),$$$$L_{m}(\vartheta)=\theta_{3}(\vartheta),$$$$L_{x}(\vartheta)=\theta_{4}(\vartheta),$$$$I_{s}(\vartheta)=\theta_{5}(\vartheta),$$$$I_{m}(\vartheta)=\theta_{6}(\vartheta),$$$$I_{x}(\vartheta)=\theta_{7}(\vartheta),$$$$R(\vartheta)=\theta_{8}(\vartheta),$$$$\vartheta\in[-\tau,0],$$ where, $$\theta=(\theta_{1},\theta_{2}, \cdot\cdot\cdot, \theta_{8},)^{T}$$$$\in$$$$C,$$ where $$C$$ is the Banach space $$C([0,\tau],\mathbb{R}^{8})$$. From biological meaning, we further assume that $$\theta_{i}>0$$ for $$i = 1, \cdot\cdot\cdot, 8.$$ Throughout this article, we focus on the dynamics of the solutions of (2)–(9) in the restricted region: $${\it{\Omega}}=\{(S,L_{s},L_{m},L_{x},I_{s},I_{m},I_{x},R)\in \mathbb{R}^{8}|S+L_{s}+L_{m}+L_{x}+I_{s}+I_{m}+I_{x}+R=N\}.$$ It is worth mentioning that, the usual local existence, uniqueness and continuation results hold in this region as prove in (Hale, 1977; Kuang, 1993). Hence, a unique solution $$(S(t)$$,$$L_{s}(t)$$, $$L_{m}(t)$$, $$L_{x}(t),$$$$I_{s}(t)$$, $$I_{m}(t),$$$$I_{x}(t),$$$$R(t))$$ of (2)–(9) with the given initial conditions exists for all time $$t \geq 0$$. Consider the solutions of (2)–(9) with $$(S,L_{s},L_{m},L_{x},I_{s},I_{m},I_{x},R)$$$$\in$$$${\it{\Omega}}'$$, where $${\it{\Omega}}'$$, is the interior of $${\it{\Omega}},$$ for all $$\vartheta\in[-\tau,0]$$. Then the solutions stay in the interior of the region for all time $$t\geq 0,$$ i.e., the region is positively invariant with respect to system (2)–(9) (see e.g., Kuang, 1993). It is known that, a mathematical model has a disease free equilibrium if it has an equilibrium point at which the population remains in the absence of the disease (Van den Driessche & Watmough, 2002). The model (2)–(9) has a disease free equilibrium given by $$E_{0}=(N, 0,0,0,0,0,0,0)$$. The basic reproduction number ($$R_{0}$$) represents the expected average number of new TB infections produced by a single TB active infected individual when in contact with a completely susceptible population (Van den Driessche & Watmough, 2002). For the model (2)–(9) $$R_{0}$$ is given as follows: \begin{align} R_{0}&=max(R_{0s},R_{0m},R_{0x}),\ \ \ \ \ \ where: \\ R_{0s}&=\frac{\beta_{s}^{\alpha}(\varepsilon_{s}^{\alpha}+(1-\lambda_{s}^{\alpha})(d^{\alpha}+t_{1s}^{\alpha}))}{(\varepsilon_{s}^{\alpha}+d^{\alpha}+t_{1s}^{\alpha})(t_{2s}^{\alpha}+\delta_{s}^{\alpha}+d^{\alpha})+\gamma_{s}^{\alpha}(t_{1s}^{\alpha}+d^{\alpha})},\notag\\ R_{0m}&=\frac{\beta_{m}^{\alpha}(\varepsilon_{m}^{\alpha}+(1-\lambda_{m}^{\alpha})d^{\alpha})}{(\varepsilon_{m}^{\alpha}+d^{\alpha})(t_{2m}^{\alpha}+\delta_{m}^{\alpha}+d^{\alpha})+d^{\alpha}\gamma_{m}^{\alpha}},\notag\\ R_{0x}&=\frac{\beta_{x}^{\alpha}(\varepsilon_{x}^{\alpha}+(1-\lambda_{x}^{\alpha})d^{\alpha})}{(\varepsilon_{x}^{\alpha}+d^{\alpha})(t_{2x}^{\alpha}+\delta_{x}^{\alpha}+d^{\alpha})+d^{\alpha}\gamma_{x}^{\alpha}}.\notag \end{align} (10) 3.1. Equilibrium points and their asymptotic stability To discuss the local asymptotic stability for evaluate the equilibrium points, let us consider the following (El-Saka, 2013): $$D^{\alpha}_{t}S=D^{\alpha}_{t}L_{s}=D^{\alpha}_{t}L_{m}=D^{\alpha}_{t}L_{x}=D^{\alpha}_{t}I_{s}=D^{\alpha}_{t}I_{m}=D^{\alpha}_{t}I_{x}=D^{\alpha}_{t}R=0$$. Then, from 1, $$g_{i}(\bar{S}, \bar{L_{s}},\bar{L_{m}}, \bar{L_{x}},\bar{I_{s}}, \bar{I_{m}},\bar{I_{x}},\bar{R})=0$$, $$i=1,2,3,...,8,$$ where $$(\bar{S}, \bar{L_{s}},\bar{L_{m}}, \bar{L_{x}},\bar{I_{s}}, \bar{I_{m}},\bar{I_{x}},\bar{R})$$ denotes any equilibrium point. 3.2. Stability of the disease free equilibrium point Let us consider $$I_s(t)=I_m(t)=I_x(t)=0,$$ then $$L_s(t)=L_m(t)=L_x(t)=0$$, $$R(t)=0,$$ and $$S(t)=N.$$ Then the disease free equilibrium point is $$E_{0}=\{N,0,0,0,0,0,0,0)\}$$. Let us consider the coordinate transformation: $$s(t)=S(t)-\bar{S},$$$$l_{s}(t)=L_{s}(t)-\bar{L_{s}},$$$$l_{m}(t)=L_{m}(t)-\bar{L_{m}},$$$$l_{x}(t)=L_{x}(t)-\bar{L_{x}},$$$$i_{s}(t)=I_{s}(t)-\bar{I_{s}},$$$$i_{m}(t)=I_{m}(t)-\bar{I_{m}},$$$$i_{x}(t)=I_{x}(t)-\bar{I_{x}},$$$$r(t)=R(t)-\bar{R}$$. The corresponding characteristic equation for the disease free equilibrium point is given as follows: $$\displaylines{ (J({E_0}) - \lambda I) = \cr \left( {\matrix{ {\lambda - a} & 0 & 0 & 0 & b & c & {{d_1}} & 0 \cr 0 & {\lambda - e} & 0 & 0 & f & 0 & 0 & 0 \cr 0 & g & {\lambda - h} & 0 & p & q & 0 & 0 \cr 0 & 0 & 0 & {\lambda - r} & 0 & {t_{2m}^\alpha {e^{( - \lambda {\tau _1})}} - p_3^\alpha t_{2m}^\alpha } & t & 0 \cr 0 & u & 0 & 0 & {\lambda - v} & 0 & 0 & 0 \cr 0 & 0 & w & 0 & 0 & {\lambda - x + t_{2m}^\alpha {e^{( - \lambda {\tau _1})}}} & 0 & 0 \cr 0 & 0 & 0 & y & 0 & 0 & {\lambda - z + t_{2x}^\alpha {e^{( - \lambda {\tau _2})}}} & 0 \cr 0 & m & 0 & 0 & n & j & {t_{2x}^\alpha {e^{( - \lambda {\tau _2})}}} & {\lambda - a} \cr } } \right), \cr}$$ where $$a=-d^{\alpha},$$$$b=-\beta_{s}^{\alpha},$$$$c=-\beta_{m}^{\alpha},$$$$d_{1}=-\beta_{x}^{\alpha},$$$$e=-(d^{\alpha}+\varepsilon_{s}^{\alpha}+t_{1s}^{\alpha}),$$$$f=\gamma_{s}^{\alpha}+\lambda_{s}^{\alpha}\beta_{s}^{\alpha},$$$$g=(1-p_{1}^{\alpha})t_{1s}^{\alpha},$$$$h=-(d^{\alpha}+\varepsilon_{m}^{\alpha}),$$$$p=(1-p_{2}^{\alpha})t_{2s}^{\alpha},$$$$q=\gamma_{m}^{\alpha}+\lambda_{m}^{\alpha}\beta_{m}^{\alpha},$$$$r=-(d^{\alpha}+\varepsilon_{x}^{\alpha}),$$$$t=\gamma_{x}^{\alpha}+\lambda_{x}^{\alpha}\beta_{x}^{\alpha},$$$$u=\varepsilon_{s}^{\alpha},$$$$v=-(d^{\alpha}+\delta_{s}^{\alpha}+t_{2s}^{\alpha}+\gamma_{s}^{\alpha}),$$$$w=\varepsilon_{m}^{\alpha},$$$$x=-(d^{\alpha}+\delta_{m}^{\alpha)}+\gamma_{m}^{\alpha}),$$$$y=\varepsilon_{x}^{\alpha},$$$$z=-(d^{\alpha}+\delta_{x}^{\alpha}+\gamma_{x}^{\alpha}),$$$$m=p_{1}^{\alpha}t_{1s}^{\alpha},$$$$n=p_{2}^{\alpha}t_{2s}^{\alpha},$$$$j=p_{3}^{\alpha}t_{2m}^{\alpha}$$. The characteristic equation associated with above matrix is given by El-Saka (2013) as follows: $${\it{\Delta}}(\lambda)=|J(E_{0})-\lambda I|=0,$$ then we have: \begin{align} &(a-\lambda)^{2}(\lambda^{2}-(r+z+t_{2x}^{\alpha}e^{(-\lambda\tau_{2})})\lambda-yt+(z-t_{2x}^{\alpha}e^{(-\lambda\tau_{2})})r)(\lambda^2-(h+x-t_{2m}^{\alpha}e^{(-\lambda\tau_{1})})\lambda \notag \\ &+(x+t_{2m}^{\alpha}e^{(-\lambda\tau_{2})})h-wq)(\lambda^{2}-(e+v)\lambda-uf+ve)=0 \label{RR}. \end{align} (11) Lemma 1 If $$R_{0}<1$$, then the disease free equilibrium point $$E_{0}$$ is locally asymptotically stable for $$\tau_{1}=\tau_{2}=0$$. – Proof. When $$\tau_{1}=\tau_{2}=0,$$ the associated transcendental characteristic equation of system (2)–(9) at $$E_{0}$$ become $${\it{\Delta}}(\lambda)=P(\lambda)=0.$$ Then the eigenvalues of the Jacobian matrix are: \begin{eqnarray*} \lambda_{1,2} &=& -d,\\ \lambda_{3,4}&=&\frac{r+(z-t_{2x}^{\alpha})\pm \sqrt{(r^2-2(z-t_{2x}^{\alpha})r+(z-t_{2x}^{\alpha})^2+4yt)}}{2},\\ \lambda_{5,6}&=&\frac{x-t_{2m}^{\alpha}+h\pm \sqrt{(x-t_{2m}^{\alpha})^2-2(x-t_{2m}^{\alpha})h+h^2+4wq)}}{2},\\ \lambda_{7,8}&=&\frac{v+e\pm \sqrt{(v^2-2ve+e^2+4uf)}}{2}, \end{eqnarray*} these roots are negative or have negative real parts and all eigenvalues satisfies Matignon’s conditions (Matignon, 1996), given by $$\left(|arg\lambda_{i}|>\frac{\alpha\pi}{2}\right)\!,$$ then the disease free equilibrium $$E_{0}$$ is locally asymptotically stable. □ Lemma 2 Let $$R_{0} < 1,$$ then the disease free equilibrium $$E_{0}$$ is locally asymptotically stable for $$\tau_{1}>0$$ and $$\tau_{2}>0$$. Proof. Let us consider $$\tau_{1}>0$$ and $$\tau_{2}>0$$, we noted that second and third factor of the characteristic equation (11), which are $$(\lambda^{2}-(r+z+t_{2x}^{\alpha}e^{(-\lambda\tau_{2})})\lambda-yt+(z-t_{2x}^{\alpha}e^{(-\lambda\tau_{2})})r)$$ and $$(\lambda^2-(h+x-t_{2m}^{\alpha}e^{(-\lambda\tau_{1})})\lambda \&+(x+t_{2m}^{\alpha}e^{(-\lambda\tau_{1})})h-wq),$$ have no pure imaginary roots for any value of the delays $$\tau_{1}, \tau_{2}$$ if $$R_{0} < 1$$. Hence all the roots of the characteristic equation have negative real parts and we get that the disease free equilibrium is locally asymptotically stable regardless of the value of the delay and all eigenvalues satisfies Matignon’s conditions (Matignon, 1996), given by $$\left(|arg\lambda_{i}|>\frac{\alpha\pi}{2}\right)$$ then the disease free equilibrium $$E_{0}$$ is locally asymptotically stable. □ 3.3. Stability of the endemic equilibrium In the absence of controls, i.e., $$u_{1}=u_{2}=0,$$ the system (2)–(9) has endemic equilibrium if at least one of the infected variables is not zero. The expression of analytic for this point is cumbersome and not useful for our purposes (Silva et al., 2017). However, we will determine the stability of the endemic equilibrium of the proposed system numerically. Consider the values of parameters from Table 2. Then the basic reproduction number is $$R_{0}>1$$. The endemic equilibrium $$S=1.0781\times10^{3},$$$$L_{s}=0.0134,$$$$L_{m}=2.1814\times10^{4},$$$$L_{x}=1.8583\times10^{4},$$$$I_{s}=0.0143,$$$$I_{m}=1.5252\times10^{3},$$$$I_{x}=1.1477\times10^{3},$$$$R=1.5852\times10^{4}.$$ The matrices $$A_{1}$$ and $$A_{2}$$ associated to the linearized system at the endemic equilibrium are computed as follows: $${\small A_{1}=\left(\!\!\!\!\begin{array}{cccccccc} -0.6373 & 0 & 0 & 0 & -0.2516 & -0.2516 & -0.2516& 0\\ 0 & -2.0494 & 0 & 0 & 0.5881 & 0 & 0 & 0 \\ 0.1779 & 0.2489& -0.0450 & 0 & 0.2400 & 0.3337 & -0.2545 & 0.0445 \\ 0.1339 & 0.0067 & 0.0067 & -0.0272 & 0 & 0.1200 & 0.4986 & 0.0335\\ 0 & 0.0002 & 0 & 0 &-1.4255& 0 & 0 & 0 \\ 0.0002\times10^{3} & 0 & 0 & 0 & 0.0003\times10^{3} & 1.4095\times10^{3} & 0 & 0.0678 \\ 0.1339 & 0.0067 & 0.0067 & 0.0136& 0 & 0 & 0.0814 & 0.0335 \\ 0 & 1.7600& 0 & 0 &0.8353& -0.0447 & 0.0753 & -0.1695\\ \end{array}\!\!\!\right),}$$ $$A_{2}=\left( \begin{array}{cccccccc} 0& 0 & 0 & 0 & 0& 0 & 0& 0 \\ 0 & 0 & 0 & 0 & 0& 0 & 0 & 0 \\ 0 & 0&0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0& 0 & t_{2m}^{\alpha}e^{-\tau_{1}\lambda} & 0 & 0\\ 0 & 0 & 0 & 0 &0 & 0 & 0 & 0 \\ 0 & 0 & 0& 0 & 0 & -t_{2m}^{\alpha}e^{-\tau_{1}\lambda}& 0 & 0 \\ 0 & 0 & 0 & 0& 0 & 0 & -t_{2x}^{\alpha}e^{-\tau_{2}\lambda}& 0 \\ 0 & 0& 0 & 0 & 0 &0 & t_{2x}^{\alpha}e^{-\tau_{2}\lambda} & 0 \\ \end{array} \right).$$ Table 2 All parameters and their references of the system (2)–(9) Parameters Value Reference $$N$$ 60000 Assumed $$d^{\alpha}$$ $$(1/73.45)^{\alpha}$$$$\left(\frac{1}{year}\right)^{\alpha}$$ http://www.indexmundi.com/EGYPT/life_expectancy_at_birth.html $$\beta_{s}^{\alpha}=\beta_{m}^{\alpha}=\beta_{x}^{\alpha}$$ $$14^{\alpha} \left(\frac{1}{year}\right)^{\alpha}$$ World Health Organization (2012) $$\lambda_{s}^{\alpha}=\lambda_{m}^{\alpha}=\lambda_{x}^{\alpha}$$ $0.5^{\alpha}$$\left(\frac{1}{year}\right)^{\alpha} Sweilam & AL-Mekhlafi (2016b)$$\epsilon_{s}^{\alpha}=\epsilon_{m}^{\alpha}=\epsilon_{x}^{\alpha}0.0002^{\alpha}\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$\alpha_{r1,r2}^{\alpha}0.05^{\alpha}\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$\gamma_{s}^{\alpha}$$=$$\gamma_{m}^{\alpha}$$=$$\gamma_{x}^{\alpha}0.00002^{\alpha}\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$t_{1s}^{\alpha}2^{\alpha}\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$t_{2r}^{\alpha}$$:$$r\in(s,m,x)t_{2s}^{\alpha}=2$$;$$t_{2m}^{\alpha}=t_{2x}^{\alpha}=1\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$\sigma_{r}^{\alpha}$$0.25^{\alpha}$$\left(\frac{1}{year}\right)^{\alpha}$ Sweilam & AL-Mekhlafi (2016b) $$P_{r}^{\alpha}$$ $$0.88^{\alpha}$$$$\left(\frac{1}{year}\right)^{\alpha}$$ Sweilam & AL-Mekhlafi (2016b) $$\delta_{r}^{\alpha}$$ $$0$$ Sweilam & AL-Mekhlafi (2016b) $$T$$ 10 years Assumed Parameters Value Reference $$N$$ 60000 Assumed $$d^{\alpha}$$ $$(1/73.45)^{\alpha}$$$$\left(\frac{1}{year}\right)^{\alpha}$$ http://www.indexmundi.com/EGYPT/life_expectancy_at_birth.html $$\beta_{s}^{\alpha}=\beta_{m}^{\alpha}=\beta_{x}^{\alpha}$$ $$14^{\alpha} \left(\frac{1}{year}\right)^{\alpha}$$ World Health Organization (2012) $$\lambda_{s}^{\alpha}=\lambda_{m}^{\alpha}=\lambda_{x}^{\alpha}$$ $0.5^{\alpha}$$\left(\frac{1}{year}\right)^{\alpha} Sweilam & AL-Mekhlafi (2016b)$$\epsilon_{s}^{\alpha}=\epsilon_{m}^{\alpha}=\epsilon_{x}^{\alpha}0.0002^{\alpha}\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$\alpha_{r1,r2}^{\alpha}0.05^{\alpha}\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$\gamma_{s}^{\alpha}$$=$$\gamma_{m}^{\alpha}$$=$$\gamma_{x}^{\alpha}0.00002^{\alpha}\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$t_{1s}^{\alpha}2^{\alpha}\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$t_{2r}^{\alpha}$$:$$r\in(s,m,x)t_{2s}^{\alpha}=2$$;$$t_{2m}^{\alpha}=t_{2x}^{\alpha}=1\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$\sigma_{r}^{\alpha}$$0.25^{\alpha}$$\left(\frac{1}{year}\right)^{\alpha}$ Sweilam & AL-Mekhlafi (2016b) $$P_{r}^{\alpha}$$ $$0.88^{\alpha}$$$$\left(\frac{1}{year}\right)^{\alpha}$$ Sweilam & AL-Mekhlafi (2016b) $$\delta_{r}^{\alpha}$$ $$0$$ Sweilam & AL-Mekhlafi (2016b) $$T$$ 10 years Assumed Table 2 All parameters and their references of the system (2)–(9) Parameters Value Reference $$N$$ 60000 Assumed $$d^{\alpha}$$ $$(1/73.45)^{\alpha}$$$$\left(\frac{1}{year}\right)^{\alpha}$$ http://www.indexmundi.com/EGYPT/life_expectancy_at_birth.html $$\beta_{s}^{\alpha}=\beta_{m}^{\alpha}=\beta_{x}^{\alpha}$$ $$14^{\alpha} \left(\frac{1}{year}\right)^{\alpha}$$ World Health Organization (2012) $$\lambda_{s}^{\alpha}=\lambda_{m}^{\alpha}=\lambda_{x}^{\alpha}$$ $0.5^{\alpha}$$\left(\frac{1}{year}\right)^{\alpha} Sweilam & AL-Mekhlafi (2016b)$$\epsilon_{s}^{\alpha}=\epsilon_{m}^{\alpha}=\epsilon_{x}^{\alpha}0.0002^{\alpha}\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$\alpha_{r1,r2}^{\alpha}0.05^{\alpha}\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$\gamma_{s}^{\alpha}$$=$$\gamma_{m}^{\alpha}$$=$$\gamma_{x}^{\alpha}0.00002^{\alpha}\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$t_{1s}^{\alpha}2^{\alpha}\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$t_{2r}^{\alpha}$$:$$r\in(s,m,x)t_{2s}^{\alpha}=2$$;$$t_{2m}^{\alpha}=t_{2x}^{\alpha}=1\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$\sigma_{r}^{\alpha}$$0.25^{\alpha}$$\left(\frac{1}{year}\right)^{\alpha}$ Sweilam & AL-Mekhlafi (2016b) $$P_{r}^{\alpha}$$ $$0.88^{\alpha}$$$$\left(\frac{1}{year}\right)^{\alpha}$$ Sweilam & AL-Mekhlafi (2016b) $$\delta_{r}^{\alpha}$$ $$0$$ Sweilam & AL-Mekhlafi (2016b) $$T$$ 10 years Assumed Parameters Value Reference $$N$$ 60000 Assumed $$d^{\alpha}$$ $$(1/73.45)^{\alpha}$$$$\left(\frac{1}{year}\right)^{\alpha}$$ http://www.indexmundi.com/EGYPT/life_expectancy_at_birth.html $$\beta_{s}^{\alpha}=\beta_{m}^{\alpha}=\beta_{x}^{\alpha}$$ $$14^{\alpha} \left(\frac{1}{year}\right)^{\alpha}$$ World Health Organization (2012) $$\lambda_{s}^{\alpha}=\lambda_{m}^{\alpha}=\lambda_{x}^{\alpha}$$ $0.5^{\alpha}$$\left(\frac{1}{year}\right)^{\alpha} Sweilam & AL-Mekhlafi (2016b)$$\epsilon_{s}^{\alpha}=\epsilon_{m}^{\alpha}=\epsilon_{x}^{\alpha}0.0002^{\alpha}\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$\alpha_{r1,r2}^{\alpha}0.05^{\alpha}\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$\gamma_{s}^{\alpha}$$=$$\gamma_{m}^{\alpha}$$=$$\gamma_{x}^{\alpha}0.00002^{\alpha}\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$t_{1s}^{\alpha}2^{\alpha}\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$t_{2r}^{\alpha}$$:$$r\in(s,m,x)t_{2s}^{\alpha}=2$$;$$t_{2m}^{\alpha}=t_{2x}^{\alpha}=1\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$\sigma_{r}^{\alpha}$$0.25^{\alpha}$$\left(\frac{1}{year}\right)^{\alpha}$ Sweilam & AL-Mekhlafi (2016b) $$P_{r}^{\alpha}$$ $$0.88^{\alpha}$$$$\left(\frac{1}{year}\right)^{\alpha}$$ Sweilam & AL-Mekhlafi (2016b) $$\delta_{r}^{\alpha}$$ $$0$$ Sweilam & AL-Mekhlafi (2016b) $$T$$ 10 years Assumed Then the transcendental characteristic equation $${\it{\Delta}}\lambda=(\lambda I-A_{1}-A_{2})$$ is given by: \begin{align} &\lambda^{8}-(e^{\tau_{1}\lambda}+e^{\tau_{2}\lambda}+1.4051\times10^3)e^{\tau_{1}\lambda}e^{-\tau_{2}\lambda}\lambda^{7}-(4.4353e^{\tau_{1}\lambda}-1.4052\times10^3e^{\tau_{2}\lambda}-e^{\tau_{1}\lambda}e^{\tau_{2}\lambda}\notag\\&+6.2450\times10^3) e^{-\tau_{1}\lambda}e^{-\tau_{2}\lambda}\lambda^{6}-(6.5237e^{\tau_{1}\lambda}-6.0835\times10^{3}e^{\tau_{2}\lambda}-4.3204e^{\tau_{1}\lambda}e^{\tau_{2}\lambda}+9.1912\times10^{3})\notag\\&e^{-\tau_{1}\lambda}e^{-\tau_{2}\lambda}\lambda^{5}-(3.7527e^{\tau_{1}\lambda}-8.4571\times10^{3}e^{\tau_{2}\lambda}-6.0014e^{\tau_{1}\lambda}e^{\tau_{2}\lambda}+5.2875\times10^{3})e^{-\tau_{1}\lambda}e^{-\tau_{2}\lambda}\lambda^{4}\notag\\ &-(0.8647e^{\tau_{1}\lambda}-4.1905\times10^{3}e^{\tau_{2}\lambda}-2.9696e^{\tau_{1}\lambda}e^{\tau_{2}\lambda}+1.2155\times10^{3})e^{-\tau_{1}\lambda}e^{-\tau_{2}\lambda}\lambda^{3}\notag\\ &+(0.0776e^{\tau_{1}\lambda}-621.5397e^{\tau_{2}\lambda}-0.4361e^{\tau_{1}\lambda}e^{\tau_{2}\lambda}+106.2920)e^{-\tau_{1}\lambda}e^{-\tau_{2}\lambda}\lambda^{2}\notag\\&-(0.0019e^{\tau_{1}\lambda}-27.9671e^{\tau_{2}\lambda}-0.0179e^{\tau_{1}\lambda}e^{\tau_{2}\lambda}+2.3034)e^{-\tau_{1}\lambda}e^{-\tau_{2}\lambda}\lambda\notag\\&+(0.3511e^{\tau_{2}\lambda}+7.1078\times10^{-6}e^{\tau_{1}\lambda}+1.7185\times10^{-4}e^{\tau_{1}\lambda}e^{\tau_{1}\lambda}+0.0233)=0.\label{2a2} \end{align} (12) Consider now the case $$\tau_{1} > 0$$ and $$\tau_{2} > 0,$$ we noted that, the roots of the characteristic equation (12) have no pure imaginary roots for any value of the delay $$\tau_{1}, \tau_{2}$$ if $$R_{0} > 1$$. Hence all the roots of the characteristic equation have negative real parts and all eigenvalues satisfies Matignon’s conditions (Matignon, 1996). Therefore, the endemic equilibrium is locally asymptotically stable. 4. Formulation of the optimal control problem Let us consider the state system (2)–(9), in $$R^{8}$$ with the set of admissible control functions: $$U=\{(u_{1}(\cdot),u_{2}(\cdot))\in (L^{\infty}(0,T)^{2}) \mid 0\leq u_{1}(\cdot),u_{2}(\cdot)\leq 1,\forall t\in[0,T]\}.$$ The objective functional is given by : \begin{align} J(u_{k})=\int_{0}^{T}\eta(S,L_{s},L_{m},L_{x},I_{s},I_{m},I_{x}, R,u_{k},t)dt,\notag \end{align} where $$k=1,2,$$ and \begin{align} J(u_{1}(\cdot),u_{2}(\cdot))&=\int_{0}^{T}[I_{m}(t)+I_{x}(t)+\frac{1}{2}B_{1}u^{2}_{1}(t)+\frac{1}{2}B_{2}u^{2}_{2}(t)] dt\label{22}, \end{align} (13) \begin{align} J(u^{\ast}_{1}(\cdot),u^{\ast}_{2}(\cdot))&= \min_{{\it{\Omega}}} J(u_{1}(\cdot),u_{2}(\cdot)), \end{align} (14) where the constants $$B_{1}$$, $$B_{2}$$, are the measure of the relative cost of the interventions associated with the controls $$u_{1}$$, $$u_{2}$$, respectively. The main point in FDOCP is to find the optimal controls $$u_{k}(t),$$ where $$k=1,2,$$ which minimizes the objective function (13), subject to the following dynamic system: \begin{align} ^{c}_{a}D^{\alpha}_{t}S&=\xi_{1}(S,L_{s},L_{m},L_{x},I_{s},I_{m\tau},I_{m},I_{x},I_{x\tau}, R,u_{k}, t),\notag\\ ^{c}_{a}D^{\alpha}_{t}L_{s}&=\xi_{2}(S,L_{s},L_{m},L_{x},I_{s},I_{m\tau},I_{m},I_{x},I_{x\tau}, R,u_{k},t),\notag\\ ^{c}_{a}D^{\alpha}_{t}L_{m}&=\xi_{3}(S,L_{s},L_{m},L_{x},I_{s},I_{m\tau},I_{m},I_{x},I_{x\tau}, R,u_{k},t),\notag\\ ^{c}_{a}D^{\alpha}_{t}L_{x}&=\xi_{4}(S,L_{s},L_{m},L_{x},I_{s},I_{m\tau},I_{m},I_{x},I_{x\tau}, R,u_{k},t),\notag\\ ^{c}_{a}D^{\alpha}_{t}I_{s}&=\xi_{5}(S,L_{s},L_{m},L_{x},I_{s},I_{m\tau},I_{m},I_{x},I_{x\tau}, R,u_{k},t),\notag\\ ^{c}_{a}D^{\alpha}_{t}I_{m}&=\xi_{6}(S,L_{s},L_{m},L_{x},I_{s},I_{m\tau},I_{m},I_{x},I_{x\tau}, R,u_{k},t),\notag\\ ^{c}_{a}D^{\alpha}_{t}I_{x}&=\xi_{7}(S,L_{s},L_{m},L_{x},I_{s},I_{m\tau},I_{m},I_{x},I_{x\tau}, R,u_{k},t),\notag\\ ^{c}_{a}D^{\alpha}_{t}R&=\xi_{8}(S,L_{s},L_{m},L_{x},I_{s},I_{m\tau},I_{m},I_{x},I_{x\tau}, R,u_{k},t),\notag \end{align} and satisfying the inital conditions: $$\displaylines{ S(\vartheta ) = {\theta _1}(\vartheta ),\,\,{L_s}(\vartheta ) = {\theta _2}(\vartheta ),\,\,{L_m}(\vartheta ) = {\theta _3}(\vartheta ),\,\,{L_x}(\vartheta ) = {\theta _4}(\vartheta ),\,\,{I_s}(\vartheta ) = {\theta _5}(\vartheta ),\,\,{I_m}(\vartheta ) = {\theta _6}(\vartheta ), \cr {I_x}(\vartheta ) = = {\theta _7}(\vartheta ),\quad R(\vartheta ) = {\theta _8}(\vartheta ), \cr}$$ (15)$$\vartheta\in[-\tau,0],$$ where, $$\theta=(\theta_{1},\theta_{2}, \cdot\cdot\cdot \theta_{8})^{'}.$$ The following expression defines a modified objective function: \begin{align} \tilde{J}&=\int_{0}^{T}[H(S(t),L_{s}(t),L_{m}(t),L_{x}(t),I_{s}(t),I_{m}(t), I_{m}(t-\tau),I_{x}(t),I_{x}(t-\tau), R(t),u_{k}(t))\notag\\ &\quad -\sum^{8}_{i=1}\lambda_{i} \xi_{i}(S(t),L_{s}(t),L_{m}(t),L_{x}(t),I_{s}(t),I_{m}(t),I_{m}(t-\tau),I_{x}(t),I_{x}(t-\tau), R(t),u_{k}(t))]dt\label{m2}, \end{align} (16) where $$H(S,L_{s},L_{m},L_{x},I_{s},I_{m},I_{x}, R,u_{k},t)$$ is the following Hamiltonian: \begin{eqnarray} &&H(S,L_{s}(t),L_{m}(t),L_{x}(t),I_{s}(t),I_{m}(t), I_{m}(t-\tau),I_{x}(t), I_{x}(t-\tau), R(t),u_{k}(t),\lambda_{i})\nonumber\\ &&\quad =\eta(S(t),L_{s}(t),L_{m}(t),L_{x}(t),I_{s}(t),I_{m}(t),I_{m}(t-\tau),I_{x}(t), I_{x}(t-\tau) R(t),u_{k}(t))\nonumber\\ &&\qquad +\sum^{8}_{i=1}\lambda_{i}(t) \xi_{i}(S(t),L_{s}(t),L_{m}(t),L_{x}(t),I_{s}(t),I_{m}(t),I_{m}(t-\tau),I_{x}(t),I_{x}(t-\tau), R(t),u_{k}(t)). \label{m3} \end{eqnarray} (17) From (16) and (17), we can derive the necessary and sufficient conditions for FDOPC (Agarwal, 2008) as follows: \begin{eqnarray} &&^{c}_{t}D^{\alpha}_{t_{f}}\lambda_{1}=\frac{\partial H}{\partial S}, \ \ ^{c}_{t}D^{\alpha}_{t_{f}}\lambda_{2}=\frac{\partial H}{\partial L_{s}},\nonumber\\ &&^{c}_{t}D^{\alpha}_{t_{f}}\lambda_{3}=\frac{\partial H}{\partial L_{m}},\ \ ^{c}_{t}D^{\alpha}_{t_{f}}\lambda_{4} =\frac{\partial H}{\partial L_{x}},\label{mh}\\ &&^{c}_{t}D^{\alpha}_{t_{f}}\lambda_{5}= \frac{\partial H}{\partial I_{s}},\ \ ^{c}_{t}D^{\alpha}_{t_{f}}\lambda_{6}=\frac{\partial H}{\partial I_{m}} +X_{[0,t_{f}-\tau]}\frac{\partial H}{\partial I_{m}(t-\tau)}(t+\tau),\nonumber\\ &&^{c}_{t}D^{\alpha}_{t_{f}}\lambda_{7}=\frac{\partial H}{\partial I_{x}(t)}+X_{[0,t_{f}-\tau]}\frac{\partial H}{\partial I_{x}(t-\tau)} (t+\tau),\ \ ^{c}_{t}D^{\alpha}_{t_{f}}\lambda_{8}=\frac{\partial H}{\partial R},\nonumber \end{eqnarray} (18) ${{X}_{[0,{{t}_{f}}-\tau ]}}=\left\{ \begin{array}{*{35}{l}} 1, & t\in \left[ 0,{{t}_{f}}-\tau \right], \\ 0, & otherwise. \\ \end{array} \right.$ (19) \begin{gather} 0=\frac{\partial H}{\partial u_{k}} \label{m4},\\ ^{c}_{0}D^{\alpha}_{t}S=\frac{\partial H}{\partial \lambda_{1}},\ \ ^{c}_{0}D^{\alpha}_{t}L_{s}=\frac{\partial H}{\partial \lambda_{2}},\nonumber\\ ^{c}_{0}D^{\alpha}_{t}L_{m}=\frac{\partial H}{\partial \lambda_{3}},\ \ ^{c}_{0}D^{\alpha}_{t}L_{x}=\frac{\partial H}{\partial \lambda_{4}},\\ ^{c}_{0}D^{\alpha}_{t}I_{s}=\frac{\partial H}{\partial \lambda_{5}}, \ \ ^{c}_{0}D^{\alpha}_{t}I_{m}=\frac{\partial H}{\partial \lambda_{6}},\nonumber\\ ^{c}_{0}D^{\alpha}_{t}I_{x}=\frac{\partial H}{\partial \lambda_{7}}, \ \ ^{c}_{0}D^{\alpha}_{t}R=\frac{\partial H}{\partial \lambda_{8}},\nonumber \end{gather} (20) and it is also required that: $$\lambda_{i}(T)=0, where \ \ i=1,2,3,...,8 \label{m5}.$$ (21) Equations (19) and (21) describe the necessary conditions in terms of a Hamiltonian for the optimal control problem defined above. Theorem 4.1 There exists an optimal pair $$u^{\ast}_{1}(\cdot),$$$$u^{\ast}_{2}(\cdot),$$ corresponding solutions $$S^{\ast}(\cdot),$$$$L^{\ast}_{s}(\cdot),$$$$L^{\ast}_{m}(\cdot),$$$$L^{\ast}_{x}(\cdot),$$$$I^{\ast}_{s}(\cdot),$$$$I^{\ast}_{m}(\cdot),$$$$I^{\ast}_{x}(\cdot),$$ and $$R^{\ast}(\cdot)$$ and that minimizes $$J(u_{1}(\cdot),u_{2}(\cdot))$$ over $${\it{\Omega}}$$. The explicit optimal controls are connected to the existence of continuous specific functions $$\lambda_{i}$$ for $$i=1,2,3,...,8$$ satisfying the follows: (i) adjoint equations: \begin{align} ^{c}_{t}D^{\alpha}_{t_{f}}\lambda_{1}^{\ast}=&-(\lambda_{1}^{\ast}(t)(\frac{\beta_{s}^{\alpha}}{N}I_{s}^{\ast}(t) +\frac{\beta_{m}^{\alpha}}{N}I_{m}^{\ast}(t)+\frac{\beta_{x}^{\alpha}}{N}I_{x}^{\ast}(t) +d^{\alpha})-\lambda_{2}^{\ast}(t)\frac{\lambda_{s}^{\alpha}\beta_{s}^{\alpha}}{N}I_{s}^{\ast}(t) -\lambda_{3}^{\ast}(t)\frac{\lambda_{m}^{\alpha}\beta_{m}^{\alpha}}{N}I_{m}^{\ast}(t) \notag\\&-\lambda_{4}^{\ast}(t)\frac{\lambda_{x}^{\alpha}\beta_{x}^{\alpha}}{N}I_{s}^{\ast}(t) -\lambda_{5}^{\ast}(t)\frac{(1-\lambda_{s}^{\alpha})\beta_{s}^{\alpha}}{N}I_{s}^{\ast}(t) -\lambda_{6}^{\ast}(t)\frac{(1-\lambda_{m}^{\alpha})\beta_{m}^{\alpha}}{N}I_{m}^{\ast}(t)\notag\\ &-\lambda_{7}^{\ast}(t)\frac{(1-\lambda_{x}^{\alpha})\beta_{x}}{N}I_{x}^{\ast}(t))\label{39},\\ \end{align} (22) \begin{align} ^{c}_{t}D^{\alpha}_{t_{f}}\lambda_{2}^{\ast}=&-(-1+\lambda_{2}^{\ast}(t)(\frac{\beta_{s}^{\alpha}\alpha_{ss}^{\alpha}I_{s}^{\ast}(t)}{N} +\frac{\beta_{m}^{\alpha}\alpha_{sm}^{\alpha}I_{m}^{\ast}(t)}{N} +\frac{\beta_{x}^{\alpha}\alpha_{sx}^{\alpha}I_{x}^{\ast}(t)}{N}) +\lambda_{2}^{\ast}(t)(d^{\alpha}+\epsilon_{s}^{\alpha}+t_{1s}^{\alpha}\notag\\ &+\epsilon_{1}u_{1}^{\ast}(t))-\lambda_{3}^{\ast}(t) \frac{\lambda_{m}^{\alpha}\beta_{m}^{\alpha}\alpha_{sm}^{\alpha}}{N}I_{m}^{\ast}(t) -(t_{1s}^{\alpha}+\epsilon_{1}u_{1}^{\ast}(t)-P_{1}^{\alpha}t_{1s}^{\alpha})\lambda_{3}^{\ast}(t)\notag\\ & -\lambda_{4}^{\ast}(t)\frac{\lambda_{x}^{\alpha}\beta_{x}^{\alpha}\alpha_{sx}^{\alpha}}{N}I_{s}^{\ast}(t) -\lambda_{5}^{\ast}(t)\frac{\beta_{s}^{\alpha}\alpha_{ss}^{\alpha}}{N}I_{s}^{\ast}(t)-\lambda_{5}^{\ast}(t)\epsilon_{s}^{\alpha}\notag\\ &-\lambda_{6}^{\ast}(1-\lambda_{m}^{\alpha})\frac{\beta_{m}^{\alpha}\alpha_{sm}^{\alpha}}{N}I_{m}^{\ast}(t)-\lambda_{7}^{\ast}(1-\lambda_{x}^{\alpha})\frac{\beta_{x}^{\alpha}\alpha_{sx}^{\alpha}}{N}I_{x}^{\ast}(t)-P_{1}^{\alpha}t_{1s}^{\alpha}\lambda_{8}^{\ast}(t)),\\ \end{align} (23) \begin{align} ^{c}_{t}D^{\alpha}_{t_{f}}\lambda_{3}^{\ast}=&-\left(\lambda_{3}^{\ast}(t) (\alpha_{mm}^{\alpha}\beta_{m}^{\alpha}\frac{I_{m}^{\ast}(t)}{N}+\alpha_{mx}^{\alpha}\beta_{x}^{\alpha} \frac{I_{x}^{\ast}(t)}{N}+d^{\alpha}+\epsilon_{m}^{\alpha})-\lambda_{4}^{\ast}(t)\alpha_{mx}^{\alpha} \beta_{x}^{\alpha}\lambda_{x}^{\alpha}\frac{I_{x}^{\ast}(t)}{N}\right.\nonumber\\ &\left.-\lambda_{6}^{\ast}(t)\left(\alpha_{mm}^{\alpha}\beta_{m}^{\alpha}\frac{I_{m}^{\ast}(t)}{N}+\epsilon_{m}^{\alpha}\right)-\lambda_{7}^{\ast}(t) (1-\lambda_{x}^{\alpha})\alpha_{mx}^{\alpha}\beta_{^{\alpha}}\frac{I_{x}^{\ast}(t)}{N}\right),\\ \end{align} (24) \begin{align} ^{c}_{t}D^{\alpha}_{t_{f}}\lambda_{4}^{\ast}=&-(\lambda_{4}^{\ast}(t)(\alpha_{xx}^{\alpha}\beta_{x}^{\alpha}\frac{I_{x}^{\ast}(t)}{N} +d^{\alpha}+\epsilon_{x}^{\alpha})-\lambda_{7}^{\ast}(t)(\alpha_{xx}^{\alpha}\beta_{x}^{\alpha}\frac{I_{x}^{\ast}(t)}{N} -\epsilon_{x}^{\alpha}),\\ \end{align} (25) \begin{align} ^{c}_{t}D^{\alpha}_{t_{f}}\lambda_{5}^{\ast}=&-(-1+\lambda_{1}^{\ast}(t)\beta_{s}^{\alpha}\frac{S^{\ast}(t)}{N}+\lambda_{2}^{\ast}(t)(\alpha_{ss}^{\alpha}\beta_{s}^{\alpha}\frac{L_{s}^{\ast}(t)}{N}-\lambda_{s}^{\alpha}\beta_{s}^{\alpha}\frac{S^{\ast}(t)}{N}-\lambda_{s}^{\alpha}\sigma_{s}^{\alpha}\beta_{s}^{\alpha}\frac{R^{\ast}(t)}{N}-\gamma_{s}^{\alpha})\notag\\&+\lambda_{5}^{\ast}(t)((d^{\alpha}+\delta_{s}^{\alpha}+t_{2s}^{\alpha}+\gamma_{s}^{\alpha})-(1-\lambda_{s}^{\alpha})\beta_{s}^{\alpha}\frac{S^{\ast}(t)}{N}-\sigma_{s}^{\alpha}(1-\lambda_{s}^{\alpha})\beta_{s}^{\alpha}\frac{R^{\ast}(t)}{N}\notag\\&+\frac{\beta_{s}^{\alpha}\alpha_{ss}^{\alpha}}{N}L_{s}^{\ast}(t))-\lambda_{3}^{\ast}(t)(t_{2s}^{\alpha}-P_{2}^{\alpha}t_{2s}^{\alpha})+\lambda_{8}^{\ast}(t)(\sigma_{s}^{\alpha}\beta_{s}^{\alpha}\frac{R^{\ast}(t)}{N}-P_{2}^{\alpha}t_{2s}^{\alpha})),\\ \end{align} (26) \begin{align} ^{c}_{t}D^{\alpha}_{t_{f}}\lambda_{6}=&- \left(-1+\lambda_{1}^{\ast}(t)\beta_{m}^{\alpha}\frac{S^{\ast}(t)}{N}+\lambda_{2}^{\ast}(t)\left(\alpha_{sm}^{\alpha}\beta_{m}^{\alpha}\frac{L_{s}^{\ast}(t)}{N}-\lambda_{3}^{\ast}(t)(\lambda_{m}^{\alpha}\beta_{m}^{\alpha}\frac{S^{\ast}(t)}{N}+\lambda_{m}^{\alpha}\sigma_{m}^{\alpha}\beta_{m}^{\alpha}\frac{R^{\ast}(t)}{N}\right.\right.\notag\\ &+\frac{\lambda_{m}^{\alpha}\beta_{m}^{\alpha}\alpha_{sm}^{\alpha}}{N}L_{s}^{\ast}(t)-\frac{\beta_{m}^{\alpha}\alpha_{mm}^{\alpha}}{N}L_{m}^{\ast}(t)+\gamma_{m}^{\alpha})-X_{[0,t_{f}-\tau_{1}]}t_{2m}^{\alpha}\lambda_{4}^{\ast}(t+\tau_{1})-\lambda_{4}^{\ast}(t)\epsilon_{1}u_{2}^{\ast}(t)\notag\\ &+\lambda_{4}^{\ast}(t)P_{3}^{\alpha}t_{2m}^{\alpha}-\lambda_{6}^{\ast}(t) (\frac{\beta_{m}^{\alpha}\alpha_{mm}^{\alpha}}{N}L_{m}^{\ast}(t)+ (1-\lambda_{m}^{\alpha})\beta_{m}^{\alpha}\frac{S^{\ast}(t)}{N}+\sigma_{m}^{\alpha}(1-\lambda_{m}^{\alpha})\beta_{m}^{\alpha} \frac{R^{\ast}(t)}{N}\notag\\ &+\alpha_{sm}^{\alpha} (1-\lambda_{m}^{\alpha})\beta_{m}^{\alpha}\frac{L_{s}^{\ast}(t)}{N}-(d^{\alpha} +\delta_{m}^{\alpha}+X_{[0,t_{f}-\tau_{1}]}t_{2m}^{\alpha}\lambda_{6}^{\ast} (t+\tau_{1})+\epsilon_{1}u_{1}^{\ast}(t)+\gamma_{m}^{\alpha}))\notag\\ &\left.\left.-\lambda_{8}^{\ast}(t)\left(P_{3}^{\alpha}t_{2m}^{\alpha}-\sigma_{m}^{\alpha}\beta_{m}^{\alpha}\frac{R^{\ast}}{N}\right)\right)\right), \end{align} (27) \begin{align} ^{c}_{t}D^{\alpha}_{t_{f}}\lambda_{7}^{\ast}=&- (-1+\lambda_{1}^{\ast}(t)\beta_{x}^{\alpha}\frac{S^{\ast}(t)}{N}+\lambda_{2}^{\ast}(t)\alpha_{sx}^{\alpha}\beta_{x}^{\alpha}\frac{L_{s}^{\ast}(t)}{N} +\lambda_{3}^{\ast}(t)\alpha_{mx}^{\alpha}\beta_{x}^{\alpha}\frac{L_{m}^{\ast}(t)}{N}-\lambda_{4}^{\ast}(t) (\lambda_{x}\sigma_{x}^{\alpha}\beta_{x}^{\alpha}\frac{R^{\ast}(t)}{N}\notag\\ &+\lambda_{x}^{\alpha}\beta_{x}^{\alpha}\frac{S^{\ast}(t)}{N}+\lambda_{x}^{\alpha}\alpha_{sx}^{\alpha}\beta_{x}^{\alpha}\frac{L_{s}^{\ast}(t)}{N}+\lambda_{x}^{\alpha}\alpha_{mx}^{\alpha}\beta_{x}^{\alpha}\frac{L_{m}^{\ast}(t)}{N}+\gamma_{x}^{\alpha}-\alpha_{xx}^{\alpha}\beta_{x}^{\alpha}\frac{L_{x}^{\ast}(t)}{N})\notag\\&-\lambda_{7}^{\ast}(t)(\alpha_{xx}^{\alpha}\beta_{x}^{\alpha}\frac{L_{x}^{\ast}(t)}{N}+(1-\lambda_{x}^{\alpha})\beta_{x}^{\alpha}\frac{S^{\ast}}{N}+(1-\lambda_{x}^{\alpha})\beta_{x}^{\alpha}\alpha_{sx}^{\alpha}\frac{L_{s}^{\ast}}{N}+(1-\lambda_{x}^{\alpha})\beta_{x}^{\alpha}\sigma_{x}^{\alpha}\frac{R^{\ast}}{N}\notag\\&+(1-\lambda_{x}^{\alpha})\beta_{x}^{\alpha}\alpha_{mx}^{\alpha}\frac{L_{m}^{\ast}}{N}+(d^{\alpha}+\delta_{x}^{\alpha}+\epsilon_{2}u_{2}^{\ast}+\gamma_{x}^{\alpha}))+X_{[0,t_{f}-\tau_{2}]}t_{2x}^{\alpha}\lambda_{8}^{\ast}(t+\tau_{2}) -X_{[0,t_{f}-\tau_{2}]}\notag\\&\times t_{2x}^{\alpha}\lambda_{7}^{\ast}(t+\tau_{2})-\lambda_{8}^{\ast}(t)(\epsilon_{2}u_{2}^{\ast}-\sigma_{x}^{\alpha}\beta_{x}^{\alpha}\frac{R}{N})),\\ \end{align} (28) \begin{align} ^{c}_{t}D^{\alpha}_{t_{f}}\lambda_{8}^{\ast}=& - \left( -\lambda_{2}^{\ast}(t)\beta_{s}^{\alpha}\lambda_{s}^{\alpha}\sigma_{s}^{\alpha}\frac{I_{s}^{\ast}(t)}{N} -\lambda_{3}^{\ast}(t)\beta_{m}^{\alpha}\lambda_{m}^{\alpha}\sigma_{m}^{\alpha}\frac{I_{m}^{\ast}(t)}{N} -\lambda_{4}^{\ast}(t)\beta_{x}^{\alpha}\lambda_{x}^{\alpha}\sigma_{x}^{\alpha}\frac{I_{x}^{\ast}(t)}{N}\right.\notag\\ &-\lambda_{5}^{\ast}(t)(1-\lambda_{s}^{\alpha})\beta_{s}^{\alpha}\sigma_{s}^{\alpha}\frac{I_{s}^{\ast}(t)}{N} -\lambda_{6}^{\ast}(t)(1-\lambda_{m}^{\alpha})\beta_{m}^{\alpha}\sigma_{m}^{\alpha}\frac{I_{m}^{\ast}(t)}{N}\notag\\ &\left.-\lambda_{7}^{\ast}(t)(1-\lambda_{x}^{\alpha})\beta_{x}^{\alpha}\sigma_{x}^{\alpha}\frac{I_{x}^{\ast}(t)}{N} +\lambda_{8}^{\ast}(t)\left(\sigma_{s}^{\alpha}\beta_{s}^{\alpha}\frac{I_{s}^{\ast}}{N} +\sigma_{m}^{\alpha}\beta_{m}^{\alpha}\frac{I_{m}^{\ast}}{N}+\sigma_{x}^{\alpha}\beta_{x}^{\alpha}\frac{I_{x}^{\ast}}{N}+d^{\alpha}\right)\right) \label{46}, \end{align} (29) (ii) with transversality conditions $$\lambda_{i}^{\ast}(T)=0$$, $$i = 1, . . . , 8$$. ${{X}_{[0,{{t}_{f}}-\tau ]}}=\left\{ \begin{array}{*{35}{l}} 1, & t\in \left[ 0,{{t}_{f}}-\tau \right], \\ 0, & otherwise. \\ \end{array} \right.$ (iii) optimality conditions: \begin{align} &H(S^{\ast}(t),L_{s}^{\ast}(t),L_{m}^{\ast}(t), L_{x}^{\ast}(t),I_{s}^{\ast}(t),I_{m}^{\ast}(t),I_{m}^{\ast}(t-\tau_{1}),I_{x}^{\ast}(t), I_{x}^{\ast}(t-\tau_{2}), R^{\ast}(t),\lambda^{\ast}(t), u_{k}^{\ast}(t) )= \notag\\ &\min _{0\leq u_{k}\leq 1}H(S^{\ast}(t),L_{s}^{\ast}(t),L_{m}^{\ast}(t),L_{x}^{\ast}(t),I_{s}^{\ast}(t),I_{m}^{\ast}(t), I_{m}^{\ast}(t-\tau_{1}),I_{x}^{\ast}(t), I_{x}^{\ast}(t-\tau_{2}) R^{\ast}(t),\lambda^{\ast}(t), u_{k}(t)),\label{7m}& \end{align} (30) \begin{align} u_{1}^{\ast}(t)&=\min\{\max\{0,\frac{\epsilon_{1}I_{m}^{\ast}(\lambda_{6}^{\ast}(t)-\lambda_{4}^{\ast}(t))}{B_{1}}\},1\}\label{10m},\\ \end{align} (31) \begin{align} u_{2}^{\ast}(t)&=\min\{\max\{0,\frac{\epsilon_{2}I_{x}^{\ast}(\lambda_{7}^{\ast}(t)-\lambda_{8}^{\ast}(t))}{B_{2}}\},1\}\label{11m}, \end{align} (32) where the stationarity condition is $$\frac{\partial H}{\partial u_{k}}=0,$$$$k=1,2.$$ Proof. Using conditions (18), we can get (22)–(29), where the Hamiltonian $$H$$ is given by: \begin{align} H&=H(S, L_{s}, L_{m}, L_{x}, I_{s},I_{m}, I_{x}, R, \lambda, u_{k} )\notag\\ &=I_{s}+I_{m}+I_{x}+L_{s}+u_{1}^{2}\frac{B_{1}}{2}+u_{2}^{2}\frac{B_{2}}{2}\notag\\ &\quad +\lambda_{1}(b^{\alpha}-d^{\alpha}S-\beta_{s}^{\alpha}\frac{SI_{s}}{N}-\beta_{m}^{\alpha}\frac{SI_{m}}{N}-\beta_{x}^{\alpha}\frac{SI_{x}}{N})\notag\\ &\quad +\lambda_{2}(\lambda_{s}^{\alpha}\beta_{s}^{\alpha}\frac{SI_{s}}{N}+\sigma_{s}^{\alpha}\lambda_{s}^{\alpha}\beta_{s}^{\alpha}\frac{RI_{s}}{N}-\alpha_{ss}^{\alpha}\beta_{s}^{\alpha}\frac{L_{s}I_{s}}{N}-\alpha_{sm}^{\alpha}\beta_{m}^{\alpha}\frac{L_{s}I_{m}}{N}-\alpha_{sx}^{\alpha}\beta_{x}^{\alpha}\frac{L_{s}I_{x}}{N}+\gamma_{s}^{\alpha}I_{s} \notag\\ &\quad -(d^{\alpha}+\epsilon_{s}^{\alpha}+t_{1s}^{\alpha})L_{s}\notag\\ &\quad +\lambda_{3}(\lambda_{m}^{\alpha}\beta_{m}^{\alpha}\frac{SI_{m}}{N}+\sigma_{m}^{\alpha}\lambda_{m}^{\alpha}\beta_{m}^{\alpha}\frac{RI_{m}}{N}+\alpha_{sm}^{\alpha}\beta_{m}^{\alpha}\lambda_{m}^{\alpha}\frac{L_{s}I_{m}}{N}-\alpha_{mm}^{\alpha}\beta_{m}^{\alpha}\frac{L_{m}I_{m}}{N}-\alpha_{mx}^{\alpha}\beta_{x}^{\alpha}\frac{L_{m}I_{x}}{N}\notag\\ &\quad +\gamma_{m}^{\alpha}I_{m}-(d^{\alpha}+\epsilon_{m}^{\alpha})L_{m}+t_{1s}^{\alpha}L_{s} -P_{1}^{\alpha}t_{1s}^{\alpha}L_{s}+t_{2s}^{\alpha}I_{s}-P_{2}^{\alpha}t_{2s}^{\alpha}I_{s},)\notag\\ &\quad +\lambda_{4}(\lambda_{x}^{\alpha}\beta_{x}^{\alpha}\frac{SI_{x}}{N}+\sigma_{x}^{\alpha}\lambda_{x}^{\alpha}\beta_{x}^{\alpha}\frac{RI_{x}}{N}+\alpha_{sx}^{\alpha}\beta_{x}^{\alpha}\lambda_{x}^{\alpha}\frac{L_{s}I_{x}}{N} +\alpha_{mx}^{\alpha}\beta_{x}^{\alpha}\lambda_{x}^{\alpha}\frac{L_{m}I_{x}}{N}-\alpha_{xx}^{\alpha}\beta_{x}^{\alpha}\frac{L_{x}I_{x}}{N} \notag\\ &\quad -(d^{\alpha}+\epsilon_{x}^{\alpha})L_{x}+\gamma_{x}^{\alpha}I_{x}+t_{2m}^{\alpha}I_{m}(t-\tau_{1})+\epsilon_{1}u_{1}(t)I_{m}-P_{3}^{\alpha}t_{2m}^{\alpha}I_{m})\notag\\ &\quad +\lambda_{5}(\alpha_{ss}^{\alpha}\beta_{s}^{\alpha}\frac{L_{s}I_{s}}{N}+(1-\lambda_{s}^{\alpha})\beta_{s}^{\alpha}(\frac{SI_{s}}{N}+\sigma_{s}^{\alpha}\frac{RI_{s}}{N})+\epsilon_{s}^{\alpha}L_{s}-(d^{\alpha}+\delta_{s}^{\alpha}+t_{2s}^{\alpha}+\gamma_{s}^{\alpha})\notag\\ &\quad +\lambda_{6}(\alpha_{mm}^{\alpha}\beta_{m}^{\alpha}\frac{L_{m}I_{m}}{N}+(1-\lambda_{m}^{\alpha})\beta_{m}^{\alpha}(\frac{SI_{m}}{N}+\sigma_{m}^{\alpha}\frac{RI_{m}}{N}+\alpha_{sm}^{\alpha}\frac{L_{s}I_{m}}{N})+\epsilon_{m}^{\alpha}L_{m} \notag\\ &\quad -(d^{\alpha}+\delta_{m}^{\alpha}+\epsilon_{1}u_{1}(t)+\gamma_{m}^{\alpha})I_{m}-t_{2m}^{\alpha}I_{m}(t-\tau_{1}))\notag\\ &\quad +\lambda_{7}(\alpha_{xx}^{\alpha}\beta_{x}^{\alpha}\frac{L_{x}I_{x}}{N} +(1-\lambda_{x}^{\alpha})\beta_{x}^{\alpha}(\frac{SI_{x}}{N}+\sigma_{x}^{\alpha}\frac{RI_{x}}{N} +\alpha_{sx}^{\alpha}\frac{L_{s}I_{x}}{N}+\alpha_{mx}^{\alpha}\frac{L_{m}I_{x}}{N})+\epsilon_{x}^{\alpha}L_{x} \notag\\ &\quad -(d^{\alpha}+\delta_{x}^{\alpha}+\gamma_{x}^{\alpha}+\epsilon_{2}u_{2}(t))I_{x}-t_{2x}^{\alpha}I_{x}(t-\tau_{2}))\notag\\ &\quad +\lambda_{8}(P_{1}^{\alpha}t_{1s}^{\alpha}L_{s}+P_{2}^{\alpha}t_{2s}^{\alpha}I_{s}+P_{3}^{\alpha}t_{2m}^{\alpha}I_{m}+t_{2x}^{\alpha}I_{x}+t_{2x}^{\alpha}I_{x}(t-\tau_{2})+\epsilon_{2}u_{2}(t)I_{x}-\sigma_{s}^{\alpha}\beta_{s}^{\alpha}\frac{RI_{s}}{N} \notag\\ &\quad -\sigma_{m}^{\alpha}\beta_{m}^{\alpha}\frac{RI_{m}}{N}-\sigma_{x}^{\alpha}\beta_{x}^{\alpha}\frac{RI_{x}}{N} -d^{\alpha}R)\label{6m}, \end{align} (33) and $$\lambda_{i},$$$$i=1,2,3,...,8$$ are the Lagrange multipliers. It is known as a co-state or adjoint variables. Moreover, the transversality conditions $$\lambda_{i}(T ) =0,$$$$i = 1, . . . , 8,$$ hold and the optimal controls (31)–(32) can be claimed from the minimization condition (30). □ 5. Numerical method 5.1. Generalized Euler method In the following, we will use the generalization of the classical Euler’s method, for solving the optimal control delay problem. For more details on this method, see Odibat & Momani (2008). The headlines of this method is given as follows, consider the problem (2)–(9) and (15). Assume that $$\tau_1$$ and $$\tau_2$$ are positive. The general formula for GEM when $t_{k+1}$$=$$t_{k}+h,$ for $$k = 0, 1, ..., Q-1,$$ is given as follows: \begin{align} S(t_{k+1})&=S(t_{k})+\frac{h^{\alpha}}{{\it{\Gamma}}(\alpha+1)}f_{1}(S(t_{k}),L_{s}(t_{k}),L_{m}(t_{k}),L_{x}(t_{k}),I_{s}(t_{k}),I_{m}(t_{k}) ,I_{m}(t_{\tau_{1}k}), I_{x}(t_{k}),\notag\\ &\quad \times I_{x}(t_{\tau_{2} k}),R(t_{k}), u_{1}, u_{2},t_{k}),\notag\\ L_{s}(t_{k+1})&=L_{s}(t_{k})+\frac{h^{\alpha}}{{\it{\Gamma}}(\alpha+1)}f_{2}(S(t_{k}),L_{s}(t_{k}),L_{m}(t_{k}),L_{x}(t_{k}),I_{s}(t_{k}),I_{m}(t_{k}),I_{m}(t_{\tau_{1} k}), I_{x}(t_{k}), \notag\\ &\quad \times I_{x}(t_{\tau_{2} k}),R(t_{k}),u_{1}, u_{2},t_{k}),\notag\\ L_{m}(t_{k+1})&=L_{m}(t_{k})+\frac{h^{\alpha}}{{\it{\Gamma}}(\alpha+1)}f_{3}(S(t_{k}),L_{s}(t_{k}),L_{m}(t_{k}),L_{x}(t_{k}),I_{s}(t_{k}), I_{m}(t_{k}),I_{m}(t_{\tau_{1} k}), I_{x}(t_{k}),\notag\\ &\quad \times I_{x}(t_{\tau_{2} k}),R(t_{k}), u_{1}, u_{2},t_{k}),\notag\\ L_{x}(t_{k+1})&=L_{x}(t_{k})+\frac{h^{\alpha}}{{\it{\Gamma}}(\alpha+1)}f_{4}(S(t_{k}),L_{s}(t_{k}),L_{m}(t_{k}),L_{x}(t_{k}),I_{s}(t_{k}) , I_{m}(t_{k}), I_{m}(t_{\tau_{1} k}), I_{x}(t_{k}),\notag\\ &\quad \times I_{x}(t_{\tau_{2} k}),R(t_{k}),u_{1}, u_{2},t_{k}),\notag\\ I_{s}(t_{k+1})&=I_{s}(t_{k})+\frac{h^{\alpha}}{{\it{\Gamma}}(\alpha+1)}f_{5}(S(t_{k}),L_{s}(t_{k}),L_{m}(t_{k}),L_{x}(t_{k}),I_{s}(t_{k}), I_{m}(t_{k}), I_{m}(t_{\tau_{1} k}), I_{x}(t_{k}),\notag\\ &\quad \times I_{x}(t_{\tau_{2} k}),R(t_{k}),u_{1}, u_{2},t_{k}),\notag\\ I_{m}(t_{k+1})&=I_{m}(t_{k})+\frac{h^{\alpha}}{{\it{\Gamma}}(\alpha+1)}f_{6}(S(t_{k}),L_{s}(t_{k}),L_{m}(t_{k}),L_{x}(t_{k}),I_{s}(t_{k}) , I_{m}(t_{k}),I_{m}(t_{\tau_{1} k}), I_{x}(t_{k}), \notag\\ &\quad \times I_{x}(t_{\tau_{2} k}),R(t_{k}),u_{1}, u_{2},t_{k}),\notag\\ I_{x}(t_{k+1})&=I_{x}(t_{k})+\frac{h^{\alpha}}{{\it{\Gamma}}(\alpha+1)}f_{7}(S(t_{k}),L_{s}(t_{k}),L_{m}(t_{k}),L_{x}(t_{k}),I_{s}(t_{k}) , I_{m}(t_{k}), I_{m}(t_{\tau_{1} k}), I_{x}(t_{k}), \notag\\& \quad \times I_{x}(t_{\tau_{2} k}),R(t_{k}),u_{1}, u_{2},t_{k}),\notag\\ R(t_{k+1})&=R(t_{k})+\frac{h^{\alpha}}{{\it{\Gamma}}(\alpha+1)}f_{8}(S(t_{k}),L_{s}(t_{k}),L_{m}(t_{k}),L_{x}(t_{k}),I_{s}(t_{j}), I_{m}(t_{k}),I_{m}(t_{\tau_{1} k}), I_{x}(t_{k}), \notag\\ &\quad \times I_{x}(t_{\tau_{2} k}),R(t_{k}),u_{1}, u_{2},t_{k}).\notag\ \end{align} It is clear that if $$\alpha = 1,$$ then the GEM reduces to the classical Euler’s method (Baleanu et al., 2012). To solve (22)–(29) with transversality conditions $$\lambda_{i}^{\ast}(T)=0$$, $$i = 1, . . . , 8$$ and $$\tau_{1}>0$$, $$\tau_{2}>0$$ by using backward Euler method. 5.2. Iterative optimal control method This method is a simple technique to study numerically fractional delay control problems, for more details see Sweilam & AL-Mekhlafi (2017), and the references cited therein. The main advantage of this method is to avoid the expected difficulties during computations such as integration by parts rule in case of fractional calculus. IOCM could adapt the same existence and characteristic optimal control theorems as in DDE systems for the system of fractional DDEs. This method can be described as follows: Algorithm 1: Step 0. Choose a starting point $$(S_{0}, L_{s0}, L_{m0},L_{x0},I_{s0},I_{m0},I_{x0},R_{0},)$$ and choose a tolerance $$\epsilon>0,$$ for the accuracy desired. Then set $$j = 1$$. Step 1. Let $$\tau_{1}$$, $$\tau_{2},$$ positive numberes, set $$u_{1}(t)$$ and $$u_{2}(t)$$ as given in (31) and (32) respectively. Solve the co-state system (22)–(29) with transversality conditions $$\lambda_{i}(T)=0,$$$$j = 1,2, 3,...,8,$$ then go to Step 3. Step 2. Plugging the values of $$u_{1}(t)$$ and $$u_{2}(t)$$ into the system (2)–(9), then solve it with the same starting point to obtain the new starting point $$(S(t)$$, $$L_{s}(t),$$$$L_{m}(t)$$,$$L_{x}(t)$$,$$I_{s}(t)$$, $$I_{m}(t)$$,$$I_{x}(t)$$, $$R(t))$$. Step 3. The stopping criterion is as follows: if $$|u_{j} -u_{j + 1}|< \epsilon$$ stop, otherwise put $$j = j + 1$$, and go to Step 1. 5.3. Discretizations of IOCM To solve optimization problem (2)–(9) and (13), first we need to discretized the state system (2)–(9). The Grünwald–Letnikov’s definition is used to discretization the fractional derivative $$^{c}_{a}D^{\alpha}_{t}f(t)$$. Then the system (2)–(9), can be written in the following form: \begin{align} S^{n}=&\frac{b^{\alpha}-\sum^{n}_{j=1}q^{\alpha}_{j}S^{n-j}}{q^{-\alpha}_{0}+d^{\alpha}+\frac{\beta_{s}^{\alpha}I_{s}^{n}+\beta_{m}^{\alpha}I_{m}^{n}+\beta_{x}^{\alpha}I_{x}^{n}}{N^{n}}},\\ \end{align} (34) \begin{align} L_{s}^{n}=&\frac{\frac{\beta_{s}^{\alpha}I_{s}^{n}}{N^{n}}\lambda_{s}^{\alpha}(S^{n}+\sigma_{s}^{\alpha}R^{n})+\gamma_{s}^{\alpha}I_{s}^{n}-\sum^{n}_{j=1}q^{\alpha}_{j}L_{s}^{n-j}}{q^{-\alpha}_{0}+(d^{\alpha}+t_{1s}^{\alpha}+\epsilon_{s}^{\alpha})+\frac{1}{N^{n}}(\alpha_{ss}^{\alpha}\beta_{s}^{\alpha}I_{s}^{n}+\alpha_{sm}^{\alpha}\beta_{m}^{\alpha}I_{m}^{n}+\alpha_{sx}^{\alpha}\beta_{x}^{\alpha}I_{x}^n)},\\ \end{align} (35) \begin{align} L_{m}^{n}=&\frac{\frac{\beta_{m}^{\alpha}\lambda_{m}^{\alpha}I_{m}^{n}}{N^{n}}(S^{n}+\sigma_{m}^{\alpha}R^{n}+\alpha_{sm}^{\alpha}L_{s}^{n})+\gamma_{m}^{\alpha}I_{m}^{n}+t_{1s}^{\alpha}L_{s}^{n}(1-P_{1}^{\alpha})}{q^{-\alpha}_{0}+(d^{\alpha}+\epsilon_{m}^{\alpha})+\frac{1}{N^{n}}(\alpha_{mm}^{\alpha}\beta_{m}^{\alpha}I_{m}^{n}+\alpha_{mx}^{\alpha}\beta_{x}^{\alpha}I_{x}^{n})} \notag\\ &+\frac{t_{2s}^{\alpha}I_{s}^{n}(1-P_{2}^{\alpha})-\sum^{n}_{j=1}q^{\alpha}_{j}L_{m}^{n-j}}{\omega^{-\alpha}_{0}+(d^{\alpha}+\epsilon_{m}^{\alpha})+\frac{1}{N^{n}}(\alpha_{mm}^{\alpha}\beta_{m}^{\alpha}I_{m}^{n}+\alpha_{mx}^{\alpha}\beta_{x}^{\alpha}I_{x}^{n})},\\ \end{align} (36) \begin{align} L_{x}^{n}=&\frac{\frac{\beta_{x}^{\alpha}\lambda_{x}^{\alpha}I_{x}^{n}}{N^{n}}(S^{n}+\sigma_{x}^{\alpha}R^{n}+\alpha_{sx}^{\alpha}L_{s}^{n}+\alpha_{mx}^{\alpha}L_{m}^{n})+t_{2m}^{\alpha}I_{m}^{n\tau_{1}}+\epsilon_{1}u_{1}^{n}I_{m}^{n}-P_{3}^{\alpha}t_{2m}^{\alpha}I_{m}^{n}}{q^{-\alpha}_{0}+(d^{\alpha}+\epsilon_{x}^{\alpha})+\frac{1}{N^{n}}(\alpha_{xx}^{\alpha}\beta_{x}^{\alpha}I_{x}^{n})} \notag\\ &+\frac{\gamma_{x}^{\alpha}I_{x}^{n}-\sum^{n}_{j=1}q^{\alpha}_{j}L_{x}^{n-j}}{q^{-\alpha}_{0}+(d^{\alpha}+\epsilon_{x}^{\alpha})+\frac{1}{N^{n}}(\alpha_{xx}^{\alpha}\beta_{x}^{\alpha}I_{x}^{n})},\\ \end{align} (37) \begin{align} I_{s}^{n}=&\frac{\varphi_{5}(h)\beta_{s}^{\alpha}\frac{I_{s}^{n}}{N^{n}}(\alpha_{ss}^{\alpha}L_{s}^{n}+(1-\lambda_{s}^{\alpha})(S^{n}+\sigma_{s}^{\alpha}R^{n}))}{q^{-\alpha}_{0}+(d^{\alpha}+\delta_{s}^{\alpha}+\gamma_{s}^{\alpha}+t_{2s}^{\alpha}} \notag\\ &+\frac{\epsilon_{s}^{\alpha}L_{s}^{n}-\sum^{n}_{j=1}q^{\alpha}_{j}I_{s}^{n-j}}{q^{-\alpha}_{0}+(d^{\alpha}+\delta_{s}^{\alpha}+\gamma_{s}^{\alpha}+t_{2s}^{\alpha})},\\ \end{align} (38) \begin{align} I_{m}^{n}=&\frac{\beta_{m}^{\alpha}\frac{I_{m}^{n}}{N^{n}}(\alpha_{mm}^{\alpha}L_{m}^{n}+(1-\lambda_{m}^{\alpha})(S^{n}+\sigma_{m}^{\alpha}R^{n}+\alpha_{sm}^{\alpha}L_{s}^{n}))}{q^{-\alpha}_{0}+(d^{\alpha}+\delta_{m}^{\alpha}+\gamma_{m}^{\alpha}+\epsilon_{1}u_{1}^{n}))+t_{2m}^{\alpha}I_{m}^{n\tau_{1}}} \notag\\ &+\frac{\epsilon_{m}^{\alpha}L_{m}^{n}-\sum^{n}_{j=1}q^{\alpha}_{j}I_{m}^{n-j}}{q^{-\alpha}_{0}+(d^{\alpha}+\delta_{m}^{\alpha}+\gamma_{m}^{\alpha}+\epsilon_{1}u_{1}^{n})+t_{2m}^{\alpha}I_{m}^{n\tau_{1}}},\\ \end{align} (39) \begin{align} I_{x}^{n}=&\frac{\beta_{x}^{\alpha}\frac{I_{x}^{n}}{N^{n}}(\alpha_{xx}^{\alpha}L_{x}^{n}+(1-\lambda_{x}^{\alpha})(S^{n}+\sigma_{x}^{\alpha}R^{n}+\alpha_{sx}^{\alpha}L_{s}^{n}+\alpha_{mx}^{\alpha}L_{m}^{n}))}{q^{-\alpha}_{0}+(d^{\alpha}+\delta_{x}^{\alpha}+\gamma_{x}^{\alpha}+\epsilon_{2}u_{2}^{n})+t_{2x}^{\alpha}I_{x}^{n\tau_{1}}} \notag\\ &+\frac{\epsilon_{x}L_{x}^{n}-\sum^{n}_{j=1}q^{\alpha}_{j}I_{x}^{n-j}}{q^{-\alpha}_{0}+(d^{\alpha}+\delta_{x}^{\alpha}+\gamma_{x}^{\alpha}+\epsilon_{2}u_{2}^{n})+t_{2x}^{\alpha}I_{x}^{n\tau_{2}}},\\ \end{align} (40) \begin{align} R^{n}=&\frac{t_{1s}^{\alpha}P_{1}^{\alpha}L_{s}^{n}+P_{2}^{\alpha}t_{2s}^{\alpha}I_{s}^{n}+t_{2m}^{\alpha}P_{3}^{\alpha}I_{m}^{n}+t_{2x}^{\alpha}I_{x}^{n\tau_{2}}+\epsilon_{2}u_{2}^{n}I_{x}^{n}-\sum^{n}_{j=1}q^{\alpha}_{j}R^{n-j}}{q^{-\alpha}_{0}+d^{\alpha}+\frac{1}{N^{n}}(\sigma_{s}^{\alpha}\beta_{s}^{\alpha}I_{s}^{n}+\sigma_{m}^{\alpha}\beta_{m}^{\alpha}I_{m}^{n}+\sigma_{x}^{\alpha}\beta_{x}^{\alpha}I_{x}^{n})}. \end{align} (41) To solve (22)–(29) with transversality conditions $$\lambda_{i}^{\ast}(T)=0$$, $$i = 1, . . . , 8,$$ we will use the backward Euler method. 6. Numerical results and discussions In the following, IOCM and GEM are used to obtain the approximate solutions for systems (2)–(9) and (22)–(29). Using the initial conditions ($$S(0)$$, $$L_{s}(0)$$, $$L_{m}(0)$$, $$L_{x}(0)$$, $$I_{s}(0),$$$$I_{m}(0),$$$$I_{x}(0),$$$$R(0))$$$$=$$ ($$\frac{76}{120}N,$$$$\frac{20}{120}N,$$$$\frac{5}{120}N,$$$$\frac{2}{120}N,$$$$\frac{8}{120}N,$$$$\frac{4}{120}N,$$$$\frac{2}{120}N,$$$$\frac{3}{120}N$$), and the parameters given in Table 2. The approximate solutions of the proposed system are given in Figs 1–9 at different values of $$\tau_{1},$$$$\tau_{2}$$ and $$0<\alpha\leq1$$. In Fig. 1, the effect of different time delays $$\tau_{1},$$ with fixed time delay $$\tau_{2},$$ on the behavior the approximate solutions for state variables $$I_{m}(t-\tau_{1}),$$$$I_{x}(t-\tau_{2})$$ in the controlled case. The delay $$\tau_{1}$$ was chosen to equal $$1, 2$$, and $$3$$. The delay $$\tau _2$$ is chosen in this case fixed and equal $$2$$. Moreover, we comper the obtianed results with the non-delay case, i.e., $$\tau_{1}=\tau_{2}=0$$. We noted that, the result of $$I_{x}(t-\tau_{2})$$ is not affected by changing the value of $$\tau_{1}$$. Also, Fig. 2, shows the effect of different time delays $$\tau_{2},$$ with fixed $$\tau_{1},$$ on the behavior of the approximate solutions for state variables $$I_{m}(t-\tau_{1}),$$$$I_{x}(t-\tau_{2})$$ in the controlled case. The delay $$\tau_{2}$$ was chosen to equal $$1, 2, 3,$$ and the delay $$\tau_{1}=1.$$ Moreover, the results are compared with non-delay case. We noted that, the results of $$I_{m}(t-\tau_{1})$$ is not affected by changing the value of $$\tau_{2}$$. Table 3, shows the values of the objective functional in both controlled and uncontrolled cases by using IOCM with $$\tau_{1}=0.5$$, $$\tau_{2}=0.2$$. Moreover, if we compute the value of the objective functional with $$\tau_{1}=0.2$$, $$\tau_{2}=0.5$$, we obtian the same result of the previous case. Fig. 3, shows that, the quantity $$S(t)+L_{s}(t)+L_{m}(t)+L_{x}(t)+I_{s}(t)+I_{m}(t)+I_{x}(t)+R(t)/N$$ is a constant in time at $$\alpha=1,$$ in controlled case, $$0\leq u_{1},u_{2}\leq 1$$. Fig. 4, shows the behavior of the approximate solutions for state variables $$I_{m},$$$$I_{x}$$ and $$R(t)$$ with $$\alpha=1$$ and $$\tau_{1}=\tau_{2}=\tau=0.1,$$ in the controlled and uncontrolled cases by using IOCM and GEM. We noted that, in the uncontrolled case, the delays in diagnosis and commencement of treatment to the individuals of $$I_{m}$$ and $$I_{x}$$ causing a shortage in the number of individuals of $$R(t)$$, but in the controlled case the number of $$I_{m}$$ and $$I_{x}$$ are decreasing with time as we expected and the number of $$R(t)$$ is increasing. Fig. 5, shows the behavior of the approximate solutions for state variables $$I_{m},$$$$I_{x}$$ with $$\alpha=0.97$$ and $$\tau=0.2$$ by using GEM and IOCM. We noted that, the results obtained by IOCM is better than the results obtained by GEM, where the number of $$I_{m}$$ and $$I_{x}$$ are decreasing with time in a rate faster than the results obtained by GEM. Fig. 6, shows the optimal controls $$u_{1}^{\ast}$$ and $$u_{2}^{\ast},$$ for $$\alpha=0.99$$, in a time units of years and $$\tau_{1}=\tau_{2}=\tau$$ by using ICOM. We noted that, when the value of $$\tau$$ is increasing the value of the objective functional is decreasing. Fig. 7, shows that the relationship between $$I_{m}(t)$$, $$I_{m}(t-\tau)$$ and $$I_{x}(t)$$, $$I_{x}(t-\tau)$$ in controlled case, where $$\tau=0.5$$, $$\alpha=0.99$$. We noted, that the number of $$I_{m}(t-\tau)$$ and $$I_{x}(t-\tau)$$ are more than the number of $$I_{m}(t)$$ and $$I_{x}(t)$$. Fig. 8, shows the behavior of the approximate solutions of the state variables $$I_{m}(t)$$, $$I_{x}(t)$$ and $$R(t)$$ with different value of $$\alpha$$, which are given to demonstrate how the fractional model is a generalization of the integer order model. Fig. 9, shows the values of $$u_{1}^{\ast}, u_{2}^{\ast}$$, in a time units of years by using ICOM with different values of $$\alpha$$ and in controlled case. Tables 4 and 5 show competitive studies between the two proposed methods to find the optimum values of the objective functional in both controlled and uncontrolled cases at different values of $$\alpha$$ and $$\tau_{1}=\tau_{2}=\tau$$. From the numerical results in these tables, it is found that IOCM can claim the minimum values in a rate faster than GEM and we can conclude that it is better than GEM, because, GEM can be failed to derive the optimum values for some values of $$\alpha$$, as given in Table 4. All results were obtained by using MATLAB (R2013a). on a computer machine with intel(R) core $$i3-3110M$$$$@$$$$2.40 GHz$$ and $$4 GB$$ RAM. Fig. 1. View largeDownload slide Plot of values of $$I_{m}(t-\tau_{1}), I_{x}(t-\tau_{2})$$ in a time units of years by using IOCM at different value of $$\tau_{1}$$ and fixed $$\tau_{2}=1$$, $$\alpha=0.99$$. Fig. 1. View largeDownload slide Plot of values of $$I_{m}(t-\tau_{1}), I_{x}(t-\tau_{2})$$ in a time units of years by using IOCM at different value of $$\tau_{1}$$ and fixed $$\tau_{2}=1$$, $$\alpha=0.99$$. Fig. 2. View largeDownload slide Plot of values of $$I_{m}(t-\tau_{1}), I_{x}(t-\tau_{2})$$ in a time units of years by using IOCM at different value of $$\tau_{2}$$ and fixed $$\tau_{1}=1$$, $$\alpha=0.99$$. Fig. 2. View largeDownload slide Plot of values of $$I_{m}(t-\tau_{1}), I_{x}(t-\tau_{2})$$ in a time units of years by using IOCM at different value of $$\tau_{2}$$ and fixed $$\tau_{1}=1$$, $$\alpha=0.99$$. Fig. 3. View largeDownload slide Plot of $$S(t)+L_{s}(t)+L_{m}(t)+L_{x}(t)+I_{s}(t)+I_{m}(t)+I_{x}(t)+R(t)/N$$ versus $$t$$ in years and $$\alpha=1$$, $$\tau=0.1,$$ using GEM. Fig. 3. View largeDownload slide Plot of $$S(t)+L_{s}(t)+L_{m}(t)+L_{x}(t)+I_{s}(t)+I_{m}(t)+I_{x}(t)+R(t)/N$$ versus $$t$$ in years and $$\alpha=1$$, $$\tau=0.1,$$ using GEM. Fig. 4. View largeDownload slide Numerical simulations of the state variable $$I_{m}$$, $$I_{x}$$ and $$R(t)$$ for the controlled case when $$0\leq u_{1}, u_{2}\leq 1$$, compared with the uncontrolled case when $$u_{1}= u_{2}=0$$ and $$\alpha=1$$, $$\tau=0.1,$$ using GEM and IOCM. Fig. 4. View largeDownload slide Numerical simulations of the state variable $$I_{m}$$, $$I_{x}$$ and $$R(t)$$ for the controlled case when $$0\leq u_{1}, u_{2}\leq 1$$, compared with the uncontrolled case when $$u_{1}= u_{2}=0$$ and $$\alpha=1$$, $$\tau=0.1,$$ using GEM and IOCM. Fig. 5. View largeDownload slide Numerical simulations of the state variables $$I_{m}$$, $$I_{x}$$ for the controlled case when $$0\leq u_{1}, u_{2}\leq 1$$, compared with the uncontrolled case when $$u_{1}= u_{2}=0$$ and $$\alpha=0.97$$, $$\tau=0.2,$$ using GEM and IOCM. Fig. 5. View largeDownload slide Numerical simulations of the state variables $$I_{m}$$, $$I_{x}$$ for the controlled case when $$0\leq u_{1}, u_{2}\leq 1$$, compared with the uncontrolled case when $$u_{1}= u_{2}=0$$ and $$\alpha=0.97$$, $$\tau=0.2,$$ using GEM and IOCM. Fig. 6. View largeDownload slide Plot of values of $$u_{1}^{\ast}, u_{2}^{\ast}$$ for $$\alpha=0.99$$, in a time units of years by using IOCM and different $$\tau$$. Fig. 6. View largeDownload slide Plot of values of $$u_{1}^{\ast}, u_{2}^{\ast}$$ for $$\alpha=0.99$$, in a time units of years by using IOCM and different $$\tau$$. Fig. 7. View largeDownload slide The relationship between $$I_{m}(t)$$, $$I_{m}(t-\tau)$$ and $$I_{x}(t)$$, $$I_{x}(t-\tau)$$ in controlled case, where $$\tau=0.5$$, $$\alpha=0.99$$. Fig. 7. View largeDownload slide The relationship between $$I_{m}(t)$$, $$I_{m}(t-\tau)$$ and $$I_{x}(t)$$, $$I_{x}(t-\tau)$$ in controlled case, where $$\tau=0.5$$, $$\alpha=0.99$$. Fig. 8. View largeDownload slide Numerical simulations of the state variable $$I_{m}, I_{x}$$ and $$R$$ for the controlled case with different values of $$\alpha$$, $$\tau=0.1$$, using IOCM. Fig. 8. View largeDownload slide Numerical simulations of the state variable $$I_{m}, I_{x}$$ and $$R$$ for the controlled case with different values of $$\alpha$$, $$\tau=0.1$$, using IOCM. Fig. 9. View largeDownload slide Plot of values of $$u_{1}^{\ast}, u_{2}^{\ast}$$ in a time units of years by using IOCM at different $$\alpha$$ and $$\tau=0.1$$. Fig. 9. View largeDownload slide Plot of values of $$u_{1}^{\ast}, u_{2}^{\ast}$$ in a time units of years by using IOCM at different $$\alpha$$ and $$\tau=0.1$$. Table 3 The values of objective functional in both controlled and uncontrolled cases by using IOCM with $$\tau_{1}=0.5$$, $$\tau_{2}=0.2$$ $$\alpha$$ $$J(u_{1}^{\ast},u_{2}^{\ast})$$ with control $$J(u_{1}^{\ast},u_{2}^{\ast})$$ without control 1 35801 53734 0.98 27062 49283 0.90 10392 35762 0.80 3184.1 26817.3 0.70 1054.2 17813.6 0.60 4580.2 10888.2 0.50 230.5 6898.6 0.40 118.4 5065.9 $$\alpha$$ $$J(u_{1}^{\ast},u_{2}^{\ast})$$ with control $$J(u_{1}^{\ast},u_{2}^{\ast})$$ without control 1 35801 53734 0.98 27062 49283 0.90 10392 35762 0.80 3184.1 26817.3 0.70 1054.2 17813.6 0.60 4580.2 10888.2 0.50 230.5 6898.6 0.40 118.4 5065.9 Table 3 The values of objective functional in both controlled and uncontrolled cases by using IOCM with $$\tau_{1}=0.5$$, $$\tau_{2}=0.2$$ $$\alpha$$ $$J(u_{1}^{\ast},u_{2}^{\ast})$$ with control $$J(u_{1}^{\ast},u_{2}^{\ast})$$ without control 1 35801 53734 0.98 27062 49283 0.90 10392 35762 0.80 3184.1 26817.3 0.70 1054.2 17813.6 0.60 4580.2 10888.2 0.50 230.5 6898.6 0.40 118.4 5065.9 $$\alpha$$ $$J(u_{1}^{\ast},u_{2}^{\ast})$$ with control $$J(u_{1}^{\ast},u_{2}^{\ast})$$ without control 1 35801 53734 0.98 27062 49283 0.90 10392 35762 0.80 3184.1 26817.3 0.70 1054.2 17813.6 0.60 4580.2 10888.2 0.50 230.5 6898.6 0.40 118.4 5065.9 Table 4 The values of objective functional in both controlled and uncontrolled cases with $$\tau_{1}=\tau_{2}=0.3$$ $$\alpha$$ Methods $$J(u_{1}^{\ast},u_{2}^{\ast})$$ with control $$J(u_{1}^{\ast},u_{2}^{\ast})$$ without control 1 GEM 33496 56227 IOCM 33496 0.98 GEM 27165 45941 IOCM 23243 0.90 GEM 10150 15402 IOCM 7269 0.80 GEM 5146.1 3179.8 IOCM 2911.3 0.70 GEM Non 995.87 IOCM 966.13 0.60 GEM Non 448.002 IOCM 380.555 0.50 GEM Non 232.32 IOCM 198.1596 $$\alpha$$ Methods $$J(u_{1}^{\ast},u_{2}^{\ast})$$ with control $$J(u_{1}^{\ast},u_{2}^{\ast})$$ without control 1 GEM 33496 56227 IOCM 33496 0.98 GEM 27165 45941 IOCM 23243 0.90 GEM 10150 15402 IOCM 7269 0.80 GEM 5146.1 3179.8 IOCM 2911.3 0.70 GEM Non 995.87 IOCM 966.13 0.60 GEM Non 448.002 IOCM 380.555 0.50 GEM Non 232.32 IOCM 198.1596 Table 4 The values of objective functional in both controlled and uncontrolled cases with $$\tau_{1}=\tau_{2}=0.3$$ $$\alpha$$ Methods $$J(u_{1}^{\ast},u_{2}^{\ast})$$ with control $$J(u_{1}^{\ast},u_{2}^{\ast})$$ without control 1 GEM 33496 56227 IOCM 33496 0.98 GEM 27165 45941 IOCM 23243 0.90 GEM 10150 15402 IOCM 7269 0.80 GEM 5146.1 3179.8 IOCM 2911.3 0.70 GEM Non 995.87 IOCM 966.13 0.60 GEM Non 448.002 IOCM 380.555 0.50 GEM Non 232.32 IOCM 198.1596 $$\alpha$$ Methods $$J(u_{1}^{\ast},u_{2}^{\ast})$$ with control $$J(u_{1}^{\ast},u_{2}^{\ast})$$ without control 1 GEM 33496 56227 IOCM 33496 0.98 GEM 27165 45941 IOCM 23243 0.90 GEM 10150 15402 IOCM 7269 0.80 GEM 5146.1 3179.8 IOCM 2911.3 0.70 GEM Non 995.87 IOCM 966.13 0.60 GEM Non 448.002 IOCM 380.555 0.50 GEM Non 232.32 IOCM 198.1596 Table 5 The values of objective functional in both controlled and uncontrolled cases with $$\tau_{1}=\tau_{2}=0$$ $$\alpha$$ Method $$J(u_{1}^{\ast},u_{2}^{\ast})$$ with control $$J(u_{1}^{\ast},u_{2}^{\ast})$$ without control 1 GEM 36626 57704 IOCM 36626 0.98 GEM 30614 48921 IOCM 25441 0.95 GEM 22791 36854 IOCM 15463 0.90 GEM 13383 21322 IOCM 6672.1 $$\alpha$$ Method $$J(u_{1}^{\ast},u_{2}^{\ast})$$ with control $$J(u_{1}^{\ast},u_{2}^{\ast})$$ without control 1 GEM 36626 57704 IOCM 36626 0.98 GEM 30614 48921 IOCM 25441 0.95 GEM 22791 36854 IOCM 15463 0.90 GEM 13383 21322 IOCM 6672.1 Table 5 The values of objective functional in both controlled and uncontrolled cases with $$\tau_{1}=\tau_{2}=0$$ $$\alpha$$ Method $$J(u_{1}^{\ast},u_{2}^{\ast})$$ with control $$J(u_{1}^{\ast},u_{2}^{\ast})$$ without control 1 GEM 36626 57704 IOCM 36626 0.98 GEM 30614 48921 IOCM 25441 0.95 GEM 22791 36854 IOCM 15463 0.90 GEM 13383 21322 IOCM 6672.1 $$\alpha$$ Method $$J(u_{1}^{\ast},u_{2}^{\ast})$$ with control $$J(u_{1}^{\ast},u_{2}^{\ast})$$ without control 1 GEM 36626 57704 IOCM 36626 0.98 GEM 30614 48921 IOCM 25441 0.95 GEM 22791 36854 IOCM 15463 0.90 GEM 13383 21322 IOCM 6672.1 Conclusion In this article, numerical solutions for the optimal control of fractional multi-strain TB model with two time delays are presented. Modified parameters are introduced to account for the fractional order model. We introduce two delays in the proposed model, representing the time delay in the diagnosis and commencement of treatment of individuals with active TB infection in two and three strains. Two control functions $$u_{1},$$$$u_{2}$$ are introduced. The controls $$u_{1},$$$$u_{2}$$ represent the effort that prevents the failure of treatment in active TB infectious individuals $$I_{m},$$$$I_{x},$$ e.g., supervising the patients, helping them to take the TB medications regularly and to complete the TB treatment. Necessary optimality conditions are derived. Two numerical methods; IOCM and GEM are introduced to study the behavior of the model problem. Comparative studies are implemented. Some figures are given to demonstrate how the fractional delay model is a generalization of the integer order model. It can be concluded from the numerical results presented in this article that, IOCM is better than GEM and it can be applied simply and effectively to solve optimal control of fractional delay multi-strain TB. Acknowledgements The authors would like to thanks the anonymous reviewers very much for their positive comments, careful reading, and useful suggestions on improving this article. References Abta, A. Laarabi, H. & Alaoui, H. T. ( 2014 ) The Hopf bifurcation analysis and optimal control of a delayed SIR epidemic model. Int. J. Anal. , 2014 , 1 – 10 . Google Scholar CrossRef Search ADS Agarwal, O. P. ( 2008 ) A formulation and numerical scheme for fractional optimal control problems. J. Vib. Control , 14 , 1291 – 1299 . Google Scholar CrossRef Search ADS Agarwal, P. & Choi, J. ( 2016 ) Fractional calculus operators and their image formulas. J. Korean Math. Soc. , 53 , 1183 – 1210 . Google Scholar CrossRef Search ADS Aldila, D. , Gãtz, T. & Soewono, E. ( 2013 ) An optimal control problem arising from a dengue disease transmission model. Math. Biosci ., 242 , 9 – 16 . Google Scholar CrossRef Search ADS PubMed Allen, L. J. S. ( 2007 ) An Introduction to Mathematical Biology . NJ : Prentice Hall . Arenas, A. J. , Gonzãlez-Parra, G. & Chen-Charpentierc, B. M. ( 2016 ) Construction of nonstandard finite difference schemes for the SI and SIR epidemic models of fractional order. Math. Comput. Simul. , 121 , 48 – 63 . Google Scholar CrossRef Search ADS Arino, J. & Soliman, I. A. ( 2015 ) A model for the spread of tuberculosis with drug-sensitive and emerging multidrug-resistant and extensively drug resistant strains. Math. Comput. Model. , Wiley, 1 – 120 , doi:10.1002/9781118853887.ch5 . Baleanu, D. , Diethelm, K. , Scalas, E. & Trujillo, J. J. ( 2012 ) Fractional Calculus Models and Numerical Methods. Series on Complexity, Nonlinearity and Chaos . Hackensack, NJ : World Scientific Publishing Co. Pte. Ltd . Google Scholar CrossRef Search ADS Bellen, A. & Zennaro, M. ( 2003 ) Numerical Methods for Delay Differential Equations, Oxford Science Publications . Oxford, UK : Clarendon Press . Google Scholar CrossRef Search ADS Bhalekar, S. & Daftardar-Gejji, V. ( 2011 ) A predictor-corrector scheme for solving non-linear delay differential equations of fractional order. J. Fract. Calc. Appl. , 5 , 1 – 9 . Bocharov, G. & A Rihan, F. ( 2000 ) Numerical modelling in biosciences using delay differential equations. J. Comput. Appl. Math. , 125 , 183 – 199 . Google Scholar CrossRef Search ADS Boklund, A. , Toft, N. , Alban, L. & Uttenthal, A. ( 2009 ) Comparing the epidemiological and economic effects of control strategies against classical swine fever in Denmark. Prev. Vet. Med. , 90 , 180 – 193 . Google Scholar CrossRef Search ADS PubMed Castillo-ćhavez, C. & Feng, Z. ( 1997 ) To treat or not to treat: the case of tuberculosis. J. Math. Biol. , 35 , 629 – 656 . Google Scholar CrossRef Search ADS PubMed Cohen, T. , & Murray, M. ( 2004 ) Modeling epidemics of multidrug-resistant M. tuberculosis of heterogeneous fitness. Nat. Med. , 10 , 1117 – 1121 . Google Scholar CrossRef Search ADS PubMed Cooke, R. & Bellman, K. L. ( 1963 ) Differential-Difference Equations . New York, NY : Academic Press. Google Scholar CrossRef Search ADS Denysiuk, R. , Silva, C. J. & Torres, D. F. M. ( 2014 ) Multiobjective approach to optimal control for a tuberculosis model. Optim. Meth. Software, 2015 , 1029 – 4937 . Driver, R. ( 1977 ) Ordinary and Delay Differential Equations . Berlin, Germany : Springer . Google Scholar CrossRef Search ADS El-Saka, H. A. A. ( 2013 ) The fractional-order SIR and SIRs epidemic models with variable population size. Math. Sci. Lett., 2 , 1 – 6 . Google Scholar CrossRef Search ADS Fowler, A. C. , Mackey, M. C. ( 2002 ) Relaxation oscillations in a class of delay differential equations. SIAM J. Appl. Math. 63 : 299 – 323 , (2002). Google Scholar CrossRef Search ADS Hale, J. ( 1977 ) Theory of Functional Differential Equations . New York, NY, USA : Springer . Google Scholar CrossRef Search ADS Hattaf, K. , Yousfi, N. ( 2012 ) Optimal control of a delayed HIV infection model with immune response using an efficient numerical method. Int. Sch. Res. Netw. , 2012 , 1 – 7 . Kuang, Y. ( 1993 ) Delay Differential Equations with Applications in Population Dynamics . San Diego : Academic Press . Laarabi, H. , Abta, A. & Hattaf, K. ( 2015 ) Optimal control of a delayed sirs epidemic model with vaccination and treatment, Acta Biotheor. , 63 , 87 – 97 . Google Scholar CrossRef Search ADS PubMed Matignon, D. ( 1996 ) Stability results for fractional differential equations with applications to control processing, Computational Engineering in Systems and Application. Multiconference , Vol. 2. Lille, France : IMACS, IEEE-SMC , pp. 963 – 968 . Nagy, A. M. , Sweilam, N. H. ( 2014 ) An efficient method for solving fractional Hodgkin Huxley model. Phys. Lett. A, 378 , 1980 – 1984 . Google Scholar CrossRef Search ADS Nelson, P. W. , Perelson, A. S. ( 2002 ) Mathematical analysis of delay differential equation models of HIV-1 infection. Math. Biosci. , 179 , 73 – 94 . Google Scholar CrossRef Search ADS PubMed Odibat, Z. & Momani, S. ( 2008 ) An algorithm for the numerical solution of differential equations of fractional order. J. Appl. Math. Inform ., 26 , 15 – 27 . Pimenov, V. G. & Tashirova, E. E. ( 2013 ) Numerical methods for solving a hereditary equation of hyperbolic type. Proc. Steklov Inst. Math. , 281 , 126 – 136 . Google Scholar CrossRef Search ADS Podlubny, I. ( 1999 ) Fractional Differential Equations . New York : Academic Press . Rihanm, F. A. ( 2000 ) Numerical Treatment of Delay Dfferential Equation in Bioscience . Manchester, UK : University of Manchester . Rihan, F. A. , Hashish, A. , Al-Maskari, F. , Hussein, M. S. and Ahmed, E. ( 2016 ) Dynamics of tumor-immune system with fractional-order. J. Tumor Res. 2 , 109 . Rihan, F. A. , Lakshmanan, S. , Hashish, A. H. , Rakkiyappan, R. & Ahmed, E. ( 2015 ) Fractional-order delayed predator-prey systems with Holling type-II functional response, Nonlinear Dynamics. An International. J. Nonlinear Dynam. Chaos Eng. Syst. , 80 , 777 – 789 . Google Scholar CrossRef Search ADS Salahshour, S. , Ahmadian, A. , Senu, N. , Baleanu, D. & Agarwal, P. ( 2015 ) On analytical solutions of the fractional differential equation with uncertainty: application to the basset problem. Entropy , 17 , 885 – 902 . Google Scholar CrossRef Search ADS Small, P. M. and Fujiwara, P. I. ( 2001 ) Management of Tuberculosis in the United States. N. Engl. J. Med. , 345 , 189 – 200 . Google Scholar CrossRef Search ADS PubMed Smith, H. ( 2011 ) An Introduction To Delay Differential Equations with Applications To The Life Sciences . New York : Springer . Google Scholar CrossRef Search ADS Silva, C. J. , Maurer, H. & Torres, D. F. ( 2017 ) Optimal control of a tuberculosis model with state and control delays. Math. Biosci. Eng. , 14 . Sreeramareddy, C. T. , Panduru, K. V. , Menten, J. & Van den Ende, J. ( 2009 ) Time delays in diagnosis of pulmonary tuberculosis: a systematic review of literature. BMC Infect. Dis. , 9 . Styblo, K. ( 1978 ) State of art: epidemiology of tuberculosis. Bull. Int. Union Tuberc., 53 , 141 – 152 . Sweilam, N. H. , AL-Mekhlafi, S. M. ( 2016 ) Comparative study for multi-strain tuberculosis (TB) model of fractional order. J. Appl. Math. Inf. Sci. , 10 , 1403 – 1413 . Google Scholar CrossRef Search ADS Sweilam, N. H. , AL-Mekhlafi, S. M. ( 2017 ) Legendre spectral-collocation method for solving fractional optimal control of HIV infection of $$Cd4^{+}T$$ cells mathematical model. J. Defense Model. Simulat. Appl. Methodol. Technol. , 14 , 273 – 284 . Google Scholar CrossRef Search ADS Sweilam, N. H. & AL-Mekhlafi, S. M. ( 2016a ) Numerical study for multi-strain tuberculosis (TB) model of variable-order fractional derivatives. J. Adv. Res. , 7 , 271 – 283 . Google Scholar CrossRef Search ADS Sweilam, N. H. , AL-Mekhlafi, S. M. ( 2016b ) On the optimal control for fractional multi-strain TB model. Optim. Contr. Appl. Meth. , 37 , 1355 – 1374 . Google Scholar CrossRef Search ADS Sweilam, N. H. , Khader, M. M. & Adel, M. ( 2012a ) Numerical studies for fractional-order logistic differential equation with two different delays. J. Appl. Math. , 14 pages, Article ID 764894, doi:10.1155/2012/764894 . Sweilam, N. H. , Khader, M. M. & Adel, M. ( 2012b ) On the stability analysis of weighted average finite difference methods for fractional wave equation. Fract. Differ. Calc. , 2 , 17 – 29 . Google Scholar CrossRef Search ADS Sweilam, N. H. , Khader, M. M. & Nagy, A. M. ( 2011 ) Numerical solution of two-sided space-fractional wave equation using finite difference method. J. Comput. Appl. Math. 235 , 2832 – 2841 . Google Scholar CrossRef Search ADS Sweilam, N. H. , Soliman, I. A. & Al-Mekhlafi, S. M. ( 2017 ) Nonstandard finite difference method for solving the multi-strain TB model. J. Egypt. Math. Soc. , 25 , 129 – 138 . Google Scholar CrossRef Search ADS Tariboon, J. , Ntouyas, S. K. & Agarwal, P. ( 2015 ) New concepts of fractional quantum calculus and applications to impulsive fractional q-difference equations. Adv. Differ. Equ. , 2015 . Toman, K. ( 1979 ) Tuberculosis Case-Finding and Chemotherapy: Questions and Answers . Geneva : WHO . Uys, P. W. , Warren, M. & van Helden, P. D. ( 2007 ) A threshold value for the time delay to TB diagnosis. PLoS One , 2 . Van den Driessche, P. & Watmough, J. ( 2002 ) Reproduction numbers and subthreshold endemic equilibria for compartmental models of disease transmission. Math. Biosc., 180 , 29 – 48 . Google Scholar CrossRef Search ADS Wang, Z. , Huang, X. & Zhou, J. ( 2013 ) A numerical method for delayed fractional order differential equations: based on G-L Definition. Appl. Math. Inf. Sci. , 7 , 525 – 529 . Google Scholar CrossRef Search ADS World Health Organization ( 2012 ) Multidrug and Extensively Drug-Resistant TB (M/XDR-TB): 2012 Global Report on Surveillance and Response . World Health Organization : Geneva, Switzerland . World Health Organization ( 2014 ) Global Tuberculosis Report 2014 . Geneva : World Health Organization , http://www.who.int/tb/publications/global_report/en/. Zhang, X. , Agarwal, P. , Liu, Z. , Ding, W. & Ciancio, A. ( 2016 ) On the fractional differential equations with not instantaneous impulses. Open Phys., 14 , 676 – 684 . Zhou, H. , Yang, L. , Agarwal, P. ( 2017 ) Solvability for fractional p-Laplacian differential equations with multipoint boundary conditions at resonance on infinite interval. J. Appl. Math. Comput , 53 , 51 – 76 . Google Scholar CrossRef Search ADS © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png IMA Journal of Mathematical Control and Information Oxford University Press

# Optimal control for a time delay multi-strain tuberculosis fractional model: a numerical approach

, Volume Advance Article – Nov 9, 2017
24 pages

/lp/ou_press/optimal-control-for-a-time-delay-multi-strain-tuberculosis-fractional-WnO7SNRLjH
Publisher
Oxford University Press
ISSN
0265-0754
eISSN
1471-6887
D.O.I.
10.1093/imamci/dnx046
Publisher site
See Article on Publisher Site

### Abstract

Abstract In this article, optimal control for a novel fractional multi-strain Tuberculosis model with time delay is presented. The proposed model is governed by a system of fractional delay differential equations, where the fractional derivative is defined in the Caputo sense. Modified parameters are introduced to account for the fractional order. Two delays in the proposed model representing the time delay on the diagnosis and commencement of treatment of individuals with active Tuberculosis infection in two and three strains. Necessary and sufficient conditions that guarantee the existence and the uniqueness of the solution of the resulting systems are given. Two control variables are proposed to minimize the cost of interventions. Two simple numerical methods are used to study the nonlinear fractional delay optimal control problem. The methods are the iterative optimal control method (IOCM) and the generalized Euler method (GEM). Comparative studies are implemented, it is found that the IOCM is better than the GEM. 1. Introduction Mathematical modelling with delay differential equations (DDEs) is widely used for analysis and predictions in epidemiology, immunology and physiology (Bocharov & Rihan, 2000; Rihan, 2000; Fowler & Mackey, 2002; Nelson & Perelson, 2002; Smith, 2011). The time delays in these models can be considered as dependence of the present state of the modelled system on its history. The general theory of differential equations with delays is discussed in the literature (Cooke & Bellman, 1963; Driver, 1977; Hale, 1977; Van den Driessche & Watmough, 2002; Bellen & Zennaro, 2003; Pimenov & Tashirova, 2013). In reality, some issues are rarely so instantaneous; but there is usually a propagation delay before the effects are felt. Recently, many optimal control models pertaining to epidemic diseases have appeared in the literature. They include, but not limited to, delayed SIRS epidemic model (Laarabi et al., 2015, where the authors proposed an optimal control model governed by a delayed SIRS model with a single time delay representing the incubation period), delayed SIR model (Abta et al., 2014, where the authors formulated an optimal control model to optimize the costs of vaccination and treatment on a delayed SIR model with a single delay representing the incubation period and a saturated incidence rate) and tuberculosis model (Silva et al., 2017), where the authors proposed optimal control of a Tuberculosis model with state and control delays using data from Angola). Sweilam & AL-Mekhlafi (2016b), studied optimal control of a general nonlinear multi-strain tuberculosis (TB) model. Four controls are used for preventing the failure of treatment in active TB. In Sweilam & AL-Mekhlafi (2017), we have also performed numerical studies of optimal control of HIV: HIV model (Hattaf & Yousfi, 2012, where the authors used optimal controls to represent the efficiency of drug treatment in inhibiting viral production and preventing new infections), swine flu (Boklund et al., 2009, where the authors studied epidemiological and economic consequences of some control strategies in a classical swine flu epidemic under Danish conditions with respect to herd demographics and geography and investigated the effect of extra biosecurity measures on farms and dengue fever (Aldila et al., 2013, where the authors designed an optimal control problem using four different control parameters and discussed some results about epidemic prevention). Tuberculosis (TB) can be considered as one of the most important infectious diseases, it is the second largest cause of mortality by infectious diseases and is a challenging disease to control (World Health Organization, 2014). Delays to the diagnosis of active TB present a major obstacle to the control of a TB epidemic (Uys et al., 2007); it may worsen the disease, increase the risk of death and enhance tuberculosis transmission to the community (Toman, 1979; Sreeramareddy et al., 2009). Both patient and the health system may be responsible for the diagnosis delay (Sreeramareddy et al., 2009). On the other hand, mathematical models are quite important and efficient tool to describe and investigate TB epidemics, for more details see (Styblo, 1978; Castillo-ćhavez & Feng, 1997; Cohen & Murray, 2004; Small & Fujiwara, 2001; Denysiuk et al., 2014). In this work, we consider a general model of multi-strain TB diseases with time-delay memory. We introduce a discrete time-delay in the state variable of active TB infection of two and three strain, that represents the delay in the diagnosis and commencement of treatment. The multi-strain TB model which incorporates three strains: drug-sensitive, emerging multi-drug resistant (MDR) and extensively drug-resistant (XDR) is developed by Arino & Soliman (2015). This model included several factors of spreading TB such as the fast infection, the exogenous reinfection and secondary infection along with the resistance factor. Sweilam and AL$$-$$Mekhlafi introduced some numerical studies of this model in (Sweilam & AL-Mekhlafi, 2016,a,b,c; Sweilam et al., 2017). Fractional differential equations have been the focus of many studies due to their frequent appearance in various sciences (Sweilam et al., 2011, 2012a,b; Nagy & Sweilam, 2014; Salahshour et al., 2015; Tariboon et al., 2015; Agarwal & Choi, 2016; Rihan et al., 2016; Sweilam & AL-Mekhlafi, 2016,a,b, 2017; Zhang et al., 2016; Sweilam et al., 2017; Zhou et al., 2017). Delayed fractional differential equations are used to describe dynamical systems (Allen, 2007; Rihan et al., 2016). Recently, DFDEs begin to raises the attention of many researchers (Bhalekar & Daftardar-Gejji, 2011; Wang et al., 2013; Rihan et al., 2015). Simulating these equations are an important technique in the research, accordingly, finding effective numerical methods for the DFDEs are necessary process. In this article, we introduce the fractional order multi-strain TB with memory time-delay and modified parameters, this modified parameters are introduced to account for the fractional order (Arenas et al., 2016). The aim of this article is to study optimal control of the proposed model for minimization of the number of active TB infectious of two and three strains. Two simple numerical methods are used to study the nonlinear fractional delay optimal control problem (FDOCP). The methods are the iterative optimal control method (IOCM) and the generalized Euler method (GEM). Comparative studies are implemented. This article is organized as follows: In Section 2, mathematical preliminaries of fractional calculus which are required for establishing the results are given. In Section 3, the proposed system is described by fractional order derivatives with memory time-delay is presented and the stability of equilibrium points is studied. Formulation of the optimal control problem and the necessary optimality conditions for the fractional order multi-strain TB model are derived in Section 4. Numerical methods are introduced for solving the proposed control model in Section 5. Numerical simulations are given in Section 6. In Section 6, the conclusions are given. 2. Definitions and preliminaries Let us consider the signal fractional differential equation (see Podlubny, 1999; Nagy & Sweilam, 2014), $$D_{t}^{\alpha}z(t)=g(t,z(t)),\,\,\,T \geq t \geq 0,\,\, and,\,z(t_{0})=0. \label{ss}$$ (1) Definition 2.1 Podlubny (1999) The Grünwald–Letnikovs fractional derivative is defined as: \begin{align} D^{\alpha}_{t}g(t)&=\lim_{h\rightarrow0}\frac{1}{h^{\alpha}}\sum_{i=0}^{[\frac{t}{h}]}q^{\alpha}_{i}f(t-ih),\notag \end{align} where $$[\frac{t}{h}]$$ denotes the integer part of $$\frac{t}{h}$$ and \begin{align} q^{\alpha}_{i}&=(-1)^{i}(^{\alpha}_{i})=\frac{{\it{\Gamma}}(\alpha+1)}{i!{\it{\Gamma}}(\alpha-i+1)},\notag \end{align} with $$(^{\alpha}_{i})$$ being the fractional binomial coefficients. It is clear that the coefficients $$q^{\alpha}_{i}$$ can be evaluated recursively as follows : $$q^{\alpha}_{0}=1,$$$$q^{\alpha}_{i}=(1-\frac{\alpha+1}{i})q^{\alpha}_{i-1},$$$$i\geq 1.$$ Definition 2.2 Podlubny (1999) The left Caputo fractional derivative is defined as follows: \begin{align} ^{c}_{a}D^{\alpha}_{t}g(t)&=\frac{1}{{\it{\Gamma}}(n-\alpha)}\int^{t}_{a}(t-\tau)^{n-\alpha-1}\frac{d^{n}}{d\tau^{n}}z(\tau)d\tau,\notag \end{align} the right Caputo fractional derivative is defined as follows : \begin{align} ^{c}_{a}D^{\alpha}_{t}g(t)&=\frac{(-1)^{n}}{{\it{\Gamma}}(n-\alpha)}\int^{b}_{t}(t-\tau)^{n-\alpha-1}\frac{d^{n}}{d\tau^{n}}z(\tau)d\tau.\notag \end{align} Where $${\it{\Gamma}}$$ is the Euler gamma function. 3. Mathematical model In the following, the control problem of the fractional model of multi-strain TB with time-delay memory and modified parameters are presented. This model incorporates three strains: drug-sensitive, MDR and XDR. The population is divided into eight compartments depending on their epidemiological stages as follows: susceptible $$(S)$$; latently infected with drug sensitive TB $$(L_{s})$$; latently infected with MDR-TB $$(L_{m})$$; latently infected with XDR-TB $$(L_{x})$$; sensitive drug TB infectious $$(I_{s})$$; MDR-TB infectious $$(I_{m})$$; XDR-TB infectious $$(I_{x})$$; recovered $$R$$. The new parameters of the model are described in Table 1, this modified parameters are introduced to account for the fractional order (Arenas et al., 2016). We introduce a discrete time-delay in the state variables $$I_{m}$$ and $$I_{x},$$ denoted by $$\tau_{1}$$ and $$\tau_{2}$$, these represents the delay to diagnosis and commencement of treatment of active TB infection of two and three strain. One of the main assumptions of this model is that, the total population, $$N;$$$$N = S(t)+L_{s}(t)+L_{m}(t)+L_{x}(t)+I_{s}(t)+I_{m}(t)+I_{x}(t)+R(t),$$ is a constant in time. In other words, we assume that the birth and death rates are equal and there are no TB-related deaths i.e., $$\delta_{s}^{\alpha}=\delta_{m}^{\alpha}=\delta_{x}^{\alpha}$$=0. Table 1 All adapted parameters and their interpretation of the system (2)–(9) Parameters Interpretation $$d^{\alpha}$$ Death and birth rate $$\lambda_{j}^{\alpha}$$ Rate of infected individuals move to $$L_{j}$$ with strain $$j \in\{s ,m,x\}\$$ $$1-\lambda_{j}^{\alpha}$$ Rate of newly infected individuals progressing to active TB with strain $$j \in\{s,m,x \}$$ $$\beta_{j}^{\alpha}$$ Transmission coefficient with strain $$j \in\{s,m,x \}$$ $$\varepsilon_{j}^{\alpha}$$ Rate of endogenous reactivation of $$L_{j}$$ $$\gamma_{j}^{\alpha}$$ Rate of natural recovery to the latent stage $$L_{j}$$ $$\delta_{j}^{\alpha}$$ Rate of death due to $$TB$$ of strain $$j$$ $$\alpha_{j1}^{\alpha},\alpha_{j2}^{\alpha}$$ Rate of exogenous reinfection of $$L_{j1}$$ due to contact with $$I_{j2}$$ $$1-\sigma_{j}^{\alpha}$$ Efficiency of treatment in preventing infection with strain $$j$$ $$P_{1}^{\alpha}$$ Probability of treatment success for $$L_{s}$$ $$1-P_{1}^{\alpha}$$ Proportion of treated $$L_{s}$$ moved to $$L_{m}$$ due to incomplete treatment or lack of strict compliance in the use of drugs $$P_{2}^{\alpha}$$ Probability of treatment success for $$I_{s}$$ $$1-P_{2}^{\alpha}$$ Proportion of treated $$I_{s}$$ moved to $$L_{m}$$ due to incomplete treatment or lack of strict compliance in the use of drugs $$P_{3}^{\alpha}$$ Probability of treatment success for $$I_{m}$$ $$1-P_{3}^{\alpha}$$ proportion of treated $$I_{m}$$ moved to $$L_{x}$$due to incomplete treatment or lack of strict compliance in the use of drugs $$t_{1s}^{\alpha}$$ Rate of treatment for $$L_{s}$$ $$t_{2j}^{\alpha}$$ Rate of treatment for $$I_{j}$$. Note that $$t_{2x}$$ is the rate of successful treatment of $$I_{x}, j\in\{x,m,s\}$$ Parameters Interpretation $$d^{\alpha}$$ Death and birth rate $$\lambda_{j}^{\alpha}$$ Rate of infected individuals move to $$L_{j}$$ with strain $$j \in\{s ,m,x\}\$$ $$1-\lambda_{j}^{\alpha}$$ Rate of newly infected individuals progressing to active TB with strain $$j \in\{s,m,x \}$$ $$\beta_{j}^{\alpha}$$ Transmission coefficient with strain $$j \in\{s,m,x \}$$ $$\varepsilon_{j}^{\alpha}$$ Rate of endogenous reactivation of $$L_{j}$$ $$\gamma_{j}^{\alpha}$$ Rate of natural recovery to the latent stage $$L_{j}$$ $$\delta_{j}^{\alpha}$$ Rate of death due to $$TB$$ of strain $$j$$ $$\alpha_{j1}^{\alpha},\alpha_{j2}^{\alpha}$$ Rate of exogenous reinfection of $$L_{j1}$$ due to contact with $$I_{j2}$$ $$1-\sigma_{j}^{\alpha}$$ Efficiency of treatment in preventing infection with strain $$j$$ $$P_{1}^{\alpha}$$ Probability of treatment success for $$L_{s}$$ $$1-P_{1}^{\alpha}$$ Proportion of treated $$L_{s}$$ moved to $$L_{m}$$ due to incomplete treatment or lack of strict compliance in the use of drugs $$P_{2}^{\alpha}$$ Probability of treatment success for $$I_{s}$$ $$1-P_{2}^{\alpha}$$ Proportion of treated $$I_{s}$$ moved to $$L_{m}$$ due to incomplete treatment or lack of strict compliance in the use of drugs $$P_{3}^{\alpha}$$ Probability of treatment success for $$I_{m}$$ $$1-P_{3}^{\alpha}$$ proportion of treated $$I_{m}$$ moved to $$L_{x}$$due to incomplete treatment or lack of strict compliance in the use of drugs $$t_{1s}^{\alpha}$$ Rate of treatment for $$L_{s}$$ $$t_{2j}^{\alpha}$$ Rate of treatment for $$I_{j}$$. Note that $$t_{2x}$$ is the rate of successful treatment of $$I_{x}, j\in\{x,m,s\}$$ Table 1 All adapted parameters and their interpretation of the system (2)–(9) Parameters Interpretation $$d^{\alpha}$$ Death and birth rate $$\lambda_{j}^{\alpha}$$ Rate of infected individuals move to $$L_{j}$$ with strain $$j \in\{s ,m,x\}\$$ $$1-\lambda_{j}^{\alpha}$$ Rate of newly infected individuals progressing to active TB with strain $$j \in\{s,m,x \}$$ $$\beta_{j}^{\alpha}$$ Transmission coefficient with strain $$j \in\{s,m,x \}$$ $$\varepsilon_{j}^{\alpha}$$ Rate of endogenous reactivation of $$L_{j}$$ $$\gamma_{j}^{\alpha}$$ Rate of natural recovery to the latent stage $$L_{j}$$ $$\delta_{j}^{\alpha}$$ Rate of death due to $$TB$$ of strain $$j$$ $$\alpha_{j1}^{\alpha},\alpha_{j2}^{\alpha}$$ Rate of exogenous reinfection of $$L_{j1}$$ due to contact with $$I_{j2}$$ $$1-\sigma_{j}^{\alpha}$$ Efficiency of treatment in preventing infection with strain $$j$$ $$P_{1}^{\alpha}$$ Probability of treatment success for $$L_{s}$$ $$1-P_{1}^{\alpha}$$ Proportion of treated $$L_{s}$$ moved to $$L_{m}$$ due to incomplete treatment or lack of strict compliance in the use of drugs $$P_{2}^{\alpha}$$ Probability of treatment success for $$I_{s}$$ $$1-P_{2}^{\alpha}$$ Proportion of treated $$I_{s}$$ moved to $$L_{m}$$ due to incomplete treatment or lack of strict compliance in the use of drugs $$P_{3}^{\alpha}$$ Probability of treatment success for $$I_{m}$$ $$1-P_{3}^{\alpha}$$ proportion of treated $$I_{m}$$ moved to $$L_{x}$$due to incomplete treatment or lack of strict compliance in the use of drugs $$t_{1s}^{\alpha}$$ Rate of treatment for $$L_{s}$$ $$t_{2j}^{\alpha}$$ Rate of treatment for $$I_{j}$$. Note that $$t_{2x}$$ is the rate of successful treatment of $$I_{x}, j\in\{x,m,s\}$$ Parameters Interpretation $$d^{\alpha}$$ Death and birth rate $$\lambda_{j}^{\alpha}$$ Rate of infected individuals move to $$L_{j}$$ with strain $$j \in\{s ,m,x\}\$$ $$1-\lambda_{j}^{\alpha}$$ Rate of newly infected individuals progressing to active TB with strain $$j \in\{s,m,x \}$$ $$\beta_{j}^{\alpha}$$ Transmission coefficient with strain $$j \in\{s,m,x \}$$ $$\varepsilon_{j}^{\alpha}$$ Rate of endogenous reactivation of $$L_{j}$$ $$\gamma_{j}^{\alpha}$$ Rate of natural recovery to the latent stage $$L_{j}$$ $$\delta_{j}^{\alpha}$$ Rate of death due to $$TB$$ of strain $$j$$ $$\alpha_{j1}^{\alpha},\alpha_{j2}^{\alpha}$$ Rate of exogenous reinfection of $$L_{j1}$$ due to contact with $$I_{j2}$$ $$1-\sigma_{j}^{\alpha}$$ Efficiency of treatment in preventing infection with strain $$j$$ $$P_{1}^{\alpha}$$ Probability of treatment success for $$L_{s}$$ $$1-P_{1}^{\alpha}$$ Proportion of treated $$L_{s}$$ moved to $$L_{m}$$ due to incomplete treatment or lack of strict compliance in the use of drugs $$P_{2}^{\alpha}$$ Probability of treatment success for $$I_{s}$$ $$1-P_{2}^{\alpha}$$ Proportion of treated $$I_{s}$$ moved to $$L_{m}$$ due to incomplete treatment or lack of strict compliance in the use of drugs $$P_{3}^{\alpha}$$ Probability of treatment success for $$I_{m}$$ $$1-P_{3}^{\alpha}$$ proportion of treated $$I_{m}$$ moved to $$L_{x}$$due to incomplete treatment or lack of strict compliance in the use of drugs $$t_{1s}^{\alpha}$$ Rate of treatment for $$L_{s}$$ $$t_{2j}^{\alpha}$$ Rate of treatment for $$I_{j}$$. Note that $$t_{2x}$$ is the rate of successful treatment of $$I_{x}, j\in\{x,m,s\}$$ Two control functions $$u_{1}(\cdot),$$$$u_{2}(\cdot),$$ and two real positive constants $$\epsilon_{1},$$$$\epsilon_{2},$$ will be added to the model. The controls $$u_{1},$$$$u_{2}$$ represent the effort in preventing the failure of treatment in active TB infectious individuals $$I_{m},$$$$I_{x},$$ respectively, e.g., supervising the patients, helping them to take the TB medications regularly and to complete the TB treatment. The parameters $$\epsilon_{i}$$$$\in$$$$(0, 1)$$, $$i = 1, 2,$$ measure the effectiveness of the controls $$u_{k},$$$$k = 1, 2,$$ respectively, i.e., these parameters measure the efficacy of treatment interventions for active TB infectious individuals $$I_{m},$$$$I_{x}$$. The resulting model is given by a system of nonlinear fractional order differential equations with time-delay memory as follows: \begin{align} ^{c}_{a}D^{\alpha}_{t}S=&\,d^{\alpha}N-d^{\alpha}S-\beta_{s}^{\alpha}\frac{SI_{s}}{N}-\beta_{m}^{\alpha}\frac{SI_{m}}{N}-\beta_{x}^{\alpha}\frac{SI_{x}}{N}, \label{1} \\ \end{align} (2) \begin{align} ^{c}_{a}D^{\alpha}_{t}L_{s}=&\,\lambda_{s}^{\alpha}\beta_{s}^{\alpha}\frac{SI_{s}}{N}+\sigma_{s}^{\alpha}\lambda_{s}^{\alpha}\beta_{s}^{\alpha}\frac{RI_{s}}{N}-\alpha_{ss}^{\alpha}\beta_{s}^{\alpha}\frac{L_{s}I_{s}}{N}-\alpha_{sm}^{\alpha}\beta_{m}^{\alpha}\frac{L_{s}I_{m}}{N}-\alpha_{sx}^{\alpha}\beta_{x}^{\alpha}\frac{L_{s}I_{x}}{N}+\gamma_{s}^{\alpha}I_{s}\notag\\ &-(d^{\alpha}+\epsilon_{s}^{\alpha}+t_{1s}^{\alpha})L_{s},\\ \end{align} (3) \begin{align} ^{c}_{a}D^{\alpha}_{t}L_{m}=&\,\lambda_{m}^{\alpha}\beta_{m}^{\alpha}\frac{SI_{m}}{N}+\sigma_{m}^{\alpha}\lambda_{m}^{\alpha}\beta_{m}^{\alpha}\frac{RI_{m}}{N}+\alpha_{sm}^{\alpha}\beta_{m}^{\alpha}\lambda_{m}^{\alpha}\frac{L_{s}I_{m}}{N}-\alpha_{mm}^{\alpha}\beta_{m}^{\alpha}\frac{L_{m}I_{m}}{N}-\alpha_{mx}^{\alpha}\beta_{x}^{\alpha}\frac{L_{m}I_{x}}{N}\notag\\ &+\gamma_{m}^{\alpha}I_{m}-(d^{\alpha}+\epsilon_{m}^{\alpha})L_{m}+t_{1s}^{\alpha}L_{s}-P_{1}^{\alpha}t_{1s}^{\alpha}L_{s}+t_{2s}^{\alpha}I_{s}%\notag\\& -P_{2}^{\alpha}t_{2s}^{\alpha}I_{s},\\ \end{align} (4) \begin{align} ^{c}_{a}D^{\alpha}_{t}L_{x}=&\,\lambda_{x}^{\alpha}\beta_{x}^{\alpha}\frac{SI_{x}}{N}+\sigma_{x}^{\alpha}\lambda_{x}^{\alpha}\beta_{x}^{\alpha}\frac{RI_{x}}{N}+\alpha_{sx}^{\alpha}\beta_{x}^{\alpha}\lambda_{x}^{\alpha}\frac{L_{s}I_{x}}{N} +\alpha_{mx}^{\alpha}\beta_{x}^{\alpha}\lambda_{x}^{\alpha}\frac{L_{m}I_{x}}{N}-\alpha_{xx}^{\alpha}\beta_{x}^{\alpha}\frac{L_{x}I_{x}}{N}\notag\\ &-(d^{\alpha}+\epsilon_{x}^{\alpha})L_{x}+\gamma_{x}^{\alpha}I_{x}+t_{2m}^{\alpha}I_{m}(t-\tau_{1})+\epsilon_{1}u_{1}(t)I_{m}-P_{3}^{\alpha}t_{2m}^{\alpha}I_{m},\\ \end{align} (5) \begin{align} ^{c}_{a}D^{\alpha}_{t}I_{s}=&\,\alpha_{ss}^{\alpha}\beta_{s}^{\alpha}\frac{L_{s}I_{s}}{N}+(1-\lambda_{s}^{\alpha})\beta_{s}^{\alpha} \left(\frac{SI_{s}}{N}+\sigma_{s}^{\alpha}\frac{RI_{s}}{N}\right)+\epsilon_{s}^{\alpha}L_{s}-(d^{\alpha}+\delta_{s}^{\alpha}+t_{2s}^{\alpha}+\gamma_{s}^{\alpha}),\\ \end{align} (6) \begin{align} ^{c}_{a}D^{\alpha}_{t}I_{m}=&\,\alpha_{mm}^{\alpha}\beta_{m}^{\alpha}\frac{L_{m}I_{m}}{N}+(1-\lambda_{m}^{\alpha})\beta_{m}^{\alpha} \left(\frac{SI_{m}}{N}+\sigma_{m}^{\alpha}\frac{RI_{m}}{N}+\alpha_{sm}^{\alpha}\frac{L_{s}I_{m}}{N}\right)+\epsilon_{m}^{\alpha}L_{m} \notag\\& -(d^{\alpha}+\delta_{m}^{\alpha}+\epsilon_{1}u_{1}(t)+\gamma_{m}^{\alpha})I_{m}-t_{2m}^{\alpha}I_{m}(t-\tau_{1}),\\ \end{align} (7) \begin{align} ^{c}_{a}D^{\alpha}_{t}I_{x}=&\,\alpha_{xx}^{\alpha}\beta_{x}^{\alpha}\frac{L_{x}I_{x}}{N}+(1-\lambda_{x}^{\alpha})\beta_{x}^{\alpha} \left(\frac{SI_{x}}{N}+\sigma_{x}^{\alpha}\frac{RI_{x}}{N}+\alpha_{sx}^{\alpha}\frac{L_{s}I_{x}}{N}+\alpha_{mx}^{\alpha}\frac{L_{m}I_{x}}{N}\right)+\epsilon_{x}^{\alpha}L_{x} \notag\\&-(d^{\alpha}+\delta_{x}^{\alpha}+\gamma_{x}^{\alpha}+\epsilon_{2}u_{2}(t))I_{x}-t_{2x}^{\alpha}I_{x}(t-\tau_{2}),\\ \end{align} (8) \begin{align} ^{c}_{a}D^{\alpha}_{t}R=&\,P_{1}^{\alpha}t_{1s}^{\alpha}L_{s}+P_{2}^{\alpha}t_{2s}^{\alpha}I_{s}+P_{3}^{\alpha}t_{2m}^{\alpha}I_{m}+t_{2x}^{\alpha}I_{x}(t-\tau_{2})+\epsilon_{2}u_{2}(t)I_{x}-\sigma_{s}^{\alpha}\beta_{s}^{\alpha}\frac{RI_{s}}{N} \notag\\&-\sigma_{m}^{\alpha}\beta_{m}^{\alpha}\frac{RI_{m}}{N}-\sigma_{x}^{\alpha}\beta_{x}^{\alpha}\frac{RI_{x}}{N} -d^{\alpha}R \label{8}. \end{align} (9) Where $$^{c}_{a}D^{\alpha}_{t}$$ is the Caputo fractional derivative (Sweilam & AL-Mekhlafi, 2016b). Because the model (2)–(9) monitors the dynamics of human populations, all the parameters are assumed to be nonnegative. The initial conditions for system (2)–(9) are: $$S(\vartheta)=\theta_{1}(\vartheta),$$$$L_{s}(\vartheta)=\theta_{2}(\vartheta),$$$$L_{m}(\vartheta)=\theta_{3}(\vartheta),$$$$L_{x}(\vartheta)=\theta_{4}(\vartheta),$$$$I_{s}(\vartheta)=\theta_{5}(\vartheta),$$$$I_{m}(\vartheta)=\theta_{6}(\vartheta),$$$$I_{x}(\vartheta)=\theta_{7}(\vartheta),$$$$R(\vartheta)=\theta_{8}(\vartheta),$$$$\vartheta\in[-\tau,0],$$ where, $$\theta=(\theta_{1},\theta_{2}, \cdot\cdot\cdot, \theta_{8},)^{T}$$$$\in$$$$C,$$ where $$C$$ is the Banach space $$C([0,\tau],\mathbb{R}^{8})$$. From biological meaning, we further assume that $$\theta_{i}>0$$ for $$i = 1, \cdot\cdot\cdot, 8.$$ Throughout this article, we focus on the dynamics of the solutions of (2)–(9) in the restricted region: $${\it{\Omega}}=\{(S,L_{s},L_{m},L_{x},I_{s},I_{m},I_{x},R)\in \mathbb{R}^{8}|S+L_{s}+L_{m}+L_{x}+I_{s}+I_{m}+I_{x}+R=N\}.$$ It is worth mentioning that, the usual local existence, uniqueness and continuation results hold in this region as prove in (Hale, 1977; Kuang, 1993). Hence, a unique solution $$(S(t)$$,$$L_{s}(t)$$, $$L_{m}(t)$$, $$L_{x}(t),$$$$I_{s}(t)$$, $$I_{m}(t),$$$$I_{x}(t),$$$$R(t))$$ of (2)–(9) with the given initial conditions exists for all time $$t \geq 0$$. Consider the solutions of (2)–(9) with $$(S,L_{s},L_{m},L_{x},I_{s},I_{m},I_{x},R)$$$$\in$$$${\it{\Omega}}'$$, where $${\it{\Omega}}'$$, is the interior of $${\it{\Omega}},$$ for all $$\vartheta\in[-\tau,0]$$. Then the solutions stay in the interior of the region for all time $$t\geq 0,$$ i.e., the region is positively invariant with respect to system (2)–(9) (see e.g., Kuang, 1993). It is known that, a mathematical model has a disease free equilibrium if it has an equilibrium point at which the population remains in the absence of the disease (Van den Driessche & Watmough, 2002). The model (2)–(9) has a disease free equilibrium given by $$E_{0}=(N, 0,0,0,0,0,0,0)$$. The basic reproduction number ($$R_{0}$$) represents the expected average number of new TB infections produced by a single TB active infected individual when in contact with a completely susceptible population (Van den Driessche & Watmough, 2002). For the model (2)–(9) $$R_{0}$$ is given as follows: \begin{align} R_{0}&=max(R_{0s},R_{0m},R_{0x}),\ \ \ \ \ \ where: \\ R_{0s}&=\frac{\beta_{s}^{\alpha}(\varepsilon_{s}^{\alpha}+(1-\lambda_{s}^{\alpha})(d^{\alpha}+t_{1s}^{\alpha}))}{(\varepsilon_{s}^{\alpha}+d^{\alpha}+t_{1s}^{\alpha})(t_{2s}^{\alpha}+\delta_{s}^{\alpha}+d^{\alpha})+\gamma_{s}^{\alpha}(t_{1s}^{\alpha}+d^{\alpha})},\notag\\ R_{0m}&=\frac{\beta_{m}^{\alpha}(\varepsilon_{m}^{\alpha}+(1-\lambda_{m}^{\alpha})d^{\alpha})}{(\varepsilon_{m}^{\alpha}+d^{\alpha})(t_{2m}^{\alpha}+\delta_{m}^{\alpha}+d^{\alpha})+d^{\alpha}\gamma_{m}^{\alpha}},\notag\\ R_{0x}&=\frac{\beta_{x}^{\alpha}(\varepsilon_{x}^{\alpha}+(1-\lambda_{x}^{\alpha})d^{\alpha})}{(\varepsilon_{x}^{\alpha}+d^{\alpha})(t_{2x}^{\alpha}+\delta_{x}^{\alpha}+d^{\alpha})+d^{\alpha}\gamma_{x}^{\alpha}}.\notag \end{align} (10) 3.1. Equilibrium points and their asymptotic stability To discuss the local asymptotic stability for evaluate the equilibrium points, let us consider the following (El-Saka, 2013): $$D^{\alpha}_{t}S=D^{\alpha}_{t}L_{s}=D^{\alpha}_{t}L_{m}=D^{\alpha}_{t}L_{x}=D^{\alpha}_{t}I_{s}=D^{\alpha}_{t}I_{m}=D^{\alpha}_{t}I_{x}=D^{\alpha}_{t}R=0$$. Then, from 1, $$g_{i}(\bar{S}, \bar{L_{s}},\bar{L_{m}}, \bar{L_{x}},\bar{I_{s}}, \bar{I_{m}},\bar{I_{x}},\bar{R})=0$$, $$i=1,2,3,...,8,$$ where $$(\bar{S}, \bar{L_{s}},\bar{L_{m}}, \bar{L_{x}},\bar{I_{s}}, \bar{I_{m}},\bar{I_{x}},\bar{R})$$ denotes any equilibrium point. 3.2. Stability of the disease free equilibrium point Let us consider $$I_s(t)=I_m(t)=I_x(t)=0,$$ then $$L_s(t)=L_m(t)=L_x(t)=0$$, $$R(t)=0,$$ and $$S(t)=N.$$ Then the disease free equilibrium point is $$E_{0}=\{N,0,0,0,0,0,0,0)\}$$. Let us consider the coordinate transformation: $$s(t)=S(t)-\bar{S},$$$$l_{s}(t)=L_{s}(t)-\bar{L_{s}},$$$$l_{m}(t)=L_{m}(t)-\bar{L_{m}},$$$$l_{x}(t)=L_{x}(t)-\bar{L_{x}},$$$$i_{s}(t)=I_{s}(t)-\bar{I_{s}},$$$$i_{m}(t)=I_{m}(t)-\bar{I_{m}},$$$$i_{x}(t)=I_{x}(t)-\bar{I_{x}},$$$$r(t)=R(t)-\bar{R}$$. The corresponding characteristic equation for the disease free equilibrium point is given as follows: $$\displaylines{ (J({E_0}) - \lambda I) = \cr \left( {\matrix{ {\lambda - a} & 0 & 0 & 0 & b & c & {{d_1}} & 0 \cr 0 & {\lambda - e} & 0 & 0 & f & 0 & 0 & 0 \cr 0 & g & {\lambda - h} & 0 & p & q & 0 & 0 \cr 0 & 0 & 0 & {\lambda - r} & 0 & {t_{2m}^\alpha {e^{( - \lambda {\tau _1})}} - p_3^\alpha t_{2m}^\alpha } & t & 0 \cr 0 & u & 0 & 0 & {\lambda - v} & 0 & 0 & 0 \cr 0 & 0 & w & 0 & 0 & {\lambda - x + t_{2m}^\alpha {e^{( - \lambda {\tau _1})}}} & 0 & 0 \cr 0 & 0 & 0 & y & 0 & 0 & {\lambda - z + t_{2x}^\alpha {e^{( - \lambda {\tau _2})}}} & 0 \cr 0 & m & 0 & 0 & n & j & {t_{2x}^\alpha {e^{( - \lambda {\tau _2})}}} & {\lambda - a} \cr } } \right), \cr}$$ where $$a=-d^{\alpha},$$$$b=-\beta_{s}^{\alpha},$$$$c=-\beta_{m}^{\alpha},$$$$d_{1}=-\beta_{x}^{\alpha},$$$$e=-(d^{\alpha}+\varepsilon_{s}^{\alpha}+t_{1s}^{\alpha}),$$$$f=\gamma_{s}^{\alpha}+\lambda_{s}^{\alpha}\beta_{s}^{\alpha},$$$$g=(1-p_{1}^{\alpha})t_{1s}^{\alpha},$$$$h=-(d^{\alpha}+\varepsilon_{m}^{\alpha}),$$$$p=(1-p_{2}^{\alpha})t_{2s}^{\alpha},$$$$q=\gamma_{m}^{\alpha}+\lambda_{m}^{\alpha}\beta_{m}^{\alpha},$$$$r=-(d^{\alpha}+\varepsilon_{x}^{\alpha}),$$$$t=\gamma_{x}^{\alpha}+\lambda_{x}^{\alpha}\beta_{x}^{\alpha},$$$$u=\varepsilon_{s}^{\alpha},$$$$v=-(d^{\alpha}+\delta_{s}^{\alpha}+t_{2s}^{\alpha}+\gamma_{s}^{\alpha}),$$$$w=\varepsilon_{m}^{\alpha},$$$$x=-(d^{\alpha}+\delta_{m}^{\alpha)}+\gamma_{m}^{\alpha}),$$$$y=\varepsilon_{x}^{\alpha},$$$$z=-(d^{\alpha}+\delta_{x}^{\alpha}+\gamma_{x}^{\alpha}),$$$$m=p_{1}^{\alpha}t_{1s}^{\alpha},$$$$n=p_{2}^{\alpha}t_{2s}^{\alpha},$$$$j=p_{3}^{\alpha}t_{2m}^{\alpha}$$. The characteristic equation associated with above matrix is given by El-Saka (2013) as follows: $${\it{\Delta}}(\lambda)=|J(E_{0})-\lambda I|=0,$$ then we have: \begin{align} &(a-\lambda)^{2}(\lambda^{2}-(r+z+t_{2x}^{\alpha}e^{(-\lambda\tau_{2})})\lambda-yt+(z-t_{2x}^{\alpha}e^{(-\lambda\tau_{2})})r)(\lambda^2-(h+x-t_{2m}^{\alpha}e^{(-\lambda\tau_{1})})\lambda \notag \\ &+(x+t_{2m}^{\alpha}e^{(-\lambda\tau_{2})})h-wq)(\lambda^{2}-(e+v)\lambda-uf+ve)=0 \label{RR}. \end{align} (11) Lemma 1 If $$R_{0}<1$$, then the disease free equilibrium point $$E_{0}$$ is locally asymptotically stable for $$\tau_{1}=\tau_{2}=0$$. – Proof. When $$\tau_{1}=\tau_{2}=0,$$ the associated transcendental characteristic equation of system (2)–(9) at $$E_{0}$$ become $${\it{\Delta}}(\lambda)=P(\lambda)=0.$$ Then the eigenvalues of the Jacobian matrix are: \begin{eqnarray*} \lambda_{1,2} &=& -d,\\ \lambda_{3,4}&=&\frac{r+(z-t_{2x}^{\alpha})\pm \sqrt{(r^2-2(z-t_{2x}^{\alpha})r+(z-t_{2x}^{\alpha})^2+4yt)}}{2},\\ \lambda_{5,6}&=&\frac{x-t_{2m}^{\alpha}+h\pm \sqrt{(x-t_{2m}^{\alpha})^2-2(x-t_{2m}^{\alpha})h+h^2+4wq)}}{2},\\ \lambda_{7,8}&=&\frac{v+e\pm \sqrt{(v^2-2ve+e^2+4uf)}}{2}, \end{eqnarray*} these roots are negative or have negative real parts and all eigenvalues satisfies Matignon’s conditions (Matignon, 1996), given by $$\left(|arg\lambda_{i}|>\frac{\alpha\pi}{2}\right)\!,$$ then the disease free equilibrium $$E_{0}$$ is locally asymptotically stable. □ Lemma 2 Let $$R_{0} < 1,$$ then the disease free equilibrium $$E_{0}$$ is locally asymptotically stable for $$\tau_{1}>0$$ and $$\tau_{2}>0$$. Proof. Let us consider $$\tau_{1}>0$$ and $$\tau_{2}>0$$, we noted that second and third factor of the characteristic equation (11), which are $$(\lambda^{2}-(r+z+t_{2x}^{\alpha}e^{(-\lambda\tau_{2})})\lambda-yt+(z-t_{2x}^{\alpha}e^{(-\lambda\tau_{2})})r)$$ and $$(\lambda^2-(h+x-t_{2m}^{\alpha}e^{(-\lambda\tau_{1})})\lambda \&+(x+t_{2m}^{\alpha}e^{(-\lambda\tau_{1})})h-wq),$$ have no pure imaginary roots for any value of the delays $$\tau_{1}, \tau_{2}$$ if $$R_{0} < 1$$. Hence all the roots of the characteristic equation have negative real parts and we get that the disease free equilibrium is locally asymptotically stable regardless of the value of the delay and all eigenvalues satisfies Matignon’s conditions (Matignon, 1996), given by $$\left(|arg\lambda_{i}|>\frac{\alpha\pi}{2}\right)$$ then the disease free equilibrium $$E_{0}$$ is locally asymptotically stable. □ 3.3. Stability of the endemic equilibrium In the absence of controls, i.e., $$u_{1}=u_{2}=0,$$ the system (2)–(9) has endemic equilibrium if at least one of the infected variables is not zero. The expression of analytic for this point is cumbersome and not useful for our purposes (Silva et al., 2017). However, we will determine the stability of the endemic equilibrium of the proposed system numerically. Consider the values of parameters from Table 2. Then the basic reproduction number is $$R_{0}>1$$. The endemic equilibrium $$S=1.0781\times10^{3},$$$$L_{s}=0.0134,$$$$L_{m}=2.1814\times10^{4},$$$$L_{x}=1.8583\times10^{4},$$$$I_{s}=0.0143,$$$$I_{m}=1.5252\times10^{3},$$$$I_{x}=1.1477\times10^{3},$$$$R=1.5852\times10^{4}.$$ The matrices $$A_{1}$$ and $$A_{2}$$ associated to the linearized system at the endemic equilibrium are computed as follows: $${\small A_{1}=\left(\!\!\!\!\begin{array}{cccccccc} -0.6373 & 0 & 0 & 0 & -0.2516 & -0.2516 & -0.2516& 0\\ 0 & -2.0494 & 0 & 0 & 0.5881 & 0 & 0 & 0 \\ 0.1779 & 0.2489& -0.0450 & 0 & 0.2400 & 0.3337 & -0.2545 & 0.0445 \\ 0.1339 & 0.0067 & 0.0067 & -0.0272 & 0 & 0.1200 & 0.4986 & 0.0335\\ 0 & 0.0002 & 0 & 0 &-1.4255& 0 & 0 & 0 \\ 0.0002\times10^{3} & 0 & 0 & 0 & 0.0003\times10^{3} & 1.4095\times10^{3} & 0 & 0.0678 \\ 0.1339 & 0.0067 & 0.0067 & 0.0136& 0 & 0 & 0.0814 & 0.0335 \\ 0 & 1.7600& 0 & 0 &0.8353& -0.0447 & 0.0753 & -0.1695\\ \end{array}\!\!\!\right),}$$ $$A_{2}=\left( \begin{array}{cccccccc} 0& 0 & 0 & 0 & 0& 0 & 0& 0 \\ 0 & 0 & 0 & 0 & 0& 0 & 0 & 0 \\ 0 & 0&0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0& 0 & t_{2m}^{\alpha}e^{-\tau_{1}\lambda} & 0 & 0\\ 0 & 0 & 0 & 0 &0 & 0 & 0 & 0 \\ 0 & 0 & 0& 0 & 0 & -t_{2m}^{\alpha}e^{-\tau_{1}\lambda}& 0 & 0 \\ 0 & 0 & 0 & 0& 0 & 0 & -t_{2x}^{\alpha}e^{-\tau_{2}\lambda}& 0 \\ 0 & 0& 0 & 0 & 0 &0 & t_{2x}^{\alpha}e^{-\tau_{2}\lambda} & 0 \\ \end{array} \right).$$ Table 2 All parameters and their references of the system (2)–(9) Parameters Value Reference $$N$$ 60000 Assumed $$d^{\alpha}$$ $$(1/73.45)^{\alpha}$$$$\left(\frac{1}{year}\right)^{\alpha}$$ http://www.indexmundi.com/EGYPT/life_expectancy_at_birth.html $$\beta_{s}^{\alpha}=\beta_{m}^{\alpha}=\beta_{x}^{\alpha}$$ $$14^{\alpha} \left(\frac{1}{year}\right)^{\alpha}$$ World Health Organization (2012) $$\lambda_{s}^{\alpha}=\lambda_{m}^{\alpha}=\lambda_{x}^{\alpha}$$ $0.5^{\alpha}$$\left(\frac{1}{year}\right)^{\alpha} Sweilam & AL-Mekhlafi (2016b)$$\epsilon_{s}^{\alpha}=\epsilon_{m}^{\alpha}=\epsilon_{x}^{\alpha}0.0002^{\alpha}\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$\alpha_{r1,r2}^{\alpha}0.05^{\alpha}\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$\gamma_{s}^{\alpha}$$=$$\gamma_{m}^{\alpha}$$=$$\gamma_{x}^{\alpha}0.00002^{\alpha}\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$t_{1s}^{\alpha}2^{\alpha}\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$t_{2r}^{\alpha}$$:$$r\in(s,m,x)t_{2s}^{\alpha}=2$$;$$t_{2m}^{\alpha}=t_{2x}^{\alpha}=1\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$\sigma_{r}^{\alpha}$$0.25^{\alpha}$$\left(\frac{1}{year}\right)^{\alpha}$ Sweilam & AL-Mekhlafi (2016b) $$P_{r}^{\alpha}$$ $$0.88^{\alpha}$$$$\left(\frac{1}{year}\right)^{\alpha}$$ Sweilam & AL-Mekhlafi (2016b) $$\delta_{r}^{\alpha}$$ $$0$$ Sweilam & AL-Mekhlafi (2016b) $$T$$ 10 years Assumed Parameters Value Reference $$N$$ 60000 Assumed $$d^{\alpha}$$ $$(1/73.45)^{\alpha}$$$$\left(\frac{1}{year}\right)^{\alpha}$$ http://www.indexmundi.com/EGYPT/life_expectancy_at_birth.html $$\beta_{s}^{\alpha}=\beta_{m}^{\alpha}=\beta_{x}^{\alpha}$$ $$14^{\alpha} \left(\frac{1}{year}\right)^{\alpha}$$ World Health Organization (2012) $$\lambda_{s}^{\alpha}=\lambda_{m}^{\alpha}=\lambda_{x}^{\alpha}$$ $0.5^{\alpha}$$\left(\frac{1}{year}\right)^{\alpha} Sweilam & AL-Mekhlafi (2016b)$$\epsilon_{s}^{\alpha}=\epsilon_{m}^{\alpha}=\epsilon_{x}^{\alpha}0.0002^{\alpha}\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$\alpha_{r1,r2}^{\alpha}0.05^{\alpha}\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$\gamma_{s}^{\alpha}$$=$$\gamma_{m}^{\alpha}$$=$$\gamma_{x}^{\alpha}0.00002^{\alpha}\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$t_{1s}^{\alpha}2^{\alpha}\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$t_{2r}^{\alpha}$$:$$r\in(s,m,x)t_{2s}^{\alpha}=2$$;$$t_{2m}^{\alpha}=t_{2x}^{\alpha}=1\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$\sigma_{r}^{\alpha}$$0.25^{\alpha}$$\left(\frac{1}{year}\right)^{\alpha}$ Sweilam & AL-Mekhlafi (2016b) $$P_{r}^{\alpha}$$ $$0.88^{\alpha}$$$$\left(\frac{1}{year}\right)^{\alpha}$$ Sweilam & AL-Mekhlafi (2016b) $$\delta_{r}^{\alpha}$$ $$0$$ Sweilam & AL-Mekhlafi (2016b) $$T$$ 10 years Assumed Table 2 All parameters and their references of the system (2)–(9) Parameters Value Reference $$N$$ 60000 Assumed $$d^{\alpha}$$ $$(1/73.45)^{\alpha}$$$$\left(\frac{1}{year}\right)^{\alpha}$$ http://www.indexmundi.com/EGYPT/life_expectancy_at_birth.html $$\beta_{s}^{\alpha}=\beta_{m}^{\alpha}=\beta_{x}^{\alpha}$$ $$14^{\alpha} \left(\frac{1}{year}\right)^{\alpha}$$ World Health Organization (2012) $$\lambda_{s}^{\alpha}=\lambda_{m}^{\alpha}=\lambda_{x}^{\alpha}$$ $0.5^{\alpha}$$\left(\frac{1}{year}\right)^{\alpha} Sweilam & AL-Mekhlafi (2016b)$$\epsilon_{s}^{\alpha}=\epsilon_{m}^{\alpha}=\epsilon_{x}^{\alpha}0.0002^{\alpha}\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$\alpha_{r1,r2}^{\alpha}0.05^{\alpha}\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$\gamma_{s}^{\alpha}$$=$$\gamma_{m}^{\alpha}$$=$$\gamma_{x}^{\alpha}0.00002^{\alpha}\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$t_{1s}^{\alpha}2^{\alpha}\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$t_{2r}^{\alpha}$$:$$r\in(s,m,x)t_{2s}^{\alpha}=2$$;$$t_{2m}^{\alpha}=t_{2x}^{\alpha}=1\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$\sigma_{r}^{\alpha}$$0.25^{\alpha}$$\left(\frac{1}{year}\right)^{\alpha}$ Sweilam & AL-Mekhlafi (2016b) $$P_{r}^{\alpha}$$ $$0.88^{\alpha}$$$$\left(\frac{1}{year}\right)^{\alpha}$$ Sweilam & AL-Mekhlafi (2016b) $$\delta_{r}^{\alpha}$$ $$0$$ Sweilam & AL-Mekhlafi (2016b) $$T$$ 10 years Assumed Parameters Value Reference $$N$$ 60000 Assumed $$d^{\alpha}$$ $$(1/73.45)^{\alpha}$$$$\left(\frac{1}{year}\right)^{\alpha}$$ http://www.indexmundi.com/EGYPT/life_expectancy_at_birth.html $$\beta_{s}^{\alpha}=\beta_{m}^{\alpha}=\beta_{x}^{\alpha}$$ $$14^{\alpha} \left(\frac{1}{year}\right)^{\alpha}$$ World Health Organization (2012) $$\lambda_{s}^{\alpha}=\lambda_{m}^{\alpha}=\lambda_{x}^{\alpha}$$ $0.5^{\alpha}$$\left(\frac{1}{year}\right)^{\alpha} Sweilam & AL-Mekhlafi (2016b)$$\epsilon_{s}^{\alpha}=\epsilon_{m}^{\alpha}=\epsilon_{x}^{\alpha}0.0002^{\alpha}\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$\alpha_{r1,r2}^{\alpha}0.05^{\alpha}\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$\gamma_{s}^{\alpha}$$=$$\gamma_{m}^{\alpha}$$=$$\gamma_{x}^{\alpha}0.00002^{\alpha}\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$t_{1s}^{\alpha}2^{\alpha}\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$t_{2r}^{\alpha}$$:$$r\in(s,m,x)t_{2s}^{\alpha}=2$$;$$t_{2m}^{\alpha}=t_{2x}^{\alpha}=1\left(\frac{1}{year}\right)^{\alpha}$$Sweilam & AL-Mekhlafi (2016b)$$\sigma_{r}^{\alpha}$$0.25^{\alpha}$$\left(\frac{1}{year}\right)^{\alpha}$ Sweilam & AL-Mekhlafi (2016b) $$P_{r}^{\alpha}$$ $$0.88^{\alpha}$$$$\left(\frac{1}{year}\right)^{\alpha}$$ Sweilam & AL-Mekhlafi (2016b) $$\delta_{r}^{\alpha}$$ $$0$$ Sweilam & AL-Mekhlafi (2016b) $$T$$ 10 years Assumed Then the transcendental characteristic equation $${\it{\Delta}}\lambda=(\lambda I-A_{1}-A_{2})$$ is given by: \begin{align} &\lambda^{8}-(e^{\tau_{1}\lambda}+e^{\tau_{2}\lambda}+1.4051\times10^3)e^{\tau_{1}\lambda}e^{-\tau_{2}\lambda}\lambda^{7}-(4.4353e^{\tau_{1}\lambda}-1.4052\times10^3e^{\tau_{2}\lambda}-e^{\tau_{1}\lambda}e^{\tau_{2}\lambda}\notag\\&+6.2450\times10^3) e^{-\tau_{1}\lambda}e^{-\tau_{2}\lambda}\lambda^{6}-(6.5237e^{\tau_{1}\lambda}-6.0835\times10^{3}e^{\tau_{2}\lambda}-4.3204e^{\tau_{1}\lambda}e^{\tau_{2}\lambda}+9.1912\times10^{3})\notag\\&e^{-\tau_{1}\lambda}e^{-\tau_{2}\lambda}\lambda^{5}-(3.7527e^{\tau_{1}\lambda}-8.4571\times10^{3}e^{\tau_{2}\lambda}-6.0014e^{\tau_{1}\lambda}e^{\tau_{2}\lambda}+5.2875\times10^{3})e^{-\tau_{1}\lambda}e^{-\tau_{2}\lambda}\lambda^{4}\notag\\ &-(0.8647e^{\tau_{1}\lambda}-4.1905\times10^{3}e^{\tau_{2}\lambda}-2.9696e^{\tau_{1}\lambda}e^{\tau_{2}\lambda}+1.2155\times10^{3})e^{-\tau_{1}\lambda}e^{-\tau_{2}\lambda}\lambda^{3}\notag\\ &+(0.0776e^{\tau_{1}\lambda}-621.5397e^{\tau_{2}\lambda}-0.4361e^{\tau_{1}\lambda}e^{\tau_{2}\lambda}+106.2920)e^{-\tau_{1}\lambda}e^{-\tau_{2}\lambda}\lambda^{2}\notag\\&-(0.0019e^{\tau_{1}\lambda}-27.9671e^{\tau_{2}\lambda}-0.0179e^{\tau_{1}\lambda}e^{\tau_{2}\lambda}+2.3034)e^{-\tau_{1}\lambda}e^{-\tau_{2}\lambda}\lambda\notag\\&+(0.3511e^{\tau_{2}\lambda}+7.1078\times10^{-6}e^{\tau_{1}\lambda}+1.7185\times10^{-4}e^{\tau_{1}\lambda}e^{\tau_{1}\lambda}+0.0233)=0.\label{2a2} \end{align} (12) Consider now the case $$\tau_{1} > 0$$ and $$\tau_{2} > 0,$$ we noted that, the roots of the characteristic equation (12) have no pure imaginary roots for any value of the delay $$\tau_{1}, \tau_{2}$$ if $$R_{0} > 1$$. Hence all the roots of the characteristic equation have negative real parts and all eigenvalues satisfies Matignon’s conditions (Matignon, 1996). Therefore, the endemic equilibrium is locally asymptotically stable. 4. Formulation of the optimal control problem Let us consider the state system (2)–(9), in $$R^{8}$$ with the set of admissible control functions: $$U=\{(u_{1}(\cdot),u_{2}(\cdot))\in (L^{\infty}(0,T)^{2}) \mid 0\leq u_{1}(\cdot),u_{2}(\cdot)\leq 1,\forall t\in[0,T]\}.$$ The objective functional is given by : \begin{align} J(u_{k})=\int_{0}^{T}\eta(S,L_{s},L_{m},L_{x},I_{s},I_{m},I_{x}, R,u_{k},t)dt,\notag \end{align} where $$k=1,2,$$ and \begin{align} J(u_{1}(\cdot),u_{2}(\cdot))&=\int_{0}^{T}[I_{m}(t)+I_{x}(t)+\frac{1}{2}B_{1}u^{2}_{1}(t)+\frac{1}{2}B_{2}u^{2}_{2}(t)] dt\label{22}, \end{align} (13) \begin{align} J(u^{\ast}_{1}(\cdot),u^{\ast}_{2}(\cdot))&= \min_{{\it{\Omega}}} J(u_{1}(\cdot),u_{2}(\cdot)), \end{align} (14) where the constants $$B_{1}$$, $$B_{2}$$, are the measure of the relative cost of the interventions associated with the controls $$u_{1}$$, $$u_{2}$$, respectively. The main point in FDOCP is to find the optimal controls $$u_{k}(t),$$ where $$k=1,2,$$ which minimizes the objective function (13), subject to the following dynamic system: \begin{align} ^{c}_{a}D^{\alpha}_{t}S&=\xi_{1}(S,L_{s},L_{m},L_{x},I_{s},I_{m\tau},I_{m},I_{x},I_{x\tau}, R,u_{k}, t),\notag\\ ^{c}_{a}D^{\alpha}_{t}L_{s}&=\xi_{2}(S,L_{s},L_{m},L_{x},I_{s},I_{m\tau},I_{m},I_{x},I_{x\tau}, R,u_{k},t),\notag\\ ^{c}_{a}D^{\alpha}_{t}L_{m}&=\xi_{3}(S,L_{s},L_{m},L_{x},I_{s},I_{m\tau},I_{m},I_{x},I_{x\tau}, R,u_{k},t),\notag\\ ^{c}_{a}D^{\alpha}_{t}L_{x}&=\xi_{4}(S,L_{s},L_{m},L_{x},I_{s},I_{m\tau},I_{m},I_{x},I_{x\tau}, R,u_{k},t),\notag\\ ^{c}_{a}D^{\alpha}_{t}I_{s}&=\xi_{5}(S,L_{s},L_{m},L_{x},I_{s},I_{m\tau},I_{m},I_{x},I_{x\tau}, R,u_{k},t),\notag\\ ^{c}_{a}D^{\alpha}_{t}I_{m}&=\xi_{6}(S,L_{s},L_{m},L_{x},I_{s},I_{m\tau},I_{m},I_{x},I_{x\tau}, R,u_{k},t),\notag\\ ^{c}_{a}D^{\alpha}_{t}I_{x}&=\xi_{7}(S,L_{s},L_{m},L_{x},I_{s},I_{m\tau},I_{m},I_{x},I_{x\tau}, R,u_{k},t),\notag\\ ^{c}_{a}D^{\alpha}_{t}R&=\xi_{8}(S,L_{s},L_{m},L_{x},I_{s},I_{m\tau},I_{m},I_{x},I_{x\tau}, R,u_{k},t),\notag \end{align} and satisfying the inital conditions: $$\displaylines{ S(\vartheta ) = {\theta _1}(\vartheta ),\,\,{L_s}(\vartheta ) = {\theta _2}(\vartheta ),\,\,{L_m}(\vartheta ) = {\theta _3}(\vartheta ),\,\,{L_x}(\vartheta ) = {\theta _4}(\vartheta ),\,\,{I_s}(\vartheta ) = {\theta _5}(\vartheta ),\,\,{I_m}(\vartheta ) = {\theta _6}(\vartheta ), \cr {I_x}(\vartheta ) = = {\theta _7}(\vartheta ),\quad R(\vartheta ) = {\theta _8}(\vartheta ), \cr}$$ (15)$$\vartheta\in[-\tau,0],$$ where, $$\theta=(\theta_{1},\theta_{2}, \cdot\cdot\cdot \theta_{8})^{'}.$$ The following expression defines a modified objective function: \begin{align} \tilde{J}&=\int_{0}^{T}[H(S(t),L_{s}(t),L_{m}(t),L_{x}(t),I_{s}(t),I_{m}(t), I_{m}(t-\tau),I_{x}(t),I_{x}(t-\tau), R(t),u_{k}(t))\notag\\ &\quad -\sum^{8}_{i=1}\lambda_{i} \xi_{i}(S(t),L_{s}(t),L_{m}(t),L_{x}(t),I_{s}(t),I_{m}(t),I_{m}(t-\tau),I_{x}(t),I_{x}(t-\tau), R(t),u_{k}(t))]dt\label{m2}, \end{align} (16) where $$H(S,L_{s},L_{m},L_{x},I_{s},I_{m},I_{x}, R,u_{k},t)$$ is the following Hamiltonian: \begin{eqnarray} &&H(S,L_{s}(t),L_{m}(t),L_{x}(t),I_{s}(t),I_{m}(t), I_{m}(t-\tau),I_{x}(t), I_{x}(t-\tau), R(t),u_{k}(t),\lambda_{i})\nonumber\\ &&\quad =\eta(S(t),L_{s}(t),L_{m}(t),L_{x}(t),I_{s}(t),I_{m}(t),I_{m}(t-\tau),I_{x}(t), I_{x}(t-\tau) R(t),u_{k}(t))\nonumber\\ &&\qquad +\sum^{8}_{i=1}\lambda_{i}(t) \xi_{i}(S(t),L_{s}(t),L_{m}(t),L_{x}(t),I_{s}(t),I_{m}(t),I_{m}(t-\tau),I_{x}(t),I_{x}(t-\tau), R(t),u_{k}(t)). \label{m3} \end{eqnarray} (17) From (16) and (17), we can derive the necessary and sufficient conditions for FDOPC (Agarwal, 2008) as follows: \begin{eqnarray} &&^{c}_{t}D^{\alpha}_{t_{f}}\lambda_{1}=\frac{\partial H}{\partial S}, \ \ ^{c}_{t}D^{\alpha}_{t_{f}}\lambda_{2}=\frac{\partial H}{\partial L_{s}},\nonumber\\ &&^{c}_{t}D^{\alpha}_{t_{f}}\lambda_{3}=\frac{\partial H}{\partial L_{m}},\ \ ^{c}_{t}D^{\alpha}_{t_{f}}\lambda_{4} =\frac{\partial H}{\partial L_{x}},\label{mh}\\ &&^{c}_{t}D^{\alpha}_{t_{f}}\lambda_{5}= \frac{\partial H}{\partial I_{s}},\ \ ^{c}_{t}D^{\alpha}_{t_{f}}\lambda_{6}=\frac{\partial H}{\partial I_{m}} +X_{[0,t_{f}-\tau]}\frac{\partial H}{\partial I_{m}(t-\tau)}(t+\tau),\nonumber\\ &&^{c}_{t}D^{\alpha}_{t_{f}}\lambda_{7}=\frac{\partial H}{\partial I_{x}(t)}+X_{[0,t_{f}-\tau]}\frac{\partial H}{\partial I_{x}(t-\tau)} (t+\tau),\ \ ^{c}_{t}D^{\alpha}_{t_{f}}\lambda_{8}=\frac{\partial H}{\partial R},\nonumber \end{eqnarray} (18) ${{X}_{[0,{{t}_{f}}-\tau ]}}=\left\{ \begin{array}{*{35}{l}} 1, & t\in \left[ 0,{{t}_{f}}-\tau \right], \\ 0, & otherwise. \\ \end{array} \right.$ (19) \begin{gather} 0=\frac{\partial H}{\partial u_{k}} \label{m4},\\ ^{c}_{0}D^{\alpha}_{t}S=\frac{\partial H}{\partial \lambda_{1}},\ \ ^{c}_{0}D^{\alpha}_{t}L_{s}=\frac{\partial H}{\partial \lambda_{2}},\nonumber\\ ^{c}_{0}D^{\alpha}_{t}L_{m}=\frac{\partial H}{\partial \lambda_{3}},\ \ ^{c}_{0}D^{\alpha}_{t}L_{x}=\frac{\partial H}{\partial \lambda_{4}},\\ ^{c}_{0}D^{\alpha}_{t}I_{s}=\frac{\partial H}{\partial \lambda_{5}}, \ \ ^{c}_{0}D^{\alpha}_{t}I_{m}=\frac{\partial H}{\partial \lambda_{6}},\nonumber\\ ^{c}_{0}D^{\alpha}_{t}I_{x}=\frac{\partial H}{\partial \lambda_{7}}, \ \ ^{c}_{0}D^{\alpha}_{t}R=\frac{\partial H}{\partial \lambda_{8}},\nonumber \end{gather} (20) and it is also required that: $$\lambda_{i}(T)=0, where \ \ i=1,2,3,...,8 \label{m5}.$$ (21) Equations (19) and (21) describe the necessary conditions in terms of a Hamiltonian for the optimal control problem defined above. Theorem 4.1 There exists an optimal pair $$u^{\ast}_{1}(\cdot),$$$$u^{\ast}_{2}(\cdot),$$ corresponding solutions $$S^{\ast}(\cdot),$$$$L^{\ast}_{s}(\cdot),$$$$L^{\ast}_{m}(\cdot),$$$$L^{\ast}_{x}(\cdot),$$$$I^{\ast}_{s}(\cdot),$$$$I^{\ast}_{m}(\cdot),$$$$I^{\ast}_{x}(\cdot),$$ and $$R^{\ast}(\cdot)$$ and that minimizes $$J(u_{1}(\cdot),u_{2}(\cdot))$$ over $${\it{\Omega}}$$. The explicit optimal controls are connected to the existence of continuous specific functions $$\lambda_{i}$$ for $$i=1,2,3,...,8$$ satisfying the follows: (i) adjoint equations: \begin{align} ^{c}_{t}D^{\alpha}_{t_{f}}\lambda_{1}^{\ast}=&-(\lambda_{1}^{\ast}(t)(\frac{\beta_{s}^{\alpha}}{N}I_{s}^{\ast}(t) +\frac{\beta_{m}^{\alpha}}{N}I_{m}^{\ast}(t)+\frac{\beta_{x}^{\alpha}}{N}I_{x}^{\ast}(t) +d^{\alpha})-\lambda_{2}^{\ast}(t)\frac{\lambda_{s}^{\alpha}\beta_{s}^{\alpha}}{N}I_{s}^{\ast}(t) -\lambda_{3}^{\ast}(t)\frac{\lambda_{m}^{\alpha}\beta_{m}^{\alpha}}{N}I_{m}^{\ast}(t) \notag\\&-\lambda_{4}^{\ast}(t)\frac{\lambda_{x}^{\alpha}\beta_{x}^{\alpha}}{N}I_{s}^{\ast}(t) -\lambda_{5}^{\ast}(t)\frac{(1-\lambda_{s}^{\alpha})\beta_{s}^{\alpha}}{N}I_{s}^{\ast}(t) -\lambda_{6}^{\ast}(t)\frac{(1-\lambda_{m}^{\alpha})\beta_{m}^{\alpha}}{N}I_{m}^{\ast}(t)\notag\\ &-\lambda_{7}^{\ast}(t)\frac{(1-\lambda_{x}^{\alpha})\beta_{x}}{N}I_{x}^{\ast}(t))\label{39},\\ \end{align} (22) \begin{align} ^{c}_{t}D^{\alpha}_{t_{f}}\lambda_{2}^{\ast}=&-(-1+\lambda_{2}^{\ast}(t)(\frac{\beta_{s}^{\alpha}\alpha_{ss}^{\alpha}I_{s}^{\ast}(t)}{N} +\frac{\beta_{m}^{\alpha}\alpha_{sm}^{\alpha}I_{m}^{\ast}(t)}{N} +\frac{\beta_{x}^{\alpha}\alpha_{sx}^{\alpha}I_{x}^{\ast}(t)}{N}) +\lambda_{2}^{\ast}(t)(d^{\alpha}+\epsilon_{s}^{\alpha}+t_{1s}^{\alpha}\notag\\ &+\epsilon_{1}u_{1}^{\ast}(t))-\lambda_{3}^{\ast}(t) \frac{\lambda_{m}^{\alpha}\beta_{m}^{\alpha}\alpha_{sm}^{\alpha}}{N}I_{m}^{\ast}(t) -(t_{1s}^{\alpha}+\epsilon_{1}u_{1}^{\ast}(t)-P_{1}^{\alpha}t_{1s}^{\alpha})\lambda_{3}^{\ast}(t)\notag\\ & -\lambda_{4}^{\ast}(t)\frac{\lambda_{x}^{\alpha}\beta_{x}^{\alpha}\alpha_{sx}^{\alpha}}{N}I_{s}^{\ast}(t) -\lambda_{5}^{\ast}(t)\frac{\beta_{s}^{\alpha}\alpha_{ss}^{\alpha}}{N}I_{s}^{\ast}(t)-\lambda_{5}^{\ast}(t)\epsilon_{s}^{\alpha}\notag\\ &-\lambda_{6}^{\ast}(1-\lambda_{m}^{\alpha})\frac{\beta_{m}^{\alpha}\alpha_{sm}^{\alpha}}{N}I_{m}^{\ast}(t)-\lambda_{7}^{\ast}(1-\lambda_{x}^{\alpha})\frac{\beta_{x}^{\alpha}\alpha_{sx}^{\alpha}}{N}I_{x}^{\ast}(t)-P_{1}^{\alpha}t_{1s}^{\alpha}\lambda_{8}^{\ast}(t)),\\ \end{align} (23) \begin{align} ^{c}_{t}D^{\alpha}_{t_{f}}\lambda_{3}^{\ast}=&-\left(\lambda_{3}^{\ast}(t) (\alpha_{mm}^{\alpha}\beta_{m}^{\alpha}\frac{I_{m}^{\ast}(t)}{N}+\alpha_{mx}^{\alpha}\beta_{x}^{\alpha} \frac{I_{x}^{\ast}(t)}{N}+d^{\alpha}+\epsilon_{m}^{\alpha})-\lambda_{4}^{\ast}(t)\alpha_{mx}^{\alpha} \beta_{x}^{\alpha}\lambda_{x}^{\alpha}\frac{I_{x}^{\ast}(t)}{N}\right.\nonumber\\ &\left.-\lambda_{6}^{\ast}(t)\left(\alpha_{mm}^{\alpha}\beta_{m}^{\alpha}\frac{I_{m}^{\ast}(t)}{N}+\epsilon_{m}^{\alpha}\right)-\lambda_{7}^{\ast}(t) (1-\lambda_{x}^{\alpha})\alpha_{mx}^{\alpha}\beta_{^{\alpha}}\frac{I_{x}^{\ast}(t)}{N}\right),\\ \end{align} (24) \begin{align} ^{c}_{t}D^{\alpha}_{t_{f}}\lambda_{4}^{\ast}=&-(\lambda_{4}^{\ast}(t)(\alpha_{xx}^{\alpha}\beta_{x}^{\alpha}\frac{I_{x}^{\ast}(t)}{N} +d^{\alpha}+\epsilon_{x}^{\alpha})-\lambda_{7}^{\ast}(t)(\alpha_{xx}^{\alpha}\beta_{x}^{\alpha}\frac{I_{x}^{\ast}(t)}{N} -\epsilon_{x}^{\alpha}),\\ \end{align} (25) \begin{align} ^{c}_{t}D^{\alpha}_{t_{f}}\lambda_{5}^{\ast}=&-(-1+\lambda_{1}^{\ast}(t)\beta_{s}^{\alpha}\frac{S^{\ast}(t)}{N}+\lambda_{2}^{\ast}(t)(\alpha_{ss}^{\alpha}\beta_{s}^{\alpha}\frac{L_{s}^{\ast}(t)}{N}-\lambda_{s}^{\alpha}\beta_{s}^{\alpha}\frac{S^{\ast}(t)}{N}-\lambda_{s}^{\alpha}\sigma_{s}^{\alpha}\beta_{s}^{\alpha}\frac{R^{\ast}(t)}{N}-\gamma_{s}^{\alpha})\notag\\&+\lambda_{5}^{\ast}(t)((d^{\alpha}+\delta_{s}^{\alpha}+t_{2s}^{\alpha}+\gamma_{s}^{\alpha})-(1-\lambda_{s}^{\alpha})\beta_{s}^{\alpha}\frac{S^{\ast}(t)}{N}-\sigma_{s}^{\alpha}(1-\lambda_{s}^{\alpha})\beta_{s}^{\alpha}\frac{R^{\ast}(t)}{N}\notag\\&+\frac{\beta_{s}^{\alpha}\alpha_{ss}^{\alpha}}{N}L_{s}^{\ast}(t))-\lambda_{3}^{\ast}(t)(t_{2s}^{\alpha}-P_{2}^{\alpha}t_{2s}^{\alpha})+\lambda_{8}^{\ast}(t)(\sigma_{s}^{\alpha}\beta_{s}^{\alpha}\frac{R^{\ast}(t)}{N}-P_{2}^{\alpha}t_{2s}^{\alpha})),\\ \end{align} (26) \begin{align} ^{c}_{t}D^{\alpha}_{t_{f}}\lambda_{6}=&- \left(-1+\lambda_{1}^{\ast}(t)\beta_{m}^{\alpha}\frac{S^{\ast}(t)}{N}+\lambda_{2}^{\ast}(t)\left(\alpha_{sm}^{\alpha}\beta_{m}^{\alpha}\frac{L_{s}^{\ast}(t)}{N}-\lambda_{3}^{\ast}(t)(\lambda_{m}^{\alpha}\beta_{m}^{\alpha}\frac{S^{\ast}(t)}{N}+\lambda_{m}^{\alpha}\sigma_{m}^{\alpha}\beta_{m}^{\alpha}\frac{R^{\ast}(t)}{N}\right.\right.\notag\\ &+\frac{\lambda_{m}^{\alpha}\beta_{m}^{\alpha}\alpha_{sm}^{\alpha}}{N}L_{s}^{\ast}(t)-\frac{\beta_{m}^{\alpha}\alpha_{mm}^{\alpha}}{N}L_{m}^{\ast}(t)+\gamma_{m}^{\alpha})-X_{[0,t_{f}-\tau_{1}]}t_{2m}^{\alpha}\lambda_{4}^{\ast}(t+\tau_{1})-\lambda_{4}^{\ast}(t)\epsilon_{1}u_{2}^{\ast}(t)\notag\\ &+\lambda_{4}^{\ast}(t)P_{3}^{\alpha}t_{2m}^{\alpha}-\lambda_{6}^{\ast}(t) (\frac{\beta_{m}^{\alpha}\alpha_{mm}^{\alpha}}{N}L_{m}^{\ast}(t)+ (1-\lambda_{m}^{\alpha})\beta_{m}^{\alpha}\frac{S^{\ast}(t)}{N}+\sigma_{m}^{\alpha}(1-\lambda_{m}^{\alpha})\beta_{m}^{\alpha} \frac{R^{\ast}(t)}{N}\notag\\ &+\alpha_{sm}^{\alpha} (1-\lambda_{m}^{\alpha})\beta_{m}^{\alpha}\frac{L_{s}^{\ast}(t)}{N}-(d^{\alpha} +\delta_{m}^{\alpha}+X_{[0,t_{f}-\tau_{1}]}t_{2m}^{\alpha}\lambda_{6}^{\ast} (t+\tau_{1})+\epsilon_{1}u_{1}^{\ast}(t)+\gamma_{m}^{\alpha}))\notag\\ &\left.\left.-\lambda_{8}^{\ast}(t)\left(P_{3}^{\alpha}t_{2m}^{\alpha}-\sigma_{m}^{\alpha}\beta_{m}^{\alpha}\frac{R^{\ast}}{N}\right)\right)\right), \end{align} (27) \begin{align} ^{c}_{t}D^{\alpha}_{t_{f}}\lambda_{7}^{\ast}=&- (-1+\lambda_{1}^{\ast}(t)\beta_{x}^{\alpha}\frac{S^{\ast}(t)}{N}+\lambda_{2}^{\ast}(t)\alpha_{sx}^{\alpha}\beta_{x}^{\alpha}\frac{L_{s}^{\ast}(t)}{N} +\lambda_{3}^{\ast}(t)\alpha_{mx}^{\alpha}\beta_{x}^{\alpha}\frac{L_{m}^{\ast}(t)}{N}-\lambda_{4}^{\ast}(t) (\lambda_{x}\sigma_{x}^{\alpha}\beta_{x}^{\alpha}\frac{R^{\ast}(t)}{N}\notag\\ &+\lambda_{x}^{\alpha}\beta_{x}^{\alpha}\frac{S^{\ast}(t)}{N}+\lambda_{x}^{\alpha}\alpha_{sx}^{\alpha}\beta_{x}^{\alpha}\frac{L_{s}^{\ast}(t)}{N}+\lambda_{x}^{\alpha}\alpha_{mx}^{\alpha}\beta_{x}^{\alpha}\frac{L_{m}^{\ast}(t)}{N}+\gamma_{x}^{\alpha}-\alpha_{xx}^{\alpha}\beta_{x}^{\alpha}\frac{L_{x}^{\ast}(t)}{N})\notag\\&-\lambda_{7}^{\ast}(t)(\alpha_{xx}^{\alpha}\beta_{x}^{\alpha}\frac{L_{x}^{\ast}(t)}{N}+(1-\lambda_{x}^{\alpha})\beta_{x}^{\alpha}\frac{S^{\ast}}{N}+(1-\lambda_{x}^{\alpha})\beta_{x}^{\alpha}\alpha_{sx}^{\alpha}\frac{L_{s}^{\ast}}{N}+(1-\lambda_{x}^{\alpha})\beta_{x}^{\alpha}\sigma_{x}^{\alpha}\frac{R^{\ast}}{N}\notag\\&+(1-\lambda_{x}^{\alpha})\beta_{x}^{\alpha}\alpha_{mx}^{\alpha}\frac{L_{m}^{\ast}}{N}+(d^{\alpha}+\delta_{x}^{\alpha}+\epsilon_{2}u_{2}^{\ast}+\gamma_{x}^{\alpha}))+X_{[0,t_{f}-\tau_{2}]}t_{2x}^{\alpha}\lambda_{8}^{\ast}(t+\tau_{2}) -X_{[0,t_{f}-\tau_{2}]}\notag\\&\times t_{2x}^{\alpha}\lambda_{7}^{\ast}(t+\tau_{2})-\lambda_{8}^{\ast}(t)(\epsilon_{2}u_{2}^{\ast}-\sigma_{x}^{\alpha}\beta_{x}^{\alpha}\frac{R}{N})),\\ \end{align} (28) \begin{align} ^{c}_{t}D^{\alpha}_{t_{f}}\lambda_{8}^{\ast}=& - \left( -\lambda_{2}^{\ast}(t)\beta_{s}^{\alpha}\lambda_{s}^{\alpha}\sigma_{s}^{\alpha}\frac{I_{s}^{\ast}(t)}{N} -\lambda_{3}^{\ast}(t)\beta_{m}^{\alpha}\lambda_{m}^{\alpha}\sigma_{m}^{\alpha}\frac{I_{m}^{\ast}(t)}{N} -\lambda_{4}^{\ast}(t)\beta_{x}^{\alpha}\lambda_{x}^{\alpha}\sigma_{x}^{\alpha}\frac{I_{x}^{\ast}(t)}{N}\right.\notag\\ &-\lambda_{5}^{\ast}(t)(1-\lambda_{s}^{\alpha})\beta_{s}^{\alpha}\sigma_{s}^{\alpha}\frac{I_{s}^{\ast}(t)}{N} -\lambda_{6}^{\ast}(t)(1-\lambda_{m}^{\alpha})\beta_{m}^{\alpha}\sigma_{m}^{\alpha}\frac{I_{m}^{\ast}(t)}{N}\notag\\ &\left.-\lambda_{7}^{\ast}(t)(1-\lambda_{x}^{\alpha})\beta_{x}^{\alpha}\sigma_{x}^{\alpha}\frac{I_{x}^{\ast}(t)}{N} +\lambda_{8}^{\ast}(t)\left(\sigma_{s}^{\alpha}\beta_{s}^{\alpha}\frac{I_{s}^{\ast}}{N} +\sigma_{m}^{\alpha}\beta_{m}^{\alpha}\frac{I_{m}^{\ast}}{N}+\sigma_{x}^{\alpha}\beta_{x}^{\alpha}\frac{I_{x}^{\ast}}{N}+d^{\alpha}\right)\right) \label{46}, \end{align} (29) (ii) with transversality conditions $$\lambda_{i}^{\ast}(T)=0$$, $$i = 1, . . . , 8$$. ${{X}_{[0,{{t}_{f}}-\tau ]}}=\left\{ \begin{array}{*{35}{l}} 1, & t\in \left[ 0,{{t}_{f}}-\tau \right], \\ 0, & otherwise. \\ \end{array} \right.$ (iii) optimality conditions: \begin{align} &H(S^{\ast}(t),L_{s}^{\ast}(t),L_{m}^{\ast}(t), L_{x}^{\ast}(t),I_{s}^{\ast}(t),I_{m}^{\ast}(t),I_{m}^{\ast}(t-\tau_{1}),I_{x}^{\ast}(t), I_{x}^{\ast}(t-\tau_{2}), R^{\ast}(t),\lambda^{\ast}(t), u_{k}^{\ast}(t) )= \notag\\ &\min _{0\leq u_{k}\leq 1}H(S^{\ast}(t),L_{s}^{\ast}(t),L_{m}^{\ast}(t),L_{x}^{\ast}(t),I_{s}^{\ast}(t),I_{m}^{\ast}(t), I_{m}^{\ast}(t-\tau_{1}),I_{x}^{\ast}(t), I_{x}^{\ast}(t-\tau_{2}) R^{\ast}(t),\lambda^{\ast}(t), u_{k}(t)),\label{7m}& \end{align} (30) \begin{align} u_{1}^{\ast}(t)&=\min\{\max\{0,\frac{\epsilon_{1}I_{m}^{\ast}(\lambda_{6}^{\ast}(t)-\lambda_{4}^{\ast}(t))}{B_{1}}\},1\}\label{10m},\\ \end{align} (31) \begin{align} u_{2}^{\ast}(t)&=\min\{\max\{0,\frac{\epsilon_{2}I_{x}^{\ast}(\lambda_{7}^{\ast}(t)-\lambda_{8}^{\ast}(t))}{B_{2}}\},1\}\label{11m}, \end{align} (32) where the stationarity condition is $$\frac{\partial H}{\partial u_{k}}=0,$$$$k=1,2.$$ Proof. Using conditions (18), we can get (22)–(29), where the Hamiltonian $$H$$ is given by: \begin{align} H&=H(S, L_{s}, L_{m}, L_{x}, I_{s},I_{m}, I_{x}, R, \lambda, u_{k} )\notag\\ &=I_{s}+I_{m}+I_{x}+L_{s}+u_{1}^{2}\frac{B_{1}}{2}+u_{2}^{2}\frac{B_{2}}{2}\notag\\ &\quad +\lambda_{1}(b^{\alpha}-d^{\alpha}S-\beta_{s}^{\alpha}\frac{SI_{s}}{N}-\beta_{m}^{\alpha}\frac{SI_{m}}{N}-\beta_{x}^{\alpha}\frac{SI_{x}}{N})\notag\\ &\quad +\lambda_{2}(\lambda_{s}^{\alpha}\beta_{s}^{\alpha}\frac{SI_{s}}{N}+\sigma_{s}^{\alpha}\lambda_{s}^{\alpha}\beta_{s}^{\alpha}\frac{RI_{s}}{N}-\alpha_{ss}^{\alpha}\beta_{s}^{\alpha}\frac{L_{s}I_{s}}{N}-\alpha_{sm}^{\alpha}\beta_{m}^{\alpha}\frac{L_{s}I_{m}}{N}-\alpha_{sx}^{\alpha}\beta_{x}^{\alpha}\frac{L_{s}I_{x}}{N}+\gamma_{s}^{\alpha}I_{s} \notag\\ &\quad -(d^{\alpha}+\epsilon_{s}^{\alpha}+t_{1s}^{\alpha})L_{s}\notag\\ &\quad +\lambda_{3}(\lambda_{m}^{\alpha}\beta_{m}^{\alpha}\frac{SI_{m}}{N}+\sigma_{m}^{\alpha}\lambda_{m}^{\alpha}\beta_{m}^{\alpha}\frac{RI_{m}}{N}+\alpha_{sm}^{\alpha}\beta_{m}^{\alpha}\lambda_{m}^{\alpha}\frac{L_{s}I_{m}}{N}-\alpha_{mm}^{\alpha}\beta_{m}^{\alpha}\frac{L_{m}I_{m}}{N}-\alpha_{mx}^{\alpha}\beta_{x}^{\alpha}\frac{L_{m}I_{x}}{N}\notag\\ &\quad +\gamma_{m}^{\alpha}I_{m}-(d^{\alpha}+\epsilon_{m}^{\alpha})L_{m}+t_{1s}^{\alpha}L_{s} -P_{1}^{\alpha}t_{1s}^{\alpha}L_{s}+t_{2s}^{\alpha}I_{s}-P_{2}^{\alpha}t_{2s}^{\alpha}I_{s},)\notag\\ &\quad +\lambda_{4}(\lambda_{x}^{\alpha}\beta_{x}^{\alpha}\frac{SI_{x}}{N}+\sigma_{x}^{\alpha}\lambda_{x}^{\alpha}\beta_{x}^{\alpha}\frac{RI_{x}}{N}+\alpha_{sx}^{\alpha}\beta_{x}^{\alpha}\lambda_{x}^{\alpha}\frac{L_{s}I_{x}}{N} +\alpha_{mx}^{\alpha}\beta_{x}^{\alpha}\lambda_{x}^{\alpha}\frac{L_{m}I_{x}}{N}-\alpha_{xx}^{\alpha}\beta_{x}^{\alpha}\frac{L_{x}I_{x}}{N} \notag\\ &\quad -(d^{\alpha}+\epsilon_{x}^{\alpha})L_{x}+\gamma_{x}^{\alpha}I_{x}+t_{2m}^{\alpha}I_{m}(t-\tau_{1})+\epsilon_{1}u_{1}(t)I_{m}-P_{3}^{\alpha}t_{2m}^{\alpha}I_{m})\notag\\ &\quad +\lambda_{5}(\alpha_{ss}^{\alpha}\beta_{s}^{\alpha}\frac{L_{s}I_{s}}{N}+(1-\lambda_{s}^{\alpha})\beta_{s}^{\alpha}(\frac{SI_{s}}{N}+\sigma_{s}^{\alpha}\frac{RI_{s}}{N})+\epsilon_{s}^{\alpha}L_{s}-(d^{\alpha}+\delta_{s}^{\alpha}+t_{2s}^{\alpha}+\gamma_{s}^{\alpha})\notag\\ &\quad +\lambda_{6}(\alpha_{mm}^{\alpha}\beta_{m}^{\alpha}\frac{L_{m}I_{m}}{N}+(1-\lambda_{m}^{\alpha})\beta_{m}^{\alpha}(\frac{SI_{m}}{N}+\sigma_{m}^{\alpha}\frac{RI_{m}}{N}+\alpha_{sm}^{\alpha}\frac{L_{s}I_{m}}{N})+\epsilon_{m}^{\alpha}L_{m} \notag\\ &\quad -(d^{\alpha}+\delta_{m}^{\alpha}+\epsilon_{1}u_{1}(t)+\gamma_{m}^{\alpha})I_{m}-t_{2m}^{\alpha}I_{m}(t-\tau_{1}))\notag\\ &\quad +\lambda_{7}(\alpha_{xx}^{\alpha}\beta_{x}^{\alpha}\frac{L_{x}I_{x}}{N} +(1-\lambda_{x}^{\alpha})\beta_{x}^{\alpha}(\frac{SI_{x}}{N}+\sigma_{x}^{\alpha}\frac{RI_{x}}{N} +\alpha_{sx}^{\alpha}\frac{L_{s}I_{x}}{N}+\alpha_{mx}^{\alpha}\frac{L_{m}I_{x}}{N})+\epsilon_{x}^{\alpha}L_{x} \notag\\ &\quad -(d^{\alpha}+\delta_{x}^{\alpha}+\gamma_{x}^{\alpha}+\epsilon_{2}u_{2}(t))I_{x}-t_{2x}^{\alpha}I_{x}(t-\tau_{2}))\notag\\ &\quad +\lambda_{8}(P_{1}^{\alpha}t_{1s}^{\alpha}L_{s}+P_{2}^{\alpha}t_{2s}^{\alpha}I_{s}+P_{3}^{\alpha}t_{2m}^{\alpha}I_{m}+t_{2x}^{\alpha}I_{x}+t_{2x}^{\alpha}I_{x}(t-\tau_{2})+\epsilon_{2}u_{2}(t)I_{x}-\sigma_{s}^{\alpha}\beta_{s}^{\alpha}\frac{RI_{s}}{N} \notag\\ &\quad -\sigma_{m}^{\alpha}\beta_{m}^{\alpha}\frac{RI_{m}}{N}-\sigma_{x}^{\alpha}\beta_{x}^{\alpha}\frac{RI_{x}}{N} -d^{\alpha}R)\label{6m}, \end{align} (33) and $$\lambda_{i},$$$$i=1,2,3,...,8$$ are the Lagrange multipliers. It is known as a co-state or adjoint variables. Moreover, the transversality conditions $$\lambda_{i}(T ) =0,$$$$i = 1, . . . , 8,$$ hold and the optimal controls (31)–(32) can be claimed from the minimization condition (30). □ 5. Numerical method 5.1. Generalized Euler method In the following, we will use the generalization of the classical Euler’s method, for solving the optimal control delay problem. For more details on this method, see Odibat & Momani (2008). The headlines of this method is given as follows, consider the problem (2)–(9) and (15). Assume that $$\tau_1$$ and $$\tau_2$$ are positive. The general formula for GEM when $t_{k+1}$$=$$t_{k}+h,$ for $$k = 0, 1, ..., Q-1,$$ is given as follows: \begin{align} S(t_{k+1})&=S(t_{k})+\frac{h^{\alpha}}{{\it{\Gamma}}(\alpha+1)}f_{1}(S(t_{k}),L_{s}(t_{k}),L_{m}(t_{k}),L_{x}(t_{k}),I_{s}(t_{k}),I_{m}(t_{k}) ,I_{m}(t_{\tau_{1}k}), I_{x}(t_{k}),\notag\\ &\quad \times I_{x}(t_{\tau_{2} k}),R(t_{k}), u_{1}, u_{2},t_{k}),\notag\\ L_{s}(t_{k+1})&=L_{s}(t_{k})+\frac{h^{\alpha}}{{\it{\Gamma}}(\alpha+1)}f_{2}(S(t_{k}),L_{s}(t_{k}),L_{m}(t_{k}),L_{x}(t_{k}),I_{s}(t_{k}),I_{m}(t_{k}),I_{m}(t_{\tau_{1} k}), I_{x}(t_{k}), \notag\\ &\quad \times I_{x}(t_{\tau_{2} k}),R(t_{k}),u_{1}, u_{2},t_{k}),\notag\\ L_{m}(t_{k+1})&=L_{m}(t_{k})+\frac{h^{\alpha}}{{\it{\Gamma}}(\alpha+1)}f_{3}(S(t_{k}),L_{s}(t_{k}),L_{m}(t_{k}),L_{x}(t_{k}),I_{s}(t_{k}), I_{m}(t_{k}),I_{m}(t_{\tau_{1} k}), I_{x}(t_{k}),\notag\\ &\quad \times I_{x}(t_{\tau_{2} k}),R(t_{k}), u_{1}, u_{2},t_{k}),\notag\\ L_{x}(t_{k+1})&=L_{x}(t_{k})+\frac{h^{\alpha}}{{\it{\Gamma}}(\alpha+1)}f_{4}(S(t_{k}),L_{s}(t_{k}),L_{m}(t_{k}),L_{x}(t_{k}),I_{s}(t_{k}) , I_{m}(t_{k}), I_{m}(t_{\tau_{1} k}), I_{x}(t_{k}),\notag\\ &\quad \times I_{x}(t_{\tau_{2} k}),R(t_{k}),u_{1}, u_{2},t_{k}),\notag\\ I_{s}(t_{k+1})&=I_{s}(t_{k})+\frac{h^{\alpha}}{{\it{\Gamma}}(\alpha+1)}f_{5}(S(t_{k}),L_{s}(t_{k}),L_{m}(t_{k}),L_{x}(t_{k}),I_{s}(t_{k}), I_{m}(t_{k}), I_{m}(t_{\tau_{1} k}), I_{x}(t_{k}),\notag\\ &\quad \times I_{x}(t_{\tau_{2} k}),R(t_{k}),u_{1}, u_{2},t_{k}),\notag\\ I_{m}(t_{k+1})&=I_{m}(t_{k})+\frac{h^{\alpha}}{{\it{\Gamma}}(\alpha+1)}f_{6}(S(t_{k}),L_{s}(t_{k}),L_{m}(t_{k}),L_{x}(t_{k}),I_{s}(t_{k}) , I_{m}(t_{k}),I_{m}(t_{\tau_{1} k}), I_{x}(t_{k}), \notag\\ &\quad \times I_{x}(t_{\tau_{2} k}),R(t_{k}),u_{1}, u_{2},t_{k}),\notag\\ I_{x}(t_{k+1})&=I_{x}(t_{k})+\frac{h^{\alpha}}{{\it{\Gamma}}(\alpha+1)}f_{7}(S(t_{k}),L_{s}(t_{k}),L_{m}(t_{k}),L_{x}(t_{k}),I_{s}(t_{k}) , I_{m}(t_{k}), I_{m}(t_{\tau_{1} k}), I_{x}(t_{k}), \notag\\& \quad \times I_{x}(t_{\tau_{2} k}),R(t_{k}),u_{1}, u_{2},t_{k}),\notag\\ R(t_{k+1})&=R(t_{k})+\frac{h^{\alpha}}{{\it{\Gamma}}(\alpha+1)}f_{8}(S(t_{k}),L_{s}(t_{k}),L_{m}(t_{k}),L_{x}(t_{k}),I_{s}(t_{j}), I_{m}(t_{k}),I_{m}(t_{\tau_{1} k}), I_{x}(t_{k}), \notag\\ &\quad \times I_{x}(t_{\tau_{2} k}),R(t_{k}),u_{1}, u_{2},t_{k}).\notag\ \end{align} It is clear that if $$\alpha = 1,$$ then the GEM reduces to the classical Euler’s method (Baleanu et al., 2012). To solve (22)–(29) with transversality conditions $$\lambda_{i}^{\ast}(T)=0$$, $$i = 1, . . . , 8$$ and $$\tau_{1}>0$$, $$\tau_{2}>0$$ by using backward Euler method. 5.2. Iterative optimal control method This method is a simple technique to study numerically fractional delay control problems, for more details see Sweilam & AL-Mekhlafi (2017), and the references cited therein. The main advantage of this method is to avoid the expected difficulties during computations such as integration by parts rule in case of fractional calculus. IOCM could adapt the same existence and characteristic optimal control theorems as in DDE systems for the system of fractional DDEs. This method can be described as follows: Algorithm 1: Step 0. Choose a starting point $$(S_{0}, L_{s0}, L_{m0},L_{x0},I_{s0},I_{m0},I_{x0},R_{0},)$$ and choose a tolerance $$\epsilon>0,$$ for the accuracy desired. Then set $$j = 1$$. Step 1. Let $$\tau_{1}$$, $$\tau_{2},$$ positive numberes, set $$u_{1}(t)$$ and $$u_{2}(t)$$ as given in (31) and (32) respectively. Solve the co-state system (22)–(29) with transversality conditions $$\lambda_{i}(T)=0,$$$$j = 1,2, 3,...,8,$$ then go to Step 3. Step 2. Plugging the values of $$u_{1}(t)$$ and $$u_{2}(t)$$ into the system (2)–(9), then solve it with the same starting point to obtain the new starting point $$(S(t)$$, $$L_{s}(t),$$$$L_{m}(t)$$,$$L_{x}(t)$$,$$I_{s}(t)$$, $$I_{m}(t)$$,$$I_{x}(t)$$, $$R(t))$$. Step 3. The stopping criterion is as follows: if $$|u_{j} -u_{j + 1}|< \epsilon$$ stop, otherwise put $$j = j + 1$$, and go to Step 1. 5.3. Discretizations of IOCM To solve optimization problem (2)–(9) and (13), first we need to discretized the state system (2)–(9). The Grünwald–Letnikov’s definition is used to discretization the fractional derivative $$^{c}_{a}D^{\alpha}_{t}f(t)$$. Then the system (2)–(9), can be written in the following form: \begin{align} S^{n}=&\frac{b^{\alpha}-\sum^{n}_{j=1}q^{\alpha}_{j}S^{n-j}}{q^{-\alpha}_{0}+d^{\alpha}+\frac{\beta_{s}^{\alpha}I_{s}^{n}+\beta_{m}^{\alpha}I_{m}^{n}+\beta_{x}^{\alpha}I_{x}^{n}}{N^{n}}},\\ \end{align} (34) \begin{align} L_{s}^{n}=&\frac{\frac{\beta_{s}^{\alpha}I_{s}^{n}}{N^{n}}\lambda_{s}^{\alpha}(S^{n}+\sigma_{s}^{\alpha}R^{n})+\gamma_{s}^{\alpha}I_{s}^{n}-\sum^{n}_{j=1}q^{\alpha}_{j}L_{s}^{n-j}}{q^{-\alpha}_{0}+(d^{\alpha}+t_{1s}^{\alpha}+\epsilon_{s}^{\alpha})+\frac{1}{N^{n}}(\alpha_{ss}^{\alpha}\beta_{s}^{\alpha}I_{s}^{n}+\alpha_{sm}^{\alpha}\beta_{m}^{\alpha}I_{m}^{n}+\alpha_{sx}^{\alpha}\beta_{x}^{\alpha}I_{x}^n)},\\ \end{align} (35) \begin{align} L_{m}^{n}=&\frac{\frac{\beta_{m}^{\alpha}\lambda_{m}^{\alpha}I_{m}^{n}}{N^{n}}(S^{n}+\sigma_{m}^{\alpha}R^{n}+\alpha_{sm}^{\alpha}L_{s}^{n})+\gamma_{m}^{\alpha}I_{m}^{n}+t_{1s}^{\alpha}L_{s}^{n}(1-P_{1}^{\alpha})}{q^{-\alpha}_{0}+(d^{\alpha}+\epsilon_{m}^{\alpha})+\frac{1}{N^{n}}(\alpha_{mm}^{\alpha}\beta_{m}^{\alpha}I_{m}^{n}+\alpha_{mx}^{\alpha}\beta_{x}^{\alpha}I_{x}^{n})} \notag\\ &+\frac{t_{2s}^{\alpha}I_{s}^{n}(1-P_{2}^{\alpha})-\sum^{n}_{j=1}q^{\alpha}_{j}L_{m}^{n-j}}{\omega^{-\alpha}_{0}+(d^{\alpha}+\epsilon_{m}^{\alpha})+\frac{1}{N^{n}}(\alpha_{mm}^{\alpha}\beta_{m}^{\alpha}I_{m}^{n}+\alpha_{mx}^{\alpha}\beta_{x}^{\alpha}I_{x}^{n})},\\ \end{align} (36) \begin{align} L_{x}^{n}=&\frac{\frac{\beta_{x}^{\alpha}\lambda_{x}^{\alpha}I_{x}^{n}}{N^{n}}(S^{n}+\sigma_{x}^{\alpha}R^{n}+\alpha_{sx}^{\alpha}L_{s}^{n}+\alpha_{mx}^{\alpha}L_{m}^{n})+t_{2m}^{\alpha}I_{m}^{n\tau_{1}}+\epsilon_{1}u_{1}^{n}I_{m}^{n}-P_{3}^{\alpha}t_{2m}^{\alpha}I_{m}^{n}}{q^{-\alpha}_{0}+(d^{\alpha}+\epsilon_{x}^{\alpha})+\frac{1}{N^{n}}(\alpha_{xx}^{\alpha}\beta_{x}^{\alpha}I_{x}^{n})} \notag\\ &+\frac{\gamma_{x}^{\alpha}I_{x}^{n}-\sum^{n}_{j=1}q^{\alpha}_{j}L_{x}^{n-j}}{q^{-\alpha}_{0}+(d^{\alpha}+\epsilon_{x}^{\alpha})+\frac{1}{N^{n}}(\alpha_{xx}^{\alpha}\beta_{x}^{\alpha}I_{x}^{n})},\\ \end{align} (37) \begin{align} I_{s}^{n}=&\frac{\varphi_{5}(h)\beta_{s}^{\alpha}\frac{I_{s}^{n}}{N^{n}}(\alpha_{ss}^{\alpha}L_{s}^{n}+(1-\lambda_{s}^{\alpha})(S^{n}+\sigma_{s}^{\alpha}R^{n}))}{q^{-\alpha}_{0}+(d^{\alpha}+\delta_{s}^{\alpha}+\gamma_{s}^{\alpha}+t_{2s}^{\alpha}} \notag\\ &+\frac{\epsilon_{s}^{\alpha}L_{s}^{n}-\sum^{n}_{j=1}q^{\alpha}_{j}I_{s}^{n-j}}{q^{-\alpha}_{0}+(d^{\alpha}+\delta_{s}^{\alpha}+\gamma_{s}^{\alpha}+t_{2s}^{\alpha})},\\ \end{align} (38) \begin{align} I_{m}^{n}=&\frac{\beta_{m}^{\alpha}\frac{I_{m}^{n}}{N^{n}}(\alpha_{mm}^{\alpha}L_{m}^{n}+(1-\lambda_{m}^{\alpha})(S^{n}+\sigma_{m}^{\alpha}R^{n}+\alpha_{sm}^{\alpha}L_{s}^{n}))}{q^{-\alpha}_{0}+(d^{\alpha}+\delta_{m}^{\alpha}+\gamma_{m}^{\alpha}+\epsilon_{1}u_{1}^{n}))+t_{2m}^{\alpha}I_{m}^{n\tau_{1}}} \notag\\ &+\frac{\epsilon_{m}^{\alpha}L_{m}^{n}-\sum^{n}_{j=1}q^{\alpha}_{j}I_{m}^{n-j}}{q^{-\alpha}_{0}+(d^{\alpha}+\delta_{m}^{\alpha}+\gamma_{m}^{\alpha}+\epsilon_{1}u_{1}^{n})+t_{2m}^{\alpha}I_{m}^{n\tau_{1}}},\\ \end{align} (39) \begin{align} I_{x}^{n}=&\frac{\beta_{x}^{\alpha}\frac{I_{x}^{n}}{N^{n}}(\alpha_{xx}^{\alpha}L_{x}^{n}+(1-\lambda_{x}^{\alpha})(S^{n}+\sigma_{x}^{\alpha}R^{n}+\alpha_{sx}^{\alpha}L_{s}^{n}+\alpha_{mx}^{\alpha}L_{m}^{n}))}{q^{-\alpha}_{0}+(d^{\alpha}+\delta_{x}^{\alpha}+\gamma_{x}^{\alpha}+\epsilon_{2}u_{2}^{n})+t_{2x}^{\alpha}I_{x}^{n\tau_{1}}} \notag\\ &+\frac{\epsilon_{x}L_{x}^{n}-\sum^{n}_{j=1}q^{\alpha}_{j}I_{x}^{n-j}}{q^{-\alpha}_{0}+(d^{\alpha}+\delta_{x}^{\alpha}+\gamma_{x}^{\alpha}+\epsilon_{2}u_{2}^{n})+t_{2x}^{\alpha}I_{x}^{n\tau_{2}}},\\ \end{align} (40) \begin{align} R^{n}=&\frac{t_{1s}^{\alpha}P_{1}^{\alpha}L_{s}^{n}+P_{2}^{\alpha}t_{2s}^{\alpha}I_{s}^{n}+t_{2m}^{\alpha}P_{3}^{\alpha}I_{m}^{n}+t_{2x}^{\alpha}I_{x}^{n\tau_{2}}+\epsilon_{2}u_{2}^{n}I_{x}^{n}-\sum^{n}_{j=1}q^{\alpha}_{j}R^{n-j}}{q^{-\alpha}_{0}+d^{\alpha}+\frac{1}{N^{n}}(\sigma_{s}^{\alpha}\beta_{s}^{\alpha}I_{s}^{n}+\sigma_{m}^{\alpha}\beta_{m}^{\alpha}I_{m}^{n}+\sigma_{x}^{\alpha}\beta_{x}^{\alpha}I_{x}^{n})}. \end{align} (41) To solve (22)–(29) with transversality conditions $$\lambda_{i}^{\ast}(T)=0$$, $$i = 1, . . . , 8,$$ we will use the backward Euler method. 6. Numerical results and discussions In the following, IOCM and GEM are used to obtain the approximate solutions for systems (2)–(9) and (22)–(29). Using the initial conditions ($$S(0)$$, $$L_{s}(0)$$, $$L_{m}(0)$$, $$L_{x}(0)$$, $$I_{s}(0),$$$$I_{m}(0),$$$$I_{x}(0),$$$$R(0))$$$$=$$ ($$\frac{76}{120}N,$$$$\frac{20}{120}N,$$$$\frac{5}{120}N,$$$$\frac{2}{120}N,$$$$\frac{8}{120}N,$$$$\frac{4}{120}N,$$$$\frac{2}{120}N,$$$$\frac{3}{120}N$$), and the parameters given in Table 2. The approximate solutions of the proposed system are given in Figs 1–9 at different values of $$\tau_{1},$$$$\tau_{2}$$ and $$0<\alpha\leq1$$. In Fig. 1, the effect of different time delays $$\tau_{1},$$ with fixed time delay $$\tau_{2},$$ on the behavior the approximate solutions for state variables $$I_{m}(t-\tau_{1}),$$$$I_{x}(t-\tau_{2})$$ in the controlled case. The delay $$\tau_{1}$$ was chosen to equal $$1, 2$$, and $$3$$. The delay $$\tau _2$$ is chosen in this case fixed and equal $$2$$. Moreover, we comper the obtianed results with the non-delay case, i.e., $$\tau_{1}=\tau_{2}=0$$. We noted that, the result of $$I_{x}(t-\tau_{2})$$ is not affected by changing the value of $$\tau_{1}$$. Also, Fig. 2, shows the effect of different time delays $$\tau_{2},$$ with fixed $$\tau_{1},$$ on the behavior of the approximate solutions for state variables $$I_{m}(t-\tau_{1}),$$$$I_{x}(t-\tau_{2})$$ in the controlled case. The delay $$\tau_{2}$$ was chosen to equal $$1, 2, 3,$$ and the delay $$\tau_{1}=1.$$ Moreover, the results are compared with non-delay case. We noted that, the results of $$I_{m}(t-\tau_{1})$$ is not affected by changing the value of $$\tau_{2}$$. Table 3, shows the values of the objective functional in both controlled and uncontrolled cases by using IOCM with $$\tau_{1}=0.5$$, $$\tau_{2}=0.2$$. Moreover, if we compute the value of the objective functional with $$\tau_{1}=0.2$$, $$\tau_{2}=0.5$$, we obtian the same result of the previous case. Fig. 3, shows that, the quantity $$S(t)+L_{s}(t)+L_{m}(t)+L_{x}(t)+I_{s}(t)+I_{m}(t)+I_{x}(t)+R(t)/N$$ is a constant in time at $$\alpha=1,$$ in controlled case, $$0\leq u_{1},u_{2}\leq 1$$. Fig. 4, shows the behavior of the approximate solutions for state variables $$I_{m},$$$$I_{x}$$ and $$R(t)$$ with $$\alpha=1$$ and $$\tau_{1}=\tau_{2}=\tau=0.1,$$ in the controlled and uncontrolled cases by using IOCM and GEM. We noted that, in the uncontrolled case, the delays in diagnosis and commencement of treatment to the individuals of $$I_{m}$$ and $$I_{x}$$ causing a shortage in the number of individuals of $$R(t)$$, but in the controlled case the number of $$I_{m}$$ and $$I_{x}$$ are decreasing with time as we expected and the number of $$R(t)$$ is increasing. Fig. 5, shows the behavior of the approximate solutions for state variables $$I_{m},$$$$I_{x}$$ with $$\alpha=0.97$$ and $$\tau=0.2$$ by using GEM and IOCM. We noted that, the results obtained by IOCM is better than the results obtained by GEM, where the number of $$I_{m}$$ and $$I_{x}$$ are decreasing with time in a rate faster than the results obtained by GEM. Fig. 6, shows the optimal controls $$u_{1}^{\ast}$$ and $$u_{2}^{\ast},$$ for $$\alpha=0.99$$, in a time units of years and $$\tau_{1}=\tau_{2}=\tau$$ by using ICOM. We noted that, when the value of $$\tau$$ is increasing the value of the objective functional is decreasing. Fig. 7, shows that the relationship between $$I_{m}(t)$$, $$I_{m}(t-\tau)$$ and $$I_{x}(t)$$, $$I_{x}(t-\tau)$$ in controlled case, where $$\tau=0.5$$, $$\alpha=0.99$$. We noted, that the number of $$I_{m}(t-\tau)$$ and $$I_{x}(t-\tau)$$ are more than the number of $$I_{m}(t)$$ and $$I_{x}(t)$$. Fig. 8, shows the behavior of the approximate solutions of the state variables $$I_{m}(t)$$, $$I_{x}(t)$$ and $$R(t)$$ with different value of $$\alpha$$, which are given to demonstrate how the fractional model is a generalization of the integer order model. Fig. 9, shows the values of $$u_{1}^{\ast}, u_{2}^{\ast}$$, in a time units of years by using ICOM with different values of $$\alpha$$ and in controlled case. Tables 4 and 5 show competitive studies between the two proposed methods to find the optimum values of the objective functional in both controlled and uncontrolled cases at different values of $$\alpha$$ and $$\tau_{1}=\tau_{2}=\tau$$. From the numerical results in these tables, it is found that IOCM can claim the minimum values in a rate faster than GEM and we can conclude that it is better than GEM, because, GEM can be failed to derive the optimum values for some values of $$\alpha$$, as given in Table 4. All results were obtained by using MATLAB (R2013a). on a computer machine with intel(R) core $$i3-3110M$$$$@$$$$2.40 GHz$$ and $$4 GB$$ RAM. Fig. 1. View largeDownload slide Plot of values of $$I_{m}(t-\tau_{1}), I_{x}(t-\tau_{2})$$ in a time units of years by using IOCM at different value of $$\tau_{1}$$ and fixed $$\tau_{2}=1$$, $$\alpha=0.99$$. Fig. 1. View largeDownload slide Plot of values of $$I_{m}(t-\tau_{1}), I_{x}(t-\tau_{2})$$ in a time units of years by using IOCM at different value of $$\tau_{1}$$ and fixed $$\tau_{2}=1$$, $$\alpha=0.99$$. Fig. 2. View largeDownload slide Plot of values of $$I_{m}(t-\tau_{1}), I_{x}(t-\tau_{2})$$ in a time units of years by using IOCM at different value of $$\tau_{2}$$ and fixed $$\tau_{1}=1$$, $$\alpha=0.99$$. Fig. 2. View largeDownload slide Plot of values of $$I_{m}(t-\tau_{1}), I_{x}(t-\tau_{2})$$ in a time units of years by using IOCM at different value of $$\tau_{2}$$ and fixed $$\tau_{1}=1$$, $$\alpha=0.99$$. Fig. 3. View largeDownload slide Plot of $$S(t)+L_{s}(t)+L_{m}(t)+L_{x}(t)+I_{s}(t)+I_{m}(t)+I_{x}(t)+R(t)/N$$ versus $$t$$ in years and $$\alpha=1$$, $$\tau=0.1,$$ using GEM. Fig. 3. View largeDownload slide Plot of $$S(t)+L_{s}(t)+L_{m}(t)+L_{x}(t)+I_{s}(t)+I_{m}(t)+I_{x}(t)+R(t)/N$$ versus $$t$$ in years and $$\alpha=1$$, $$\tau=0.1,$$ using GEM. Fig. 4. View largeDownload slide Numerical simulations of the state variable $$I_{m}$$, $$I_{x}$$ and $$R(t)$$ for the controlled case when $$0\leq u_{1}, u_{2}\leq 1$$, compared with the uncontrolled case when $$u_{1}= u_{2}=0$$ and $$\alpha=1$$, $$\tau=0.1,$$ using GEM and IOCM. Fig. 4. View largeDownload slide Numerical simulations of the state variable $$I_{m}$$, $$I_{x}$$ and $$R(t)$$ for the controlled case when $$0\leq u_{1}, u_{2}\leq 1$$, compared with the uncontrolled case when $$u_{1}= u_{2}=0$$ and $$\alpha=1$$, $$\tau=0.1,$$ using GEM and IOCM. Fig. 5. View largeDownload slide Numerical simulations of the state variables $$I_{m}$$, $$I_{x}$$ for the controlled case when $$0\leq u_{1}, u_{2}\leq 1$$, compared with the uncontrolled case when $$u_{1}= u_{2}=0$$ and $$\alpha=0.97$$, $$\tau=0.2,$$ using GEM and IOCM. Fig. 5. View largeDownload slide Numerical simulations of the state variables $$I_{m}$$, $$I_{x}$$ for the controlled case when $$0\leq u_{1}, u_{2}\leq 1$$, compared with the uncontrolled case when $$u_{1}= u_{2}=0$$ and $$\alpha=0.97$$, $$\tau=0.2,$$ using GEM and IOCM. Fig. 6. View largeDownload slide Plot of values of $$u_{1}^{\ast}, u_{2}^{\ast}$$ for $$\alpha=0.99$$, in a time units of years by using IOCM and different $$\tau$$. Fig. 6. View largeDownload slide Plot of values of $$u_{1}^{\ast}, u_{2}^{\ast}$$ for $$\alpha=0.99$$, in a time units of years by using IOCM and different $$\tau$$. Fig. 7. View largeDownload slide The relationship between $$I_{m}(t)$$, $$I_{m}(t-\tau)$$ and $$I_{x}(t)$$, $$I_{x}(t-\tau)$$ in controlled case, where $$\tau=0.5$$, $$\alpha=0.99$$. Fig. 7. View largeDownload slide The relationship between $$I_{m}(t)$$, $$I_{m}(t-\tau)$$ and $$I_{x}(t)$$, $$I_{x}(t-\tau)$$ in controlled case, where $$\tau=0.5$$, $$\alpha=0.99$$. Fig. 8. View largeDownload slide Numerical simulations of the state variable $$I_{m}, I_{x}$$ and $$R$$ for the controlled case with different values of $$\alpha$$, $$\tau=0.1$$, using IOCM. Fig. 8. View largeDownload slide Numerical simulations of the state variable $$I_{m}, I_{x}$$ and $$R$$ for the controlled case with different values of $$\alpha$$, $$\tau=0.1$$, using IOCM. Fig. 9. View largeDownload slide Plot of values of $$u_{1}^{\ast}, u_{2}^{\ast}$$ in a time units of years by using IOCM at different $$\alpha$$ and $$\tau=0.1$$. Fig. 9. View largeDownload slide Plot of values of $$u_{1}^{\ast}, u_{2}^{\ast}$$ in a time units of years by using IOCM at different $$\alpha$$ and $$\tau=0.1$$. Table 3 The values of objective functional in both controlled and uncontrolled cases by using IOCM with $$\tau_{1}=0.5$$, $$\tau_{2}=0.2$$ $$\alpha$$ $$J(u_{1}^{\ast},u_{2}^{\ast})$$ with control $$J(u_{1}^{\ast},u_{2}^{\ast})$$ without control 1 35801 53734 0.98 27062 49283 0.90 10392 35762 0.80 3184.1 26817.3 0.70 1054.2 17813.6 0.60 4580.2 10888.2 0.50 230.5 6898.6 0.40 118.4 5065.9 $$\alpha$$ $$J(u_{1}^{\ast},u_{2}^{\ast})$$ with control $$J(u_{1}^{\ast},u_{2}^{\ast})$$ without control 1 35801 53734 0.98 27062 49283 0.90 10392 35762 0.80 3184.1 26817.3 0.70 1054.2 17813.6 0.60 4580.2 10888.2 0.50 230.5 6898.6 0.40 118.4 5065.9 Table 3 The values of objective functional in both controlled and uncontrolled cases by using IOCM with $$\tau_{1}=0.5$$, $$\tau_{2}=0.2$$ $$\alpha$$ $$J(u_{1}^{\ast},u_{2}^{\ast})$$ with control $$J(u_{1}^{\ast},u_{2}^{\ast})$$ without control 1 35801 53734 0.98 27062 49283 0.90 10392 35762 0.80 3184.1 26817.3 0.70 1054.2 17813.6 0.60 4580.2 10888.2 0.50 230.5 6898.6 0.40 118.4 5065.9 $$\alpha$$ $$J(u_{1}^{\ast},u_{2}^{\ast})$$ with control $$J(u_{1}^{\ast},u_{2}^{\ast})$$ without control 1 35801 53734 0.98 27062 49283 0.90 10392 35762 0.80 3184.1 26817.3 0.70 1054.2 17813.6 0.60 4580.2 10888.2 0.50 230.5 6898.6 0.40 118.4 5065.9 Table 4 The values of objective functional in both controlled and uncontrolled cases with $$\tau_{1}=\tau_{2}=0.3$$ $$\alpha$$ Methods $$J(u_{1}^{\ast},u_{2}^{\ast})$$ with control $$J(u_{1}^{\ast},u_{2}^{\ast})$$ without control 1 GEM 33496 56227 IOCM 33496 0.98 GEM 27165 45941 IOCM 23243 0.90 GEM 10150 15402 IOCM 7269 0.80 GEM 5146.1 3179.8 IOCM 2911.3 0.70 GEM Non 995.87 IOCM 966.13 0.60 GEM Non 448.002 IOCM 380.555 0.50 GEM Non 232.32 IOCM 198.1596 $$\alpha$$ Methods $$J(u_{1}^{\ast},u_{2}^{\ast})$$ with control $$J(u_{1}^{\ast},u_{2}^{\ast})$$ without control 1 GEM 33496 56227 IOCM 33496 0.98 GEM 27165 45941 IOCM 23243 0.90 GEM 10150 15402 IOCM 7269 0.80 GEM 5146.1 3179.8 IOCM 2911.3 0.70 GEM Non 995.87 IOCM 966.13 0.60 GEM Non 448.002 IOCM 380.555 0.50 GEM Non 232.32 IOCM 198.1596 Table 4 The values of objective functional in both controlled and uncontrolled cases with $$\tau_{1}=\tau_{2}=0.3$$ $$\alpha$$ Methods $$J(u_{1}^{\ast},u_{2}^{\ast})$$ with control $$J(u_{1}^{\ast},u_{2}^{\ast})$$ without control 1 GEM 33496 56227 IOCM 33496 0.98 GEM 27165 45941 IOCM 23243 0.90 GEM 10150 15402 IOCM 7269 0.80 GEM 5146.1 3179.8 IOCM 2911.3 0.70 GEM Non 995.87 IOCM 966.13 0.60 GEM Non 448.002 IOCM 380.555 0.50 GEM Non 232.32 IOCM 198.1596 $$\alpha$$ Methods $$J(u_{1}^{\ast},u_{2}^{\ast})$$ with control $$J(u_{1}^{\ast},u_{2}^{\ast})$$ without control 1 GEM 33496 56227 IOCM 33496 0.98 GEM 27165 45941 IOCM 23243 0.90 GEM 10150 15402 IOCM 7269 0.80 GEM 5146.1 3179.8 IOCM 2911.3 0.70 GEM Non 995.87 IOCM 966.13 0.60 GEM Non 448.002 IOCM 380.555 0.50 GEM Non 232.32 IOCM 198.1596 Table 5 The values of objective functional in both controlled and uncontrolled cases with $$\tau_{1}=\tau_{2}=0$$ $$\alpha$$ Method $$J(u_{1}^{\ast},u_{2}^{\ast})$$ with control $$J(u_{1}^{\ast},u_{2}^{\ast})$$ without control 1 GEM 36626 57704 IOCM 36626 0.98 GEM 30614 48921 IOCM 25441 0.95 GEM 22791 36854 IOCM 15463 0.90 GEM 13383 21322 IOCM 6672.1 $$\alpha$$ Method $$J(u_{1}^{\ast},u_{2}^{\ast})$$ with control $$J(u_{1}^{\ast},u_{2}^{\ast})$$ without control 1 GEM 36626 57704 IOCM 36626 0.98 GEM 30614 48921 IOCM 25441 0.95 GEM 22791 36854 IOCM 15463 0.90 GEM 13383 21322 IOCM 6672.1 Table 5 The values of objective functional in both controlled and uncontrolled cases with $$\tau_{1}=\tau_{2}=0$$ $$\alpha$$ Method $$J(u_{1}^{\ast},u_{2}^{\ast})$$ with control $$J(u_{1}^{\ast},u_{2}^{\ast})$$ without control 1 GEM 36626 57704 IOCM 36626 0.98 GEM 30614 48921 IOCM 25441 0.95 GEM 22791 36854 IOCM 15463 0.90 GEM 13383 21322 IOCM 6672.1 $$\alpha$$ Method $$J(u_{1}^{\ast},u_{2}^{\ast})$$ with control $$J(u_{1}^{\ast},u_{2}^{\ast})$$ without control 1 GEM 36626 57704 IOCM 36626 0.98 GEM 30614 48921 IOCM 25441 0.95 GEM 22791 36854 IOCM 15463 0.90 GEM 13383 21322 IOCM 6672.1 Conclusion In this article, numerical solutions for the optimal control of fractional multi-strain TB model with two time delays are presented. Modified parameters are introduced to account for the fractional order model. We introduce two delays in the proposed model, representing the time delay in the diagnosis and commencement of treatment of individuals with active TB infection in two and three strains. Two control functions $$u_{1},$$$$u_{2}$$ are introduced. The controls $$u_{1},$$$$u_{2}$$ represent the effort that prevents the failure of treatment in active TB infectious individuals $$I_{m},$$$$I_{x},$$ e.g., supervising the patients, helping them to take the TB medications regularly and to complete the TB treatment. Necessary optimality conditions are derived. Two numerical methods; IOCM and GEM are introduced to study the behavior of the model problem. Comparative studies are implemented. Some figures are given to demonstrate how the fractional delay model is a generalization of the integer order model. It can be concluded from the numerical results presented in this article that, IOCM is better than GEM and it can be applied simply and effectively to solve optimal control of fractional delay multi-strain TB. Acknowledgements The authors would like to thanks the anonymous reviewers very much for their positive comments, careful reading, and useful suggestions on improving this article. References Abta, A. Laarabi, H. & Alaoui, H. T. ( 2014 ) The Hopf bifurcation analysis and optimal control of a delayed SIR epidemic model. Int. J. Anal. , 2014 , 1 – 10 . Google Scholar CrossRef Search ADS Agarwal, O. P. ( 2008 ) A formulation and numerical scheme for fractional optimal control problems. J. Vib. Control , 14 , 1291 – 1299 . Google Scholar CrossRef Search ADS Agarwal, P. & Choi, J. ( 2016 ) Fractional calculus operators and their image formulas. J. Korean Math. Soc. , 53 , 1183 – 1210 . Google Scholar CrossRef Search ADS Aldila, D. , Gãtz, T. & Soewono, E. ( 2013 ) An optimal control problem arising from a dengue disease transmission model. Math. Biosci ., 242 , 9 – 16 . Google Scholar CrossRef Search ADS PubMed Allen, L. J. S. ( 2007 ) An Introduction to Mathematical Biology . NJ : Prentice Hall . Arenas, A. J. , Gonzãlez-Parra, G. & Chen-Charpentierc, B. M. ( 2016 ) Construction of nonstandard finite difference schemes for the SI and SIR epidemic models of fractional order. Math. Comput. Simul. , 121 , 48 – 63 . Google Scholar CrossRef Search ADS Arino, J. & Soliman, I. A. ( 2015 ) A model for the spread of tuberculosis with drug-sensitive and emerging multidrug-resistant and extensively drug resistant strains. Math. Comput. Model. , Wiley, 1 – 120 , doi:10.1002/9781118853887.ch5 . Baleanu, D. , Diethelm, K. , Scalas, E. & Trujillo, J. J. ( 2012 ) Fractional Calculus Models and Numerical Methods. Series on Complexity, Nonlinearity and Chaos . Hackensack, NJ : World Scientific Publishing Co. Pte. Ltd . Google Scholar CrossRef Search ADS Bellen, A. & Zennaro, M. ( 2003 ) Numerical Methods for Delay Differential Equations, Oxford Science Publications . Oxford, UK : Clarendon Press . Google Scholar CrossRef Search ADS Bhalekar, S. & Daftardar-Gejji, V. ( 2011 ) A predictor-corrector scheme for solving non-linear delay differential equations of fractional order. J. Fract. Calc. Appl. , 5 , 1 – 9 . Bocharov, G. & A Rihan, F. ( 2000 ) Numerical modelling in biosciences using delay differential equations. J. Comput. Appl. Math. , 125 , 183 – 199 . Google Scholar CrossRef Search ADS Boklund, A. , Toft, N. , Alban, L. & Uttenthal, A. ( 2009 ) Comparing the epidemiological and economic effects of control strategies against classical swine fever in Denmark. Prev. Vet. Med. , 90 , 180 – 193 . Google Scholar CrossRef Search ADS PubMed Castillo-ćhavez, C. & Feng, Z. ( 1997 ) To treat or not to treat: the case of tuberculosis. J. Math. Biol. , 35 , 629 – 656 . Google Scholar CrossRef Search ADS PubMed Cohen, T. , & Murray, M. ( 2004 ) Modeling epidemics of multidrug-resistant M. tuberculosis of heterogeneous fitness. Nat. Med. , 10 , 1117 – 1121 . Google Scholar CrossRef Search ADS PubMed Cooke, R. & Bellman, K. L. ( 1963 ) Differential-Difference Equations . New York, NY : Academic Press. Google Scholar CrossRef Search ADS Denysiuk, R. , Silva, C. J. & Torres, D. F. M. ( 2014 ) Multiobjective approach to optimal control for a tuberculosis model. Optim. Meth. Software, 2015 , 1029 – 4937 . Driver, R. ( 1977 ) Ordinary and Delay Differential Equations . Berlin, Germany : Springer . Google Scholar CrossRef Search ADS El-Saka, H. A. A. ( 2013 ) The fractional-order SIR and SIRs epidemic models with variable population size. Math. Sci. Lett., 2 , 1 – 6 . Google Scholar CrossRef Search ADS Fowler, A. C. , Mackey, M. C. ( 2002 ) Relaxation oscillations in a class of delay differential equations. SIAM J. Appl. Math. 63 : 299 – 323 , (2002). Google Scholar CrossRef Search ADS Hale, J. ( 1977 ) Theory of Functional Differential Equations . New York, NY, USA : Springer . Google Scholar CrossRef Search ADS Hattaf, K. , Yousfi, N. ( 2012 ) Optimal control of a delayed HIV infection model with immune response using an efficient numerical method. Int. Sch. Res. Netw. , 2012 , 1 – 7 . Kuang, Y. ( 1993 ) Delay Differential Equations with Applications in Population Dynamics . San Diego : Academic Press . Laarabi, H. , Abta, A. & Hattaf, K. ( 2015 ) Optimal control of a delayed sirs epidemic model with vaccination and treatment, Acta Biotheor. , 63 , 87 – 97 . Google Scholar CrossRef Search ADS PubMed Matignon, D. ( 1996 ) Stability results for fractional differential equations with applications to control processing, Computational Engineering in Systems and Application. Multiconference , Vol. 2. Lille, France : IMACS, IEEE-SMC , pp. 963 – 968 . Nagy, A. M. , Sweilam, N. H. ( 2014 ) An efficient method for solving fractional Hodgkin Huxley model. Phys. Lett. A, 378 , 1980 – 1984 . Google Scholar CrossRef Search ADS Nelson, P. W. , Perelson, A. S. ( 2002 ) Mathematical analysis of delay differential equation models of HIV-1 infection. Math. Biosci. , 179 , 73 – 94 . Google Scholar CrossRef Search ADS PubMed Odibat, Z. & Momani, S. ( 2008 ) An algorithm for the numerical solution of differential equations of fractional order. J. Appl. Math. Inform ., 26 , 15 – 27 . Pimenov, V. G. & Tashirova, E. E. ( 2013 ) Numerical methods for solving a hereditary equation of hyperbolic type. Proc. Steklov Inst. Math. , 281 , 126 – 136 . Google Scholar CrossRef Search ADS Podlubny, I. ( 1999 ) Fractional Differential Equations . New York : Academic Press . Rihanm, F. A. ( 2000 ) Numerical Treatment of Delay Dfferential Equation in Bioscience . Manchester, UK : University of Manchester . Rihan, F. A. , Hashish, A. , Al-Maskari, F. , Hussein, M. S. and Ahmed, E. ( 2016 ) Dynamics of tumor-immune system with fractional-order. J. Tumor Res. 2 , 109 . Rihan, F. A. , Lakshmanan, S. , Hashish, A. H. , Rakkiyappan, R. & Ahmed, E. ( 2015 ) Fractional-order delayed predator-prey systems with Holling type-II functional response, Nonlinear Dynamics. An International. J. Nonlinear Dynam. Chaos Eng. Syst. , 80 , 777 – 789 . Google Scholar CrossRef Search ADS Salahshour, S. , Ahmadian, A. , Senu, N. , Baleanu, D. & Agarwal, P. ( 2015 ) On analytical solutions of the fractional differential equation with uncertainty: application to the basset problem. Entropy , 17 , 885 – 902 . Google Scholar CrossRef Search ADS Small, P. M. and Fujiwara, P. I. ( 2001 ) Management of Tuberculosis in the United States. N. Engl. J. Med. , 345 , 189 – 200 . Google Scholar CrossRef Search ADS PubMed Smith, H. ( 2011 ) An Introduction To Delay Differential Equations with Applications To The Life Sciences . New York : Springer . Google Scholar CrossRef Search ADS Silva, C. J. , Maurer, H. & Torres, D. F. ( 2017 ) Optimal control of a tuberculosis model with state and control delays. Math. Biosci. Eng. , 14 . Sreeramareddy, C. T. , Panduru, K. V. , Menten, J. & Van den Ende, J. ( 2009 ) Time delays in diagnosis of pulmonary tuberculosis: a systematic review of literature. BMC Infect. Dis. , 9 . Styblo, K. ( 1978 ) State of art: epidemiology of tuberculosis. Bull. Int. Union Tuberc., 53 , 141 – 152 . Sweilam, N. H. , AL-Mekhlafi, S. M. ( 2016 ) Comparative study for multi-strain tuberculosis (TB) model of fractional order. J. Appl. Math. Inf. Sci. , 10 , 1403 – 1413 . Google Scholar CrossRef Search ADS Sweilam, N. H. , AL-Mekhlafi, S. M. ( 2017 ) Legendre spectral-collocation method for solving fractional optimal control of HIV infection of $$Cd4^{+}T$$ cells mathematical model. J. Defense Model. Simulat. Appl. Methodol. Technol. , 14 , 273 – 284 . Google Scholar CrossRef Search ADS Sweilam, N. H. & AL-Mekhlafi, S. M. ( 2016a ) Numerical study for multi-strain tuberculosis (TB) model of variable-order fractional derivatives. J. Adv. Res. , 7 , 271 – 283 . Google Scholar CrossRef Search ADS Sweilam, N. H. , AL-Mekhlafi, S. M. ( 2016b ) On the optimal control for fractional multi-strain TB model. Optim. Contr. Appl. Meth. , 37 , 1355 – 1374 . Google Scholar CrossRef Search ADS Sweilam, N. H. , Khader, M. M. & Adel, M. ( 2012a ) Numerical studies for fractional-order logistic differential equation with two different delays. J. Appl. Math. , 14 pages, Article ID 764894, doi:10.1155/2012/764894 . Sweilam, N. H. , Khader, M. M. & Adel, M. ( 2012b ) On the stability analysis of weighted average finite difference methods for fractional wave equation. Fract. Differ. Calc. , 2 , 17 – 29 . Google Scholar CrossRef Search ADS Sweilam, N. H. , Khader, M. M. & Nagy, A. M. ( 2011 ) Numerical solution of two-sided space-fractional wave equation using finite difference method. J. Comput. Appl. Math. 235 , 2832 – 2841 . Google Scholar CrossRef Search ADS Sweilam, N. H. , Soliman, I. A. & Al-Mekhlafi, S. M. ( 2017 ) Nonstandard finite difference method for solving the multi-strain TB model. J. Egypt. Math. Soc. , 25 , 129 – 138 . Google Scholar CrossRef Search ADS Tariboon, J. , Ntouyas, S. K. & Agarwal, P. ( 2015 ) New concepts of fractional quantum calculus and applications to impulsive fractional q-difference equations. Adv. Differ. Equ. , 2015 . Toman, K. ( 1979 ) Tuberculosis Case-Finding and Chemotherapy: Questions and Answers . Geneva : WHO . Uys, P. W. , Warren, M. & van Helden, P. D. ( 2007 ) A threshold value for the time delay to TB diagnosis. PLoS One , 2 . Van den Driessche, P. & Watmough, J. ( 2002 ) Reproduction numbers and subthreshold endemic equilibria for compartmental models of disease transmission. Math. Biosc., 180 , 29 – 48 . Google Scholar CrossRef Search ADS Wang, Z. , Huang, X. & Zhou, J. ( 2013 ) A numerical method for delayed fractional order differential equations: based on G-L Definition. Appl. Math. Inf. Sci. , 7 , 525 – 529 . Google Scholar CrossRef Search ADS World Health Organization ( 2012 ) Multidrug and Extensively Drug-Resistant TB (M/XDR-TB): 2012 Global Report on Surveillance and Response . World Health Organization : Geneva, Switzerland . World Health Organization ( 2014 ) Global Tuberculosis Report 2014 . Geneva : World Health Organization , http://www.who.int/tb/publications/global_report/en/. Zhang, X. , Agarwal, P. , Liu, Z. , Ding, W. & Ciancio, A. ( 2016 ) On the fractional differential equations with not instantaneous impulses. Open Phys., 14 , 676 – 684 . Zhou, H. , Yang, L. , Agarwal, P. ( 2017 ) Solvability for fractional p-Laplacian differential equations with multipoint boundary conditions at resonance on infinite interval. J. Appl. Math. Comput , 53 , 51 – 76 . Google Scholar CrossRef Search ADS © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved.

### Journal

IMA Journal of Mathematical Control and InformationOxford University Press

Published: Nov 9, 2017

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations