Fractional Order Operator for Symmetric Analysis of Cancer Model on Stem Cells with Chemotherapy

Citation: Azeem, M.; Farman, M.; Akgül, A.; De la Sen, M. Fractional Order Operator for Symmetric Analysis of Cancer Model on Stem Cells with Chemotherapy. Symmetry 2023,15, 533. https://doi.org/ 10.3390/sym15020533 Academic Editors: Cemil Tunç, Jen-Chih Yao, Mouffak Benchohra and Ahmed M. A. El-Sayed Received: 12 January 2023 Revised: 13 February 2023 Accepted: 14 February 2023 Published: 16 February 2023 Copyright: © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). symmetry S S Article Fractional Order Operator for Symmetric Analysis of Cancer Model on Stem Cells with Chemotherapy Muhammad Azeem 1,† , Muhammad Farman 2,3,4,† , Ali Akgül 3,4,5,*,† and Manuel De la Sen 6,† 1Department of Mathematics and Statistics, The University of Lahore, Lahore 54590, Pakistan 2Institute of Mathematics, Khwaja Fareed University of Engineering and Information Technology, Rahim Yar Khan 64200, Pakistan 3Department of Computer Science and Mathematics, Lebanese American University, Beirut 5053, Lebanon 4Mathematics Research Center, Department of Mathematics, Near East University, Near East Boulevard, Nicosia, Mersin 99138, Turkey 5Department of Mathematics, Art and Science Faculty, Siirt University, Siirt 56100, Turkey 6Department of Electricity and Electronics, Institute of Research and Development of Processes, Faculty of Science and Technology, University of the Basque Country, 48940 Leioa, Spain *Correspondence: [email protected] † These authors contributed equally to this work. Abstract: Cancer is dangerous and one of the major diseases affecting normal human life. In this paper, a fractional-order cancer model with stem cells and chemotherapy is analyzed to check the effects of infection in individuals. The model is investigated by the Sumudu transform and a very effective numerical method. The positivity of solutions with the ABC operator of the proposed technique is verified. Fixed point theory is used to derive the existence and uniqueness of the solutions for the fractional order cancer system. Our derived solutions analyze the actual behavior and effect of cancer disease in the human body using different fractional values. Modern mathematical control with the fractional operator has many applications including the complex and crucial study of systems with symmetry. Symmetry analysis is a powerful tool that enables the user to construct numerical solutions of a given fractional differential equation in a fairly systematic way. Such an analysis will provide a better understanding to control the of cancer disease in the human body. Keywords: cancer model; existence; uniqueness; Sumudu transform; fractional operator; Atangana– Toufik method MSC: 37M05; 92B05 1. Introduction Cancer is considered to cause different diseases with different characteristics for a complex system. To better understand the dynamics of cancer, many researchers are attempting to use a variety of methodologies to investigate the relationships between immune cells and tumor cells [ 1 , 2 ]. Diseases are categorized; cancer is one of them which represents out-of-control cell growth. When two things arise, malignant tumors or extra hazards appear in the human body. The first one is when cancer cells move through the blood or lymph system through the body and destroy flourishing tissues, which is known as invasion. The second issue is when cancerous cells try to split and proliferate to nourish themselves in a practice producing new blood vessels, which is known as angiogenesis. Hence a tumor destroys the other healthy tissues and can spread throughout the body in a process called metastasizing. This phase is more difficult to treat, and this process is called metastasis [3]. Because of its significance, fractional calculus is useful to depict dynamic processes in several fields together with physics, economics and finance [ 4 , 5 ], engineering, biology and medicine, and plenty of other current fields [ 6 ]. The inclusion of reminiscence and Symmetry 2023,15, 533. https://doi.org/10.3390/sym15020533 https://www.mdpi.com/journal/symmetry Symmetry 2023,15, 533 2 of 13 genetic markers, which give a extra rational approach to the cancer remedy epidemic version, highlights the requirement of managing fractional order outputs. This was in part assisted by Saeedian et al. who built the epidemic model and learning [ 7 ], which study the behavior and the effect of memory on the disease’s spread. Considering the ABC fractional operators, the power of the smoking model type and its community was given by Ucar et al. [ 8 ]. The famous German chemist Paul Ehelich began to develop therapeutic drugs for infectious diseases in the early 1900s. He began to use the term “chemotherapy”, which he described as the use of chemicals to treat cancer of its powerful anti-inflammatory chemicals. He was the first person to do so. Although Ehrlich was not optimistic about the future, he was particularly interested in cancer treatment, including aniline dyes and early alkylating agents. He wrote “Give up hope you who enter”. In the 1960s surgery and radiotherapy dominated in the field of cancer treatment until the standard of treatment after all local therapies had stabilized when the 33A statistical fractional-order model protective test method for the body became clear by using the framework of fractional differential equations (FDEs) [ 9 ]. Mathematical models are used to analyze the interaction between different tumor cells and antibodies and dengue based on the system of fractional differential equations [ 10 ]. FDEs and a partial mathematical model have been used as an alternative operator to discuss the clinical effects of diabetes and the coexistence of tuberculosis [11]. The fractional derivative is initially divided into two major types. The fractionals with a singular kernel are Riemann–Liouville (RL) and Caputo [ 12 ]. The fractionals without singular kernels are ABC (Mittag–Leffler) and Caputo–Fabrizio (exponential) [ 13 , 14 ]. Fractional calculus is used in finance, chemical, biological, pharmaceutical, physical, and engineering fields as it has many application in our daily life [ 15 , 16 ]. Furthermore, several applications are given in [ 17 – 22 ] for fractional order versions. With the aid of Caputo–Fabrizio, fractional-order equal-width equations were solved using the homotopy perturbation transform approach in [ 23 , 24 ]. At both of the fractional-order epidemic model’s steady states, local and global stability are examined. After analytical treatment, the fractional-order epidemic model is numerically solved using a template that preserves the structure [ 25 , 26 ]. This paper aims to examine the fractional order cancer model with stem cells and chemotherapy. Furthermore, we verified the results using the advanced technique of Atangana–Toufik. Positivity of the proposed model was also derived. In terms of uniqueness of our study, one important point was the stability analysis of a scheme furnished with the aid of fixed point theorem. 2. Basic Concepts Definition 1 (Ref. [ 27 ]) . Atangana–Baleanu in the Liouville–Caputo sense (ABC) derivative is as follows ABC γDγ t{f(t)}=AB(γ) 1−γZt γ d dw f(w)Eγ[−γ(t−w)γ 1−γ]dw,n−1<γ<n, (1) where Eγ represents the function as Mittag–Leffler, and AB(γ) represents the function as normalization and AB(0) = AB(1) = 1. The Laplace transform is given for above as [ABC γDγ tf(t)](s) = AB(γ) 1−γ sγL[f(t)](s)−sγ−1f(0) sγ+γ 1−γ . (2) By applying the Sumudu transform (ST), we obtain the following result: ST[ABC γDγ tf(t)](s) = AB(γ) 1−γ(γΓ(γ+1)Eγ(−1 1−γνγ)) ×[ST(f(t)) −f(0)]. (3) Symmetry 2023,15, 533 3 of 13 3. Materials and Method In order to treat the majority of cancers by using the handiest stem cellular therapy, we developed a fractional order mathematical model that considered three populations: T(t) tumor cells, E(t) effector cells, and S(t) stem cells. Furthermore, we added an amplification A at the realization of the simplified ODEs that described the interaction among all three populations. M(t) is a chemotherapeutic attention medication, and specifics of the parameters and their values are given in [ 28 , 29 ]. The subsequent equations deliver an ABC derivative-based fractional order version for cancer. ABC 0Dγ tS(t) = γ1S−kSMS, ABC 0Dγ tE(t) = α−µE+p1ES (S+1)−p2(T+M)E, ABC 0Dγ tT(t) = r(1−bT)T−(p3E+kTM)T, ABC 0Dγ tM(t) = −γ2M+V(t), (4) with beginning conditions as S0(t) = S(0),E0(t) = E(0),T0(t) = T(0),M0(t) = M(0). (5) By applying Sumudu transform operator on both sides, we obtain OγEγ(−1 1−γωγ)ST[S(t)) −S(0)]=ST[γ1S−kSMS], OγEγ(−1 1−γωγ)ST[E(t)) −E(0)]=ST[α−µE+p1ES (S+1)−p2(T+M)E], OγEγ(−1 1−γωγ)ST[T(t)) −T(0)]=ST[r(1−bT)T−(p3E+kTM)T], OγEγ(−1 1−γωγ)ST[M(t)) −M(0)]=ST[−γ2M+V(t)], (6) where Oγ=B(γ)γΓ(γ+1) 1−γsystem (7) becomes ST[S(t)] = S(0) + 1 OγEγ(−1 1−γωγ)×ST[γ1S−kSMS], ST[E(t)] = E(0) + 1 OγEγ(−1 1−γωγ)×ST[α−µE+p1ES (S+1)−p2(T+M)E], ST[T(t)] = T(0) + 1 OγEγ(−1 1−γωγ)×ST[r(1−bT)T−(p3E+kTM)T], ST[M(t)] = M(0) + 1 OγEγ(−1 1−γωγ)×ST[−γ2M+V(t)]. (7) Using inverse Sumudu Transform, we have S(t) = S(0) + ST−11 OγEγ(−1 1−γωγ)×ST[γ1S−kSMS], E(t) = E(0) + ST−11 OγEγ(−1 1−γωγ)×ST[r(1−bT)T−(p3E+kTM)T], Symmetry 2023,15, 533 4 of 13 T(t) = T(0) + ST−11 OγEγ(−1 1−γωγ)×ST[r(1−bT)T−(p3E+kTM)T], M(t) = M(0) + ST−11 OγEγ(−1 1−γωγ)×ST[−γ2M+V(t)]. (8) Therefore, we obtain S(j+1)(t) = Sj(0) + ST−1{1 OγEγ(−1 1−γωγ)×ST[γ1Sj(t)−kSj(t)Mj(t)Sj(t)]}, E(j+1)(t) = Ej(0) + ST−1{1 OγEγ(−1 1−γωγ)×ST[α−µEj(t) + p1Em(t)Sj(t) (Sj(t) + 1) −p2(Tj(t) + Mj(t))Ej(t)]}, T(j+1)(t) = Tj(0) + ST−1{1 OγEγ(−1 1−γωγ)×ST[r(1−bTj(t))Tj(t) −(p3Ej(t) + kTj(t)Mj(t))Tj(t)]}, M(j+1)(t) = Mj(0) + ST−1{1 OγEγ(−1 1−γωγ)×ST[−γ2Mj(t) + V(t)]}. (9) The solution of system (4) is represented as S=lim j→∞ Sj;E=lim j→∞ Ej;T=lim j→∞ Tj(t);M=lim j→∞ Mj. (10) Positivity of Solutions with ABC Operator All solutions are positive if all of the beginning conditions are true for nonlocal operators. We need to define the norm kΠk∞=Supt∈DΠ|Π(t)|, (11) such that DΠ is the domain of Π . By applying this norm, we obtain for the Atangana– Baleanu derivative S(t)≥S0Eγ−γkSkMk∞−γ1 AB(γ)−(1−γ)(kSkMk∞−γ1)t,∀t>0 (12) E(t)≥E0Eγ−γµ −p1kSk∞ kSk∞+1+p2(kTk∞+kMk∞ AB(γ)−(1−γ)(µ−p1kSk∞ kSk∞+1+p2(kTk∞+kMk∞)t,∀t>0 (13) T(t)≥T0Eγ−γp3kEk∞+kTkMk∞ AB(γ)−(1−γ)(p3kEk∞+kTkMk∞)t,∀t>0 (14) M(t)≥M0Eγ−γγ2 AB(γ)−(1−γ)(γ2)t,∀t>0 (15) Theorem 1. Assume (X , |·|) to be a Banach space and consider H to be a self-map of X satisfying kHr1−Hxk ≤ θkX−Hr1k+θkr1−xk, (16) for every r1,xeX, where 0≤θ<1. Suppose that the system (4), and we obtain the following result as 1−γ B(γ)γΓ(γ+1)Eγ(−1 1−γωγ). (17) Symmetry 2023,15, 533 5 of 13 Proof. Defining Kas a self-map, we may then write this as K[S(j+1)] = S(j+1)=Sj(0) + ST−1[1 OγEγ(−1 1−γωγ)×ST[γ1Sj(t)−kSj(t)Mj(t)Sj(t)], K[E(j+1)] = E(j+1)=Ej(0) + ST−1[1 OγEγ(−1 1−γωγ) ×ST[α−µEj(t) + p1Ej(t)Sj(t) (Sj(t) + 1)−p2(Tj(t) + Mj(t))Ej(t)], K[T(j+1)(t)] = T(j+1)(t) = Tj(0) + ST−1[1 OγEγ(−1 1−γωγ) ×ST[r(1−bTj(t))Tj(t)−(p3Ej(t) + kTj(t)Mj(t))Tj(t)], K[M(j+1)(t)] = M(j+1)(t) = Mj(0) + ST−1[1 OγEγ(−1 1−γωγ) ×ST[−γ2Mj(t) + V(t)]. (18) Using the norm’s aspects along with triangular inequality, kK[Sj(t)] −K[Si(t)]k ≤ kSj(t)−Si(t)k+kST−1{1−γ BEγ(−1 1−γωγ) ×ST[γ1Sj(t)−KSj(t)Mj(t)Sj(t)]} − ST−1{1−γ BEγ(−1 1−γωγ) ×ST[γ1Si(t)−KSi(t)Mi(t)Si(t)]}k, kK[Ej(t)] −K[Ei(t)]k ≤ kEj(t)−Ei(t)k+kST−1[1−γ BEγ(−1 1−γωγ) ×ST[α−µEj(t) + p1Ej(t)Sj(t) (Sj(t) + 1)−p2(Tj(t) + Mj)(t)Ej(t)] −ST−1[1−γ BEγ(−1 1−γωγ)×ST[α−µEi(t) + p1Ei(t)Si(t) (Si(t) + 1)−p2(Ti(t) + Mi(t))Ei(t)]}k, kK[Tj(t)] −K[Ti(t)]k ≤ kTj(t)−Ti(t)k+kST−1[1−γ BEγ(−1 1−γωγ) ×ST[r(1−BTj(t))Tj(t)−(p3Ej(t) + KTj(t)Mj(t))Tj(t)]} − ST−1[1−γ BEγ(−1 1−γωγ) ×ST[r(1−BTi(t))Ti(t)−(p3Ei(t) + KTi(t)Mi(t))Ti(T)]}k, kK[Mj(t)] −K[Mi(t)]k ≤ kMj(t)−Mi(t)k+kST−1[1−γ BEγ(−1 1−γωγ) ×ST[−γ2Mj(t) + V(t)]} − ST−1[1−γ BEγ(−1 1−γωγ) ×ST[−γ2Mi(t) + V(t)]}k, (19) where B=B(γ)γΓ(γ+1)K satisfied when Symmetry 2023,15, 533 6 of 13 θ= (0, 0, 0, 0) =                                            kSj(t)−Si(t)k × k − Sj(t) + Si(t)k +γ1kSj(t)−Si(t)k − kkSj(t)−Si(t)kkMj(t)−Mi(t)kkSj(t)−Si(t)k, ×kEj(t)−Ei(t)k × k − Ej(t) + Ei(t)k +α−µkEj(t)−Ei(t)k+p1kEj(t)−Ei(t)kkSj(t)−Si(t)k (kSj(t)−Si(t)k+1), −p2(kTj(t)−Ti(t)k+kMj(t)−Mi(t)k)kEj(t)−Ei(t)k ×kTj(t)−Ti(t)k × k − Tj(t) + Ti(t)k +r(1−bkTj(t)−Ti(t)k)kTj(t)−Ti(t)k − (p3k − Ej(t) + Ei(t)k +kkTj(t)−Ti(t)kkMj(t)−Mi(t)k)kTj(t)−Ti(t)k, ×kMj(t)−Mi(t)k × k − Mj(t) + Mi(t)k −γ2kMj(t)−Mi(t)k+V(t), (20) and we find that Kis Picard K-stable. Theorem 2. System (9) is a unique and distinct solution found by utilizing the iteration method Proof. Let us consider the Hilbert space, H=L2((q,p)×(0, T)) h:(q,p)×[0, T]→R,Z Z ghdgdh <∞. (21) For this purpose, the following operators are used θ= (0, 0, 0, 0) =            γ1S−kSMS, α−µE+p1ES (S+1)−p2(T+M)E, r(1−bT)T−(p3E+kTM)T, −γ2M+V(t). (22) We have T(S11(t)−S12(t),E21(t)−E22(t),T31(t)−T32(t),M41(t)−M42(t),(v1,v2,v3,v4), (23) where ( S 11(t)− S 12(t) ,E 21(t)− E 22(t) ,T 31(t)− T 32(t) ,M 41(t)− M 42(t) , which displays the system’ Sspecial solutions. We may obtain this by using the inner function and norm. {γ1A−kADA,v1} ≤ γ1kAkkv1k − kkAkkv1kkDkkv1kkAkkv1k, {α−µB+p1A (A+1)−p2(C+D)B,v2} ≤ α−µkBkkv2k +p1kBkkAkkv2k (kAkkv2k+1)−p2(kCkkv2k+kDkkv2k)kBkkv2k, {r(1−bC)C−(p3(B+kCkD)C,v3} ≤ r(1−bkCkkv3k)Ckv2k − (p3kBkkv3k +kkBkkv3kkDkkv2k)kCkkv3k, {−γ2D+V(t),v4} ≤ −γ2kDkkv4k+V(t)kv4k, where A= S 11(t)− S 12(t) , B= E 21(t)− E 22(t) , C= T 31(t)− T 22(t) and D= M 41(t)− M 42(t) . In the case of a large number E 1 ,E 2 ,E 3 , and E 4 , all of the results converge to an exact solution. By utilizing four positive and very small parameters and the topology concept, we have χE1,χE2,χE3,χE4. Symmetry 2023,15, 533 7 of 13 kS(t)−S11(t)k,kS(t)−S12(t)k<χE1 ω, kE(t)−E21(t)k,kE(t)−E22(t)k<χE2 ς, kT(t)−T31(t)k,kT(t)−T32(t)k<χE3 ϑ, and kM(t)−M41(t)k,kM(t)−M32(t)k<χE4 $, where ω=4(γ1kAk − kkAkkv1kkDkkAk)kv1k, ς=4(α−µkBk+p1kBkkAk (kAkkv2k+1)−p2(kCk+kDk)kBk)kv2k, ϑ=4(r(1−bkCkkv3k)kCk − (p3kBk+kkBkkDk)kCk)kv3k, $=4(−γ2kDk+V(t))kv4k, where γ1kAk − kkAkkv1kkDkkAk 6=0, α−µkBk+p1kBkkAk (kAkkv2k+1)−p2(kCk+kDk)kBk 6=0, r(1−bkCkkv3k)kCk − (p3kBk+kkBkkCk)kDk 6=0, −γ2kDkk +V(t)) 6=0, where kv1k,kv2k,kv3k,kv4k 6=0; kS11(t)−S12(t)k, kE21(t)−E22(t)k,kT31(t)−T32(t)k,kM41(t)−M42(t)k=0. S11(t) = S12(t);E21(t) = E22(t);T31(t) = T32(t);M41(t) = M42(t). (24) The uniqueness proof is now complete. 4. Numerical Scheme with Atangana–Toufik In this article, an advanced scheme was applied for nonlinear FD equations on account of FD with nonsingular kernel and non-nearby fractional derivative. For this purpose, recollect the nonlinear system given in (5) and apply the technique we have used. S(t)−S(0) = (1−γ) ABC(γ){γ1S(t)−kS(t)M(t)S(t)} +γ Γ(γ)×ABC(γ)Zt 0{γ1S(τ1)−kS(τ1)M(τ1)S(τ1)}(t−τ1)γ−1dτ1, E(t)−E(0) = (1−γ) ABC(γ){α−µE(t) + p1E(t)S(t) (S(t) + 1)−p2(T(t) + M(t))E(t)} +γ Γ(γ)×ABC(γ)Zt 0{α−µE(τ1) + p1E(τ1)S(τ1) (S(τ1) + 1)−p2(T(τ1) + M(τ1))E(τ1)} (t−τ1)γ−1dτ1, T(t)−T(0) = (1−γ) ABC(γ){r(1−bT(t))T(t)−(p3E(t) + kT(t)M(t))T(t)} Symmetry 2023,15, 533 8 of 13 +α1 Γ(γ)×ABC(γ)Zt 0{r(1−bT(τ1))T(τ1)−(p3E(τ1) + kT(τ1)M(τ1))T(τ1)} (t−τ1)γ−1dτ1, M(t)−M(0) = (1−γ) ABC(γ){−γ2M(t) + V(t)} +γ Γ(γ)×ABC(γ)Zt 0{−γ2M(τ1) + V(t)}(t−τ1)γ−1dτ1. As given tM+1,M=0, 1, 2, 3 . . ., then the above equation can be reformulated as S(tM+1)−S(0) = (1−γ) ABC(γ){γ1S(tM)−kS(tM)M(tM)S(tM)} +γ Γ(γ)×ABC(γ) M ∑ k1=0Ztk1+1 tk1 {γ1S(τ1)−kS(τ1)M(τ1)S(τ1)} (tM+1−τ1)γ−1dτ1, E(tM+1)−E(0) = (1−γ) ABC(γ){α−µE(tM) + p1E(tM)S(tM) (S(tM) + 1)−p2(T(tM) + M(tM))E(tM)} +γ Γ(γ)×ABC(γ) M ∑ k1=0Ztk1+1 tk1 {α−µE(τ1) + p1E(τ1)S(τ1) (S(τ1) + 1)−p2(T(τ1) + M(τ1))E(τ1)} (tM+1−τ1)γ−1dτ1, T(tM+1)−T(0) = (1−γ) ABC(γ){r(1−bT(tM))T(tM)−(p3E(tM) + kT(tM)M(tM))T(tM)} +γ Γ(γ)×ABC(γ) M ∑ k1=0Ztk1+1 tk1 {r(1−bT(τ1))T(τ1)−(p3E(τ1) + kT(τ1)M(τ1))T(τ1)} (tM+1−τ1)γ−1dτ1, M(tM+1)−M(0) = (1−γ) ABC(γ){−γ2M(tM) + V(t)} +γ Γ(γ)×ABC(γ) M ∑ k1=0Ztk1+1 tk1 {−γ2M(τ1) + V(t)}(tM+1−τ1)γ−1dτ1. By using the above equation, we have SM+1=S0+(1−γ) ABC(γ){γ1S(tM)−kS(tM)M(tM)S(tM)} +γ Γ(γ)×ABC(γ) M ∑ k1=0 (γ1Sk1−kSk1 Mk1Sk1 hB1 −γ1Sk1−1−kSk1−1Mk1−1Sk1−1 hAγ,k1,2), EM+1=E0+(1−γ) ABC(γ){α−µE(tM) + p1E(tM)S(tM) (S(tM) + 1)−p2(T(tM) + M(tM))E(tM)} +γ Γ(γ)×ABC(γ) M ∑ k1=0 ( α−µEk1+p1Ek1Sk1 (Sk1+1)−p2(Tk1+Mk1)Ek1 hB1 Symmetry 2023,15, 533 9 of 13 − α−µEk1−1+p1Ek1−1Sk1−1 (Sk1−1+1)−p2(Tk1−1+Mk1−1)Ek1−1 hAγ,k1,2), TM+1=T0+(1−γ) ABC(γ){r(1−bT(tM))T(tM)−(p3E(tM) + kT(tM)M(tM))T(tM)} +γ Γ(γ)×ABC(γ) M ∑ k1=0 (r(1−bTk1)Tk1−(p3Ek1+kTk1 Mk1)Tk1 hB1 −r(1−bTk1−1)Tk1−1−(p3Ek1−1+kTk1−1Mk1−1)Tk1−1 hAγ,k1,2), MM+1=M0+(1−γ) ABC(γ){−γ2M(tM) + V(t)} +γ Γ(γ)×ABC(γ) M ∑ k1=0 (−γ2Mk1+V(t) hB1 −−γ2Mk1−1+V(t) hAγ,k1,2), where Aγ,k1,2 =Rtk1+1 tk1(τ1−tk1)(tM+1−τ1)γ−1dτ1 and B1=Rtk1+1 tk1(τ1−tk1−1)(tM+1− τ1)γ−1dτ1. By integrating the above and putting in a system of equations, we have SM+1=S0+(1−γ) ABC(γ){γ1S(tM)−kS(tM)M(tM)S(tM)} +γ ABC(γ) M ∑ k1=0 (hγ{γ1Sk1−kSk1 Mk1Sk1} hB1 −hγ{γ1Sk1−1−kSk1−1Mk1−1Sk1−1} hAγ,k1,2), EM+1=E0+(1−γ) ABC(γ){α−µE(tM) + p1E(tM)S(tM) (S(tM) + 1)−p2(T(tM) + M(tM))E(tM)} +γ ABC(γ) M ∑ k1=0 ( hγ{α−µEk1+p1Ek1Sk1 (Sk1+1)−p2(Tk1+Mk1)Ek1} hB1 − hγ{α−µEk1−1+p1Ek1−1Sk1−1 (Sk1−1+1)−p2(Tk1−1+Mk1−1)Ek1−1} hAγ,k1,2), TM+1=T0+(1−γ) ABC(γ){r(1−bT(tM))T(tM)−(p3E(tM) + kT(tM)M(tM))T(tM)} +γ ABC(γ) M ∑ k1=0 (hγ{r(1−bTk1)Tk1−(p3Ek1+kTk1 Mk1)Tk1} hB1 −hγ{r(1−bTk1−1)Tk1−1−(p3Ek1−1+kTk1−1Mk1−1)Tk1−1} hAγ,k1,2), MM+1=M0+(1−γ) ABC(γ){−γ2M(tM) + V(t)} +γ ABC(γ) M ∑ k1=0 (hγ{−γ2Mk1+V(t)} hB1