EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS Vol. 17, No. 4, 2024, 2706-2719 ISSN 1307-5543 – ejpam.com Published by New York Business Global Efficient Numerical Solutions for Breast Cancer Model Taghreed A. Assiri Mathematics Department, Faculty of Sciences, Umm Al-Qura University, Makkah, Saudi Arabia Abstract. We know that mathematical modeling is a serious agent to realize the dynamics of Breast Cancer development, spread, and educe new therapeutic approaches. The main target of this study is to use Simpson’s 1/3 rule as an efficient numerical scheme for integration to numerically treat the obtained fractional integral equations and reduce them to a set of algebraic equations. The fractional derivatives here are taken in the Caputo-Fabrizio (CF) sense. Particular assurance is located on elucidating the error analysis of the given scheme. The results acquired by implementing the Runge-Kutta method (RK4M) are compared to those obtained using the achieved results. The results explain that the implemented scheme offers a precise and effective agent for simulating this model. The primary benefit of the implemented method is that it relies on a small number of uncomplicated steps and does not have long-term effects. 2020 Mathematics Subject Classifications: 34A12, 41A30, 26C10, 47H10, 65N20 Key Words and Phrases: Breast Cancer, Caputo-Fabrizio fractional derivative, Numerical integration, RK4M 1. Introduction Breast Cancer (BC) is described and defined as a disease condition resulting from the uncontrolled proliferation of cells inside the breast tissue. Through the data supplied by the World Health Organization (WHO) on the global burden of cancer, we can find that the BC has the highest prevalence rate when assorted to other types of cancer [8]. Globally, through surveys conducted by WHO, which ranked BC in 2004 as the second most widespread form of cancer, we find that it represents a major potential threat to women, as it influenced approximately (8-9)% of women worldwide [1]. Despite all these meditations and numerous fulfillments, the exact reason of BC is still mysterious. It was also directly responsible for the deaths of 685,000 people in 2020, out of 2.3 million women affected, as indicated by the diagnosis of 7.8 million women during the previous 5-years. Breast cancer is most common in women after puberty and its incidence increases with age [1]. Based on all of the above considerations, we emphasize the need for an aggregate understanding of the epidemiology of BC and its effects on women’s health, as it is of DOI: https://doi.org/10.29020/nybg.ejpam.v17i4.5387 Email address: taassiri@uqu.edu.sa (T. A. Assiri) https://www.ejpam.com 2706 Copyright: © 2024 The Author(s). (CC BY-NC 4.0) T. A. Assiri / Eur. J. Pure Appl. Math, 17 (4) (2024), 2706-2719 2707 great importance in improving effective preventive and therapeutic methods all over the world ([11], [21]). Scientists are currently studying more effective fractional operators ([5]-[23]). In order to address the issue of singularity and achieve accurate and reliable modeling outcomes, a more effective CF fractional derivative has been developed. This derivative incorporates a non-singular kernel, as proposed by Caputo and Fabrizio, leading to improved efficiency and robustness in recent years. The use of Laplace transformation to convert it to integer power is regarded as a constructive method. Consequently, in some scenarios, we can readily compute the precise solution. The citation for this information is [6]. The paper [16] explores the examination of fractional operators and presents novel characteristics associated with them. Fractional derivatives have acquired increasing leverage in the last few years, as they have greatly contributed to presenting the most difficult mathematical models in a clearer and more comprehensive manner, which in turn has helped researchers to obtain more accurate solutions to such models by applying appropriate numerical and approximate methods. The fundamental focus of the research [13] is the monotonicity theorem for Fractional Differential Equations (FDEs) in the Caputo form. A study [18] has devised a highly effective numerical technique, known as the iterative reproducing kernel algorithm, for computing approximate solutions to fractional Riccati and Bernoulli DEs while considering the CF-derivative. Furthermore, a work conducted in [7] has examined the intricate dynamics of the Omicron variant of COVID-19 by employing CF-fractional operators, and others [25]. In addition, a numerical method is created that includes an exponential law kernel to analyze and model the spread of the infection. Many mathematicians find it challenging to create numerical and approximate solu- tions for the FDEs ([3]-[14]). The Adams-Bashforth method, incorporating the CF op- erator, is formulated in [22]. This method involves three steps and can be used to solve both linear/nonlinear FDEs. Additionally, it possesses diverse uses in resolving chaotic systems with fractional orders. The authors of [12] have developed a trapezoidal strategy to solve the FDEs efficiently. This scheme utilizes the CF operator and achieves a con- vergence order of two. Additionally, the convergence and stability of this technique have been thoroughly investigated. Inspired by this research, we devised Simpson’s 1/3 scheme for solving FDEs. This method achieves a high level of accuracy, with an order of four, as detailed in our work. The proposed fractional Simpson’s 1/3 approach offers superior accuracy compared to current methods and is straightforward to implement. The novelty in the current research is the attempt to reach numerical simulations (through a good numerical method) to study the important system under consideration which is of interest to many researchers, so we presented by shedding some light on the convergence and calculating the resulting error as well as comparisons with the same method but in a less accurate case and finally the effect of the parameters in the system on the behavior of the solution to provide recommendations that can be used by those interested in studying this model medically or industrially. This work is devoted to giving a numerical solution for the fractional BC model. This is by using the Simpson’s 1 3 rule for CF-fractional integral. The outline of the paper is given as follows: Section 2, presents some basic concepts T. A. Assiri / Eur. J. Pure Appl. Math, 17 (4) (2024), 2706-2719 2708 of fractional calculus in sense of the Caputo and Fabrizio. Section 3, describes the breast cancer system in its fractional form. Section 4, derives the Simpson’s-13 rule for CF- fractional integral. Section 5, outlines a numerical implementation for solving the model under study. Section 6, presents a numerical simulation for the model under study. Section 7, gives the conclusions and remarks. 2. Preliminaries Mathematical models using more accurate fractional derivatives (FDs) have been used by taking advantage of a non-singular kernel. This approach enhances the system’s ability to accurately represent and capture memory effects. In 2015, Caputo and Fabrizio suc- cessfully introduced CF-FD by substituting the singular kernel (t − τ)−γ with e ( −γ(t−τ) 1−γ ) in the Caputo derivative [6]. Here, we will provide a succinct overview of fundamental definitions pertaining to fractional calculus involving a non-singular kernel. Definition 1. For ψ(t) ∈ H1(0, a), 0 < γ < 1. Then the CF-FD CFDγψ(t) and CF fractional integral CF Iγψ(t), respectively are defined by: CFDγψ(t) := 1 1− γ ∫ t 0 Exp [ − γ 1− γ (t− τ) ] ψ̇(τ)dτ, CF Iγψ(t) := (1− γ)ψ(t) + γ ∫ t 0 ψ(τ)dτ. (1) 3. Formulation of the model As mentioned above, the mathematical modeling of the BC or any other biological phenomena has been an important tool for understanding the dynamics of tumor growth in the treatment process and solving epidemiological problems. This is evident from the studies previously done in [19]. Despite the large number of studies, none of the proposed mathematical models included a diet (ketogenic diet). Therefore, Oak et al. developed the model in [20], to include some of the control factors such as immune booster, anticancer drug, and ketogenic diet, in order to confirm that there is an interaction between cells due to the variation in the tumor cell DNA [9]. We study the following fractional forms of breast cancer ([2], [24]): CFDνψ1(t) = λ1 ψ1(t) (p1 − β1ψ1(t)− α1 ψ2(t))− (1− p)θ1 ψ1(t)ψ4(t), CFDνψ2(t) = λ2 ψ2(t) (p2d− β2ψ2(t)− α2ψ3(t))− κψ2(t) + (1− p)θ1ψ1(t) + ψ2(t)ψ4(t), CFDνψ3(t) = σρ+ λ1 ψ3(t) (p3 − β3 − α3ψ2(t))− (1− p)θ2ψ3(t)ψ4(t), CFDνψ4(t) = α4ψ4(t) + ϱ(1− p), (2) T. A. Assiri / Eur. J. Pure Appl. Math, 17 (4) (2024), 2706-2719 2709 the corresponding I.Cs of this model are taken as follows: ψ1(0) = ψ̂0 1, ψ2(0) = ψ̂0 2, ψ3(0) = ψ̂0 3, ψ4(0) = ψ̂0 4. The description of the meaning of the included parameters (∈ R+) of the system (2), will give in the Table 1 [20]. The stability analysis, equilibrium points, existence, & uniqueness, of the system under consideration are given in [24]. Table 1: The description of the included variables of the system (2). Symbol Description ψ1(t) Normal cell population ψ2(t) Luminal type tumor cells ψ3(t) Class of immune response ψ4(t) Estrogen compartment λ1 Growth rate of ψ1 λ2 Growth rate of ψ2 pi Carrying capacity of ψi, i = 1, 2, 3 (1− p) Effectiveness of anti-cancer drugs ϱ Process of constantly replenishing excess estrogen α1 Inhibition rate of ψ1(t) α2 Rate of the effectiveness of the immune system to the tumor cells α3 Rate of interaction between ψ2(t) and ψ3(t) θ1 Tumor formation rate resulting from DNA mutation caused by the presence of excess estrogen θ2 Immune suppression rate βi Logistic rate of ψi, i = 1, 2, 3 d Ketogenic diet 4. Derivation Simpson’s-1/3 rule for CF-fractional integral This section presents the formulation of the fractional Simpson’s-1/3 rule (FSR) for solving CF-FDEs [4]. Considering the subsequent γ-order FDE: CFDγu(t) = f(u(t)), u(0) = u0, (3) whereas f ∈ C[0, T ] and satisfies the Lipschitz condition: |f (u(t1))− f (u(t2)) | ≤ ϵ |u(t1)− u(t2)|, ϵ > 0. (4) Applying CF-fractional integral operator on the IVP (3) and applying Proposition 3 in [2] and formula (1), we get: u(t) = u0 + CF Iγf (u(t)) = u0 + (1− γ)f(u(t)) + γ ∫ t 0 f(u(s))ds. (5) T. A. Assiri / Eur. J. Pure Appl. Math, 17 (4) (2024), 2706-2719 2710 Theorem 1. [17] Let us assume a continuous function f : [0, T ] × R → R with γ ∈ (0, 1) that satisfies the criterion (4). Then the IVP (3) has a unique solution on C[0, T ] under the condition: (2(1− γ) + 2γ T )ϵ (2− γ) < 1. Initially, we will employ a quadratic polynomial P2 to estimate the integral function f in Equation (5). The function will be assessed at t0, t1, and t2, where t0 is less than t1, and t1 is less than t2. The interval is partitioned into two subintervals, denoted as t1− t0 = t2− t1 = h, resulting in a combined width of 2h. The integration of the quadratic polynomial P2 can be computed as follows: I2f(u(t)) = ∫ b a f(u(t))dt ≈ ∫ t2 t0 P2(u(t))dt = ∫ t2 t0  2∑ j=0 Lj(t) f (u (tj))  dt, where the second-order Lagrange polynomials L0(t), L1(t), and L2(t) are defined as fol- lows: Lj(t) = 2∏ i=0, i ̸= j (t− ti) (tj − ti) . Integrating the first interpolant function L0(t), by taking h = t2−t0 2 and substituting ”t = s+ t0”, gives us:∫ t2 t0 L0(t)dt = 1 2h2 ∫ t0+2h t0 (t− t1) (t− t2) dt = 1 2h2 ∫ 2h 0 (s+ t0 − t2) (s+ t0 − t1) ds = 1 2h2 ∫ 2h 0 (s− 2h)(s− h)ds = h 3 . After making some simplifications to the rest of the terms, we have the following: I2(f) = h 3 [f (u (t0)) + 4f (u (t1)) + f (u (t2))] . Using Eq.(5), we get: u (tn) = u0 + (1− γ)f (u (tn)) + h 3 [f (u (t0)) + 4f (u (t1)) + f (u (t2))] , n = 0, 1, 2. To improve the accuracy of numerical integration, we partition [a, b] to n sub-intervals as follows: For any even number n ≥ 2, we establish the following definitions: h = b− a n = tk+1 − tk, k = 0, 1, 2, . . . , n. T. A. Assiri / Eur. J. Pure Appl. Math, 17 (4) (2024), 2706-2719 2711 Now, by using the quadrature rule for each pair of subintervals and implementing the Simpson’s-1/3 rule to each of [ t2k, t2(k+1) ] , k = 0, 1, 2, . . . , n−2 2 , we can appoint the following formula: In(f) = n−2 2∑ k=0 ∫ t2k+2 t2k f(u(t))dt = n−2 2∑ k=0 ( h 3 [f (u (t2k)) + 4f (u (t2k+1)) + f (u (t2k+2))] ) . By referring un as the approximate solution of u (tn) and using Eq.(5), we can write the FSR for the CF-FDE (3), for n = 0[1](m− 1): un+1 = u0+(1−γ)f (un+1)+γ h 3 f (u (t0)) + 4 n∑ i=2,4,6 f (u (ti)) + 2 n−1∑ j=1,3,5 f (u (tj)) + f (u (tn+1))  . This formula can be rewritten in a compact form as follows: un+1 = u0 + (1− γ)f (un+1) + γh n+1∑ r=0 ξrf (ur) , n = 0, 1, 2, . . . ,m− 1, (6) where ξr are the weights of the FSR and defined as: ξr =  1/3, r = 0, n+ 1, 2/3, r = 1, 3, 5, . . . , 4/3, r = 2, 4, 6, . . . . Lemma 1. [4] Suppose that f(u(t)) ∈ C4([a, b]), then the error of numerical scheme (6) is estimated by: ∣∣∣∣∣ ∫ tn+1 t0 f(u(s))ds− γh n+1∑ i=0 ξif (u(ti)) ∣∣∣∣∣ ≤ h4, where ĉ = (b−a)f (4)(ζ) 180 , for some constant a < ζ < b, h = b−a n , and tk = a + hk, k = 0, 1, . . . , n+ 1. The stability and error analysis of the regulated numerical scheme are investigated and proved meanwhile the following theorems [4]. Theorem 2. The newly designed fractional numerical technique (6) exhibits conditional stability. Proof. The proof of this theorem can be found on [4]. T. A. Assiri / Eur. J. Pure Appl. Math, 17 (4) (2024), 2706-2719 2712 Theorem 3. The recently developed fractional numerical method exhibits conditional convergence of order 4, as stated in equation (6). ∥u (tn+1)− un+1∥ ≤ Ch4, where C = γĉ ch. Proof. The proof of this criterion can be found on [4]. 5. Numerical implementation Here in this section, we will try to treat the shortcomings of the existing numerical methods, which are represented by the slow convergence of most of them when solving this type of problem, which in turn leads to inaccurate approximations ([14], [15]). For this, we will use the FSR for numerical integration to calculate the resulting integral in the system of integral equations obtained from the same system of FDEs under study. Now, we will numerically treat the BC system in its fractional form by creating a numerical scheme of it. For this purpose, let us reformulate the system (2) in an operator form as follows: CFDν Ψ̄(t) = F(Ψ̄(t), t), (7) where Ψ̄(t) = [ψ1(t), ψ2(t), ψ3(t), ψ4(t)] T , F(Ψ̄(t), t) = [f1, f2, f3, f4] T , Ψ̄(0) = [ψ̂0 1, ψ̂ 0 2, ψ̂ 0 3, ψ̂ 0 4] T , (8) where each one of the functions fi(ψ1, ψ2, ψ3, ψ4, t), i = 1(1)4 is defined in the RHS of the four equations in (2), respectively. Applying CF-fractional integral operator on the model (7) and implementing Propo- sition 3 in [2] and formula (1), we get: Ψ̄(t) = Ψ̄(0) + CF Iν F(Ψ̄(t), t) = Ψ̄(0) + (1− ν)F(Ψ̄(t), t) + ν ∫ t 0 F(Ψ̄(s), s)ds. (9) Applying the derived FSR for the integration on the RHS of (9), to get the following numerical scheme as constructed in the formula (6): Ψ̄n+1 = Ψ̄(0) + (1− ν)F(Ψ̄n+1, tn+1) + ν h n+1∑ r=0 ξr F(Ψ̄r, tr), n = 0, 1, 2, . . . ,m− 1, (10) where the weights ξr, r = 0, 1, ..., n+ 1 of the FSR are defined in (6). So, the system given in (10) transforms to an algebraic equations system as follows: ψk,n+1 = ψk,0 + (1− ν)fk(ψ1,n+1, ψ2,n+1, ψ3,n+1, ψ4,n+1, tn+1) + ν h n+1∑ r=0 ξr fk(ψ1,r, ψ2,r, ψ3,r, ψ4,r, tr), k = 1(1)4, where the functions fk are defined in (8). T. A. Assiri / Eur. J. Pure Appl. Math, 17 (4) (2024), 2706-2719 2713 6. Experimental results Now, we are going to test the precision of the resulting numerical scheme by introducing numerical simulation on some cases in [0, 3] for the proposed model (2). The behavior of ψk(t), k = 1, 2, 3, 4 are presented in Figures 1-5 at different values of some parameters (ν, ϱ, p, m). We approach the system under study (2) with the following values for the parameters contained in it [10]: θ1 = 0.2, θ2 = 0.02, λ1 = 0.3, λ2 = 0.4, d = 0.5, α1 = 6× 10−8, α2 = 3× 10−7, α3 = 1×10−7, α4 = 0.97, ρ = 0.01, κ = 2, β1 = 0.15, β2 = 0.7, β3 = 0.1, p = 0.5, σ = 0.1, p1 = 0.1, p2 = 0.2, p3 = 0.3, ϱ = 0.8, with the initial conditions (I.Cs) ψi(0) = 0.2, i = 1, 2, 3, 4. Figures 1-5 display a numerical simulation for the system under investigation by applying the given procedure. (i) Figure 1 outlines the solution via various values of ν = 1.0, 0.95, 0.85, 0.75, at h = 0.01. (ii) Figure 2 gives the solution for distinct values of m = 5, 10, 15. (iii) Figure 3 shows the solution for distinct values of the estrogen source rate ϱ = 0.8, 1.4, 2.0, 2.6. (iv) Figure 4 depicts the effect of p (= 0.25, 0.5, 0.75, 1.0) on the approximate solution. (v) Figure 5 presents a comparison between the solution created by the obtained ap- proach with that numerical solution by implementing RK4M [14] with the same parameters and I.Cs. The behavior of the approximate solution is based on ν, m, ϱ, p, as shown in Figures 1-4, respectively. From Figures 3 and 4, we can confirm that the behavior of the solution consists with the nature effect of the parameters ϱ, and p, respectively. From Figures 2 and 5, we can point out that the proposed method has been well implemented to solve the model under study. Hence, we can verify that the expected behavior of the disease has been obtained, which means that we have provided an evident simulation of the proposed system that can be used by the relevant authorities to treat this deadly cancer. T. A. Assiri / Eur. J. Pure Appl. Math, 17 (4) (2024), 2706-2719 2714 Figure 1. The solution ψi(t) via various values of ν. Figure 2. The solution ψi(t) via various values of m. T. A. Assiri / Eur. J. Pure Appl. Math, 17 (4) (2024), 2706-2719 2715 Figure 3. The solution ψi(t) via various values of ϱ. T. A. Assiri / Eur. J. Pure Appl. Math, 17 (4) (2024), 2706-2719 2716 Figure 4. The solution ψi(t) via various values of p. Figure 5. Comparison the solution obtained by proposed method and RK4M. REFERENCES 2717 7. Conclusions and remarks This study aims to utilize an efficient and accurate method to gain numerical solutions for the CF-fractional BC mathematical system. Simpson’s rule 1/3 was applied in its fractional form in computing the resulting integral within the system of fractional integral equations corresponding to the FDEs expressing mathematically the BC model to achieve fourth-order accuracy for the resulting solutions. This study utilized several values of the fractional order ν and ϱ, m, p to derive solutions for the model being examined. Also, we have determined that the proposed approach is remarkably effective in analyzing this system. Furthermore, reducing the value of h allows us to control the accuracy of the numerical solution. From the obtained solutions, we can confirm that the offered approach is surprisingly successful in simulating the BCmodel, as well as demonstrating the accuracy and computational effectiveness of this method. Finally, the present study may contribute to providing more robust physical explanations for future theoretical and computational studies on the same topic. References [1] Breast cancer. https://www.who.int/news-room/fact-sheets/detail/ breast-cancer. Accessed: 15 July 2023. [2] T. Abdeljawad and D. Baleanu. On fractional derivatives with exponential kernel and their discrete versions. Rep. Math. Phys., 80:11–27, 2017. [3] M. Adel, M. M. Khader, T. A. Assiri, and W. Kaleel. Numerical simulation for COVID-19 model using a multidomain spectral relaxation technique. Symmetry, 15:1–15, 2023. [4] S. Arshad, I. Saleem, A. Akgul, J. Huang, Y. Tang, and S. M. Eldin. A novel numerical method for solving the Caputo-Fabrizio fractional differential equation. AIMS Mathematics, 8:9535–9556, 2022. [5] A. Ait Brahim, J. El Ghordaf, A. El Hajaji, K. Hilal, and J. E. N. Valdes. A com- parative analysis of conformable, non-conformable, Riemann-Liouville, and Caputo fractional derivatives. Eur. J. of Pure and Appl. Maths., 17:1842–1854, 2024. [6] M. Caputo and M. Fabrizio. A new definition of fractional derivative without singular kernel. Progr. Fract. Differ. Appl., 1:73–85, 2015. [7] M. Farman, H. Besbes, K. S. Nisar, and M. Omri. Analysis and dynamical trans- mission of COVID-19 model by using Caputo-Fabrizio derivative. Alex. Eng. J., 66:597–606, 2023. [8] M. Fathoni, G. Gunardi, F. A. Kusumo, and S. H. Hutajulu. Mathematical model analysis of breast cancer stages with side effects on the Heart in chemotherapy pa- REFERENCES 2718 tients. In Proceedings of the AIP Conference Proceedings, Yogyakarta, Indonesia, 29 July–1 August 2019, 2019. AIP Publishing. [9] M. A. Hajji and Q. Al-Mdallal. Numerical simulations of a delay model for immune system-tumor interaction. Sultan Qaboos Univ. J. Sci., 23:19–31, 2018. [10] M. A. Hammad, I. H. Jebril, S. Alshorm, I. M. Batiha, and N. A. Hammad. Numer- ical solution for a fractional-order mathematical model of immune-chemotherapeutic treatment for breast cancer using the modified fractional formula. Int. J. Anal. Appl., 21:1–9, 2023. [11] M. Idrees, A. S. Alnahdi, and M. B. Jeelani. Mathematical modeling of breast cancer based on the Caputo-Fabrizio fractal-fractional derivative. Fractal and Fractional, 7:1–14, 2023. [12] A. Jajarmi, S. Arshad, and D. Baleanu. A new fractional modeling and control strategy for the outbreak of Dengue Fever. Physica A, 535:1–10, 2019. [13] T. Jin and X. Yang. Monotonicity theorem for the uncertain fractional differential equation and application to uncertain financial market. Math. Comput. Simulat., 190:03–221, 2021. [14] M. M. Khader and M. Adel. Modeling and numerical simulation for covering the frac- tional COVID-19 model using spectral collocation-optimization algorithms. Fractal Fract., 6:1–19, 2022. [15] M. M. Khader, N. H. Sweilam, A. M. S. Mahdy, and N. K. A. Moniem. Numerical simulation for the fractional SIRC model and influenza A. Appl. Math. Inf. Sci., 8:1029–1036, 2014. [16] J. Rong Loh, A. Isah, C. Phang, and Y. T. Toh. On the new properties of Caputo- Fabrizio operator and its application in deriving shifted Legendre operational matrix. Appl. Numer. Math., 132:138–153, 2018. [17] J. Losada and J. J. Nieto. Properties of a new fractional derivative without singular kernel. Progr. Fract. Differ. Appl., 1:87–92, 2015. [18] S. Momani, N. Djeddi, M. Al-Smadi, and S. Al-Omari. Numerical investigation for Caputo-Fabrizio fractional Riccati and Bernoulli equations using iterative reproducing kernel method. Appl. Numer. Math., 170:418–434, 2021. [19] C. Mufudza, S. Walter, and E. T. Chiyaka. Assessing the effects of estrogen on the dynamics of chemo-virotherapy cancer. Comput. Math. Methods Med., 473572:1–15, 2012. [20] S. I. Oke, M. B. Matadi, and S. S. Xulu. Optimal control analysis of a mathematical model for breast cancer. Math. Comput. Appl., 21:1–16, 2018. REFERENCES 2719 [21] F. Okose, M. T. Senel, and R. Habbireeh. Fractional-order mathematical modeling of cancer cells-cancer stem cells-immune system interaction with chemotherapy. Math. Model. Numer. Simul. Appl., 1:67–83, 2021. [22] K. M. Owolabi and A. Atangana. Analysis and application of new fractional Adams- Bashforth scheme with Caputo-Fabrizio derivative. Chaos Soliton Fract., 105:111– 119, 2017. [23] M. A. Rois, Fatmawati, and C. Alfiniyah. Optimal control of COVID-19 model with partial comorbidity sub-populations and two isolation treatments in Indonesia. Eur. J. of Pure and Appl. Maths., 16:523–537, 2023. [24] A . Yousef, F. Bozkurt, and T. Abdeljawad. Mathematical modeling of the immune- chemotherapeutic treatment of breast cancer under some control parameters. Adva. Differ. Equs., 696:1–25, 2020. [25] T. Zhang and Y. Li. Exponential Euler scheme of multi-delay Caputo-Fabrizio fractional-order differential equations. Appl. Math. Lett., 124:1–15, 2022.