# On a Mathematical Model of Tumor-Immune Interaction with a Piecewise Differential and Integral Operator

Authors: Shahram Rezapour; Chernet Tuge Deressa; Robert G. Mukharlyamov; Sina Etemad
Journal: Journal of Mathematics
 (2022)

---

## Abstract

The representation of mathematical models via piecewise differential and integral operators for dynamic systems has this potential to capture cross-over behaviors such as a passage from deterministic to randomness which can be exhibited by different systems. A 3D mathematical model, similar to the prey-predator system, of tumor-immune interaction with piecewise differential and integral operators is developed and analyzed. Three different scenarios, namely, cross-overs from deterministic to randomness, the Mittag-Leffler law to randomness, and a cross-over behavior from fading memory to the power-law and a random process, are considered. The existence, uniqueness, positivity, and boundedness of the solutions of the systems are proved via the linear growth and Lipschitz conditions. The numerical approximations by Toufik, Atangana, and Araz are used for approximation of solutions and simulation of the piecewise models in different scenarios. From the nondimensionalized version of the 3D model representation, it is shown that the parameter values have an impact on the growth of tumor cells, and activating the proliferation of the resting cells has negatively affected the development of tumor cells. Moreover, the dynamics of tumor-immune interaction exhibited a cross-over behavior, and this behavior is exposed by the piecewise modeling approach used for the representations.


---

## Body

## 1. Introduction

The classical differential equations (ordinary and partial differential equations) were developed based on the concept of rate of changes, and it has been used and investigated for several decades. They have been used in developing several mathematical models representing real-world problems and are effectively used to make their analysis. Nevertheless, some drawbacks with the classical differential equations were observed. The classical differential equations are not efficient in replicating observed realities. For instance, some cases require randomness and cannot be captured by the classical differential equations, and as a result, the stochastic differential equations came into being and have been used successfully.In the same way, there are some problems in real world that cannot be captured by stochastic differential equations, and this led to the development of different concepts of fractional derivatives and integrals. Different concepts of fractional operators have been used to capture trends including nonlocalities, power-law processes, memory effects, fractal processes, and some other real-world problems. Though there are different endeavors made by mathematicians to capture different real-world problems using mathematical models, the issue of capturing dynamic systems that display multiple behaviors is not fully addressed (see [1] and the references therein). With this understanding and to solve the problem of capturing dynamics of real-world problems with cross-over behaviors, Atangana and Araz [1, 2] developed a novel concept named piecewise modeling that involves differential and integral operators.In this study, the notion of piecewise modeling is considered to develop different piecewise mathematical models for tumor-immune interaction using the classical, stochastic, and fractional derivative concepts.There are several studies conducted on tumor-immune interaction by using different mathematical models based on classical and fractional derivative concepts. A few of them are reviewed as follows: Kaur and Ahmad [3] developed a mathematical model of tumor-immune interaction by including the Michaelis-Menten function in the model. The authors used the classical derivative and showed that the inclusion of the Michaelis-Menten function helped in achieving stability of the dynamical system and increased the rate of growth of resting cells. A seven-dimensional dynamical model of the tumor-immune system equipped with the Reimann-Liouville fractal-fractional operator with the Mittag-Leffler-type kernel was considered by Farman et al. [4]. The results showed that the IL-12 cytokine and anti-PD-L1 inhibitor increased the immune system and decreased the cancer cells. Ahmed et al. [5] applied ABC fractal-fractional operators to develop a mathematical model and visualize the tumor-immune relationship. Various mathematical models addressing different cancer treatments such as cytotoxic chemotherapy, immunotherapy, and their combination are investigated by Depillis et al. [6].Wilkie [7] discussed different mathematical models with the classical derivatives of tumor growth in the presence of an immune response, and the findings suggest that feedback from the tumor to the immune response induces the existence of dormant cancer cells. A mathematical model incorporating three types of immunotherapy and focusing on the inhibitory role of Tregs in the tumor-immune system is developed and analyzed by Zhongtao et al. [8]. A mathematical model describing how cancer cell progresses and survives an encounter with the immune cell population is developed and discussed in [9]. The chaotic dynamics of a tumor-immune interaction model with delay are considered in [10]. Kuznetsov et al. [11] developed a tumor-immune interaction mathematical model and described the response of effector cells to the growth of an immunogenic tumor. There are many other studies made on the tumor-immune interaction (see for instance [12–16]).As the concept of piecewise differential and integral operators is relatively new, hence, little literature is available. The concept of the piecewise derivative and integral operator is used with different fractional derivatives to investigate an SIR mathematical model of COVID-19 in [17]. Zeb et al. [18] investigated a five-dimensional compartmental model of COVID-19 with the concept of piecewise derivative and integral operators combining the Caputo-Fabrizio, classical, and stochastic differential equations. A mathematical model representing an interaction in a food web is considered, and the concept of piecewise differential and integral operator is imposed for investigation in [19]. A mathematical model of the third wave of COVID-19 is developed and considered with the concept of piecewise differential and integral operators [2]. The concept of piecewise differentiation and integral operator is applied to a CAR-T cells-SARS-2 virus model [20]. Of course, different concepts of fractional derivatives have been used by several authors for investigation of different dynamical systems and different applications (see for instance [21–28] and the references therein).To the best of the researchers’ knowledge, there is no single study conducted on a Tumor-immune interaction mathematical model in the sense of piecewise differential and integral operators. The concept of piecewise derivatives and integrals in capturing real-world problems with multiple behaviors is a novel result as it empowers researchers in the area to use different concepts of derivative and integral operators at the same time to study multiple behaviors of a given dynamic system which may not otherwise be possible.This study, therefore, focuses on discovering different cross-over behaviors in a mathematical model of tumor-immune interaction in the sense of the piecewise differential and integral operators developed by Atangana and Seda. The classical differential equations, stochastic differential equations, and different concepts of fractional operators are included in the formation of the piecewise differential and integral operators. Accordingly, three different scenarios of the mathematical models were developed: a cross-over from deterministic to randomness, a cross-over from Mittag-Leffler law to randomness, and a cross-over from exponential decay to power-law and random process.The remaining part of this paper is organized as follows: The formulation and description of the model, the parameters and their description, and the formulation of three piecewise models representing the system are considered in Section2. Section 3 is devoted to the existence, uniqueness, positivity, and boundedness of the piecewise models. Numerical approximations of the piecewise models are considered in Section 4 followed by simulations in Section 5. The conclusion is provided in Section 6 followed by the list of references in the last section.


## 2. Formulation of Models

In this section, the tumor-immune mathematical model used in this study is described. The system comprises three nonlinear differential equations that modify to different concepts of fractional operators and stochastic differential equations. The model involves the concentration of tumor cells at timet, represented by Xt, the concentration of hunting predator cells at timet, represented by Yt, and the concentration of resting predator cells at timet, represented by Zt. The hunting and resting predator cells are normal tissue cells. The model is similar to the prey-predator system and originally developed by Kaur and Ahmad [3] describing the growth, death, and interaction among this population and is given as shown in:
(1)Ẋ=Λ+r1X1−Xk1−α1XY,Ẏ=βYZ−d1Y−α2YX,Ż=r2Z1−Zk2−βYZ−d2Z+ρXZX+η,where all the parameters in (1) are nonnegative and the initial conditions are
(2)X0>0,Y0>0,Z0>0,andk1>k2.The first equation of (1) describes the rate of growth of concentration of tumor cells. It is assumed that tumor cells follow the logistic growth in the absence of any immune intervention (hunting and resting predator CTL cells). The extinction of hunting cells and tumor cells is proportional to the densities of both the cells, in line with the law of mass action. The proliferation of the resting cells is also assumed to follow the logistic growth function in the absence of tumor cells. The multiplying of the resting cells is enhanced by the tumor cells characterized by the term ρXZ/X+η called the Michael-Menten function which indicates the saturation effect of the resting predator cells, with a rate of proliferation ρ and a half-saturation constant η. The resting cells are converted to hunting cells by direct contact with them or by fast diffusing substance (cytokines) produced by hunting cells at the rate β. It is worth mentioning that inactivated hunting cells will not get back to the resting stage once over. The parameters of the model and their description used in this study are summarized in Table 1.Table 1
Parameters of the model and their descriptions.
ParameterDescriptionDimensional valuesSourceΛThe rate of conversion of normal tissue cells to malignant (cancerous) cells (fixed input)1000000/cell/dayAssumedr1The growth rate of the tumor cellsX.0.18/day[14]r2The growth rate of the resting cellsZ0.1045/day[29]βThe rate of conversion of resting cellZ to hunting cell Y4.32×10−8/cell/day[14]k1Carrying capacity of tumor cellsX5×106/cellk2Carrying capacity of resting cellsZ3×106 cellρThe proliferation rate of resting cellZ0.49545/dayηThe value at which the growth rate of resting immune cellsZ gets half of its maximum value.1000000d1Apoptotic or natural death rate of the hunting cellsY.0.0412/day[29]d2Apoptotic or natural death rate of the resting cellsZ.0.0412/dayα1The rate of inactivation of tumor cellT by hunter cells Y.1.101×10−7/cell/dayα2The rate of inactivation of hunting cellY by tumor cell X2.2×10−8 /cell/dayFollowing the work in [1], let us use the following dimensionless variables in the system:
(3)t∗=Λt/k1,X∗=X/k1,Y∗=αk1Y/Λ,Z∗=Z/k2.After applying the dimensionless variables in (2) to the systems (1), we obtain
(4)Ẋ=1+c1X1−X−XY,Ẏ=c2YZ−c3Y−c4YX,Ż=c5Z1−Z−c6YZ−c7Z+c8XZX+K,where
(5)c1=r1k1Λ,c2=βk2k1Λ,c3=d1k1Λ,c4=α2k12Λc5=r2k1Λ,c6=βα1,c7=d2k1Λ,c8=ρk1Λ,K=ηk1.In this study, the nondimensionalized mathematical model (4) is considered in the sense of piecewise differential and integral operators using classical, stochastic differential equations and different concepts of fractional and integral operators.


### 2.1. Preliminaries

The basic definitions of different fractional derivatives and integral operators used in the study are recalled as follows.Definition 1 (see [30–32]).
Letμ∈0,1, and f∈C10,t. The fractional ABC (with Mittag-Leffler kernel), Caputo (with the power-law kernel), and Caputo-Fabrizio (with exponential decay kernel) derivatives are, respectively, defined as follows:
(6)D0ABCtμft=Gμ1−μ∫0tddτfτEμ−μ1−μt−τμdτ,Dtμft=1Γ1−μ∫0tddτ0cfτt−τ−μdτ,D0CFtμft=Gμ1−μ∫0tddτfτexp−μ1−μt−τdτ,where Gμ=1−μ+μ/Γμ is the normal operator,Eμ. is the Mittag-Leffler function, and Γ. is the Euler Gamma function.
The fractional integrals of the Caputo, Caputo-Fabrizio, and ABC types are, respectively, given by(7)I0Ctμft=1Γμ∫0tt−ρμ−1fρdρ,I0CFtμft=1−μGμft+μGμ∫0tfρdρ,Itμ0ABCft=1−μGμft+μGμΓμ∫0tfρt−ρμ−1dρ.


### 2.2. Equilibrium Points

Four equilibrium points of the system (4) are given below:
(8)E1=121+1+4c1,0,0,E2=X2,0,Z2,whereX2=121+1+1c1,Z2=1c5c5−c7+c8X2X2+K,E3=X3,Y3,0,whereX3=−c3c4<0,(9)E4=X4,Y4,Z4istheinteriorequilibriumpoint,where
(10)C3X43+C2X42+C1X4+C0=0,and C1=c5K−c7K+c5c3K/c2+c6c1K+1,
(11)C2=c5−c7+c8−c5c2c3+c4K−c1c61−K,(12)C3=c1c6−c4c5c2,C0=−c6K<0.


### 2.3. Formulation of Piecewise Model

In this subsection, we formulate three different scenarios of piecewise differential operator representations for the tumor-immune interaction model given in (4).Scenario 1.
In this case, we consider a piecewise model that involves a passage from the deterministic to a random process:
For0≤t≤t1,X0=X11,Y0=Y12,Z0=Z13,(13)Ẋ=1+c1X1−X−XY,Ẏ=c2YZ−c3Y−c4YX,Ż=c5Z1−Z−c6YZ−c7Z+c8XZX+K,for t1≤t≤T,Xt1=X21,Yt1=Y22,Zt1=Z23,(14)dX=1+c1X1−X−XYdt+σ1XdB1t,dY=c2YZ−c3Y−c4YXdt+σ2YdB2t,dZ=c5Z1−Z−c6YZ−c7Z+c8XZX+Kdt+σ3ZdB3t.The mathematical model (10) has a deterministic character, and it is extended to the stochastic model described in (11) by adding a white noise-type perturbation to the system. The parameters σ1,σ2,σ3 are positive constants and are the intensities of the random disturbances. Bt=B1t,B2t,B3t is the white noise process.Scenario 2.
In this case, we have considered a piecewise model involving a passage from the Mittag-Leffler law to a random process:
For0≤t≤t1,X0=X11,Y0=Y12,Z0=Z13,(15)D0ABCtαX=1+c1X1−X−XY,D0ABCtαY=c2YZ−c3Y−c4YX,D0ABCtαZ=c5Z1−Z−c6YZ−c7Z+c8XZX+K,for t1≤t≤t2,Xt1=X21,Yt1=Y22,Zt1=Z23,(16)dX=1+c1X1−X−XYdt+σ1X−X4dB1t,dY=c2YZ−c3Y−c4YXdt+σ2Y−Y4dB2t,dZ=c5Z1−Z−c6YZ−c7Z+c8XZX+Kdt+σ3Z−Z4dB3t.Scenario 3.
In this case, we consider a piecewise model from fading memory to the power-law and then to a random process:
For0≤t≤t1,X0=X11,Y0=Y12,Z0=Z13(17)D0CFtαX=1+c1X1−X−XY,D0CFtαY=c2YZ−c3Y−c4YX,D0CFtαZ=c5Z1−Z−c6YZ−c7Z+c8XZX+K,for t1≤t≤t2,Xt1=X21,Yt1=Y22,Zt1=Z23,(18)D0CtαX=1+c1X1−X−XY,D0CtαY=c2YZ−c3Y−c4YX,D0CtαZ=c5Z1−Z−c6YZ−c7Z+c8XZX+K,for t2≤t≤T,Xt2=X31,Yt2=Y32,Zt2=Z33,(19)dX=1+c1X1−X−XYdt+σ1XdB1t,dY=c2YZ−c3Y−c4YXdt+σ2YdB2t,dZ=c5Z1−Z−c6YZ−c7Z+c8XZX+Kdt+σ3ZdB3t.


## 2.1. Preliminaries

The basic definitions of different fractional derivatives and integral operators used in the study are recalled as follows.Definition 1 (see [30–32]).
Letμ∈0,1, and f∈C10,t. The fractional ABC (with Mittag-Leffler kernel), Caputo (with the power-law kernel), and Caputo-Fabrizio (with exponential decay kernel) derivatives are, respectively, defined as follows:
(6)D0ABCtμft=Gμ1−μ∫0tddτfτEμ−μ1−μt−τμdτ,Dtμft=1Γ1−μ∫0tddτ0cfτt−τ−μdτ,D0CFtμft=Gμ1−μ∫0tddτfτexp−μ1−μt−τdτ,where Gμ=1−μ+μ/Γμ is the normal operator,Eμ. is the Mittag-Leffler function, and Γ. is the Euler Gamma function.
The fractional integrals of the Caputo, Caputo-Fabrizio, and ABC types are, respectively, given by(7)I0Ctμft=1Γμ∫0tt−ρμ−1fρdρ,I0CFtμft=1−μGμft+μGμ∫0tfρdρ,Itμ0ABCft=1−μGμft+μGμΓμ∫0tfρt−ρμ−1dρ.


## 2.2. Equilibrium Points

Four equilibrium points of the system (4) are given below:
(8)E1=121+1+4c1,0,0,E2=X2,0,Z2,whereX2=121+1+1c1,Z2=1c5c5−c7+c8X2X2+K,E3=X3,Y3,0,whereX3=−c3c4<0,(9)E4=X4,Y4,Z4istheinteriorequilibriumpoint,where
(10)C3X43+C2X42+C1X4+C0=0,and C1=c5K−c7K+c5c3K/c2+c6c1K+1,
(11)C2=c5−c7+c8−c5c2c3+c4K−c1c61−K,(12)C3=c1c6−c4c5c2,C0=−c6K<0.


## 2.3. Formulation of Piecewise Model

In this subsection, we formulate three different scenarios of piecewise differential operator representations for the tumor-immune interaction model given in (4).Scenario 1.
In this case, we consider a piecewise model that involves a passage from the deterministic to a random process:
For0≤t≤t1,X0=X11,Y0=Y12,Z0=Z13,(13)Ẋ=1+c1X1−X−XY,Ẏ=c2YZ−c3Y−c4YX,Ż=c5Z1−Z−c6YZ−c7Z+c8XZX+K,for t1≤t≤T,Xt1=X21,Yt1=Y22,Zt1=Z23,(14)dX=1+c1X1−X−XYdt+σ1XdB1t,dY=c2YZ−c3Y−c4YXdt+σ2YdB2t,dZ=c5Z1−Z−c6YZ−c7Z+c8XZX+Kdt+σ3ZdB3t.The mathematical model (10) has a deterministic character, and it is extended to the stochastic model described in (11) by adding a white noise-type perturbation to the system. The parameters σ1,σ2,σ3 are positive constants and are the intensities of the random disturbances. Bt=B1t,B2t,B3t is the white noise process.Scenario 2.
In this case, we have considered a piecewise model involving a passage from the Mittag-Leffler law to a random process:
For0≤t≤t1,X0=X11,Y0=Y12,Z0=Z13,(15)D0ABCtαX=1+c1X1−X−XY,D0ABCtαY=c2YZ−c3Y−c4YX,D0ABCtαZ=c5Z1−Z−c6YZ−c7Z+c8XZX+K,for t1≤t≤t2,Xt1=X21,Yt1=Y22,Zt1=Z23,(16)dX=1+c1X1−X−XYdt+σ1X−X4dB1t,dY=c2YZ−c3Y−c4YXdt+σ2Y−Y4dB2t,dZ=c5Z1−Z−c6YZ−c7Z+c8XZX+Kdt+σ3Z−Z4dB3t.Scenario 3.
In this case, we consider a piecewise model from fading memory to the power-law and then to a random process:
For0≤t≤t1,X0=X11,Y0=Y12,Z0=Z13(17)D0CFtαX=1+c1X1−X−XY,D0CFtαY=c2YZ−c3Y−c4YX,D0CFtαZ=c5Z1−Z−c6YZ−c7Z+c8XZX+K,for t1≤t≤t2,Xt1=X21,Yt1=Y22,Zt1=Z23,(18)D0CtαX=1+c1X1−X−XY,D0CtαY=c2YZ−c3Y−c4YX,D0CtαZ=c5Z1−Z−c6YZ−c7Z+c8XZX+K,for t2≤t≤T,Xt2=X31,Yt2=Y32,Zt2=Z33,(19)dX=1+c1X1−X−XYdt+σ1XdB1t,dY=c2YZ−c3Y−c4YXdt+σ2YdB2t,dZ=c5Z1−Z−c6YZ−c7Z+c8XZX+Kdt+σ3ZdB3t.


## 3. Existence, Uniqueness, Positivity, and Boundedness of the Solutions

This section proves the existence and uniqueness of a solution to the system in Scenario3.Theorem 2.
LetΩ=X,Y,Z∈ℝ3:maxX,Y,Z≤L. For each initial condition U0=X11,Y12,Z13∈Ω, there exists a unique solution of (17) for all t≥0.Proof.
Let0<t1<∞. We want to find a sufficient condition for the existence and uniqueness of the solution of (17) in the domain Ω×0,t1.
Suppose thatG is a mapping such that GW=H1W,H2W,H3W, where W=X,Y,ZT,W′=X/,Y/,Z/T, and
(20)H1X,Y,Z=1+c1X1−X−XY,H2X,Y,Z=c2YZ−c3Y−c4YX,H3X,Y,Z=c5Z1−Z−c6YZ−c7Z+c8XZX+K.
Now for anyW,W′∈Ω, we have
(21)GW−GW′≤XY−X′Y′+c2YZ−Y′Z′+c3Y′−Y+c4Y′X′−YX+c6Y′Z′−YZ+c7Z′−Z+c8Z−Z′+c8XZX+K−c8X′Z′X′+K≤X′Y′−XY+c2YZ−Y′Z′+c3Y′−Y+c4Y′X′−YX+c6Y′Z′−YZ+c7Z′−Z+c8Z−Z′+c8XX+KZ−Z′+Kc8Z′X′+KX+KX−X′≤L+KLc8X′−X+c2L+c3+c4LY−Y′+c6L+c7+2c8Z−Z′≤YW−W′,≤YW−W′,where Y=maxL+KLc8,c2L+c3+c4L,c6L+c7+2c8.
Thus, the mappingG satisfies the Lipschitz condition with respect to W. This proves that the system (17) has a unique solution.Theorem 3.
The solution of (17) is invariant in the set ℝ+3=W∈ℝ3:W≥0andWt=X,Y,ZT..Proof.
By referring to the system (17), we can observe that D0CFtαXtX=0≥0,D0CFtαYtY=0≥0,D0CFtαZtZ=0≥0, and then the mentioned system is nonreducing which proves the invariance of the system in ℝ+3 and the feasible region is Ω=X,Y,Z∈ℝ3:maxX,Y,Z≤L>0.
Thus, in Theorems2 and 3, we find that the system (17) has a unique solution for each initial condition in the feasible set Ω. The existence and uniqueness of the solution of the system (18) can be shown similarly for the domain Ω×t1,t2.
We shall now show the existence and uniqueness of the solution for the stochastic differential equation given in (19) based on Theorem 3. Let us first write the system (19) in a Volterra type of integrals for all t in t2,T.
For simplicity, let us write (19) in the form
(22)dX=H1X,Y,Zdt+G1t,XdB1t,dY=H2X,Y,Zdt+G2t,YdB2t,dZ=H3X,Y,Zdt+G3t,ZdB3t,where t2≤t≤T, and the initial condition is given by Xt2=X31,Yt2=Y32,Zt2=Z33, Hi,i=1,2,3 are defined in (20) and G3t,Z=σ1Z,G2t,Y=σ2Y,G1t,X=σ1X,andE4=X4,Y4,Z4 is the interior equilibrium point of (4).
We shall now show the existence and uniqueness of the solution for allt such that t∈t2,T. Indeed, by referring to (19), we want to show the following items:
(I)
Lipschitz condition: for allW,W/∈ℝ+3andt∈t2,T(23)maxHiW,t−HiW′,t2,GiW,t−GiW′,t2k¯W−W′2,i=1,2,3for some positive constant k¯(II)
Linear growth condition: for allW,t∈ℝ+3×t2,T(24)maxHiW,t2,GiW,t2≤k1+W2,i=1,2,3,for some positive constant k,where W=X,Y,ZT,W/=X/,Y/,Z/T.Proof.
From the inequalities(25)G1W−G1W′2=X−X/2≤σ12X−X′2,G2W−G2W′2=Y−Y/2≤σ22Y−Y′2,G3W−G3W′2=Z−Z/2≤σ32Z−Z′2.
We have(26)H1X−H1X′2=c1X1−X−XY−c1X′1−X′−X/Y2=c1+Y+c1X′+X2X−X′2≤c12+2c1Y∞2+c1X∞+c1X′∞+XY∞+YX′∞+c12X∞2+X/2∞2+2XX′∞2X−X/2,where Φ∞2=supt∈t2,TΦ2.
Thus, we have(27)H1X−H1X′2≤k¯1X−X′2,where
(28)k¯1=c12+2c1Y∞2+c1X∞+c1X/∞+XY∞+YX/∞+c12X∞2+X/2∞2+2XX/∞2.
Similarly(29)k¯2=c12+2c2c4+2c2c3Z∞2+2c3c4X∞,k¯3=c5+c7+c81+c62Y∞2+2c6Y∞+c52Z∞2+Z/∞2+2Z∞Z/∞+2c5c5+c7+c8+2c5c6Y∞Z∞+Z/∞.
Now, by choosingk¯=maxσi2,k¯i,i=1,2,3, we can see that the Lipschitz condition I is satisfied. Similarly, we obtain
(30)G1W2=σ12X−X42≤σ12X421+1X42X2≤σ12X421+X2,with the condition of 1/X42≤1(31)G2W2=σ22Y−Y42≤σ22y421+1Y42Y2≤σ22Y421+Y2,
with the condition of1/Y42≤1, and
(32)G3W2=σ32Z42Z−Z42≤σ321+1Z42Z2≤σ32Z421+Z2,
with the condition ofZ42≤1.
Moreover,≤1+c1X2∞+c1X+Y2∞X2≤1+c1X2∞+c1X+Y2∞X2≤≤≤≤(33)H1X,Y,Z2=1+c1X1−X−XY2,≤1+2c1X+c12X2+c12X4+2c1X3Y+XY2≤1+c1X2∞+c1X+Y2∞X2≤k11+c1X+Y2∞1+c1X2∞X2≤k11+X2,wherek1=1+c1X2∞ with the condition c1X+Y2∞/k1<1. Also
(34)H2X,Y,Z2=c2YZ−c3Y−c4YX2≤c22Z∞+2c3c4X∞+c4X2∞+c32Y2≤k21+Y2,where k2=c22Z∞+2c3c4X∞+c4X2∞+c32, and
(35)H3X,Y,Z2=c5Z1−Z−c6YZ−c7Z+c8XZX+K2≤c5+c8−c72+c52Z∞2+c62+2c5c6+c82+2c2Z∞+2c6c8Y∞Z2≤k31+Z2,where k3=c5+c8−c72+c52Z∞2+c62+2c5c6+c82+2c2Z∞+2c6c8Y∞.
Now, by choosingk=maxσ12X42,σ22Y42,σ32Z42,k1,k2,k3, we can see that the linearity condition II is satisfied. Thus, the system (19) has a unique solution with the initial conditions Xt2=X31,Yt2=Y32,Zt2=Z33. We can then conclude that the piecewise differential equation considered in Scenario 3 has a unique solution for all t∈0,T.
By following the same method, the existence and uniqueness of solutions for the mathematical models described in Scenarios1 and 2 are proved.


## 4. Numerical Approximations

Based on the work studied in [33], we can show that the numerical approximation of the system in Scenario 1 is given as follows:For allt∈0,t1(36)Win+1=WIn+23h12HiWn+5h12HiWn−2−4h3HiWn−1,where
(37)W=X,Y,Z,h=Δt,W1=X,W2=Y,W3=Z,i=1,2,3,andH1,H2,H3are defined in (20).For allt∈t1,T, we have
(38)Xn+1=Xn+23h12H1Xn,Yn,Zn+5h12H1Xn−2,Yn−2,Zn−2−4h3H1Xn−1,Yn−1,Zn−1+512B1tn−1−B1tn−2σ1Xn−2−X4−43B1tn−B1tn−1σ1Xn−1−X4+2312B1tn+1−B1tnσ1Xn−X4,Yn+1=Yn+23h12H2Xn,Yn,Zn+5h12H2Xn−2,Yn−2,Zn−2−4h3H2Xn−1,Yn−1,Zn−1+512B2tn−1−B2tn−2σ2Yn−2−Y4−43B2tn−B2tn−1σ2Yn−1−Y4+2312B2tn+1−B2tnσ2Yn−Y4,Zn+1=Zn+23h12H3Xn,Yn,Zn+5h12H3Xn−2,Yn−2,Zn−2−4h3H3Xn−1,Yn−1,Zn−1+512B3tn−1−B3tn−2σ3Zn−2−Z4−43B3tn−B3tn−1σ3Zn−1−Z4+2312B3tn+1−B3tnσ3Zn−Z4.By referring to Scenario2 and applying the numerical method developed in [34], the approximate solution of the ABC fractional derivative part for all t∈0,t1 is given by
(39)Witn+1=Wi0+1−μFμHiWtn+μFμΓμ∑k=1nHiWtkΓμ+2×hμn+1−kμn−k+2+μ−n−kμn−k+2+2μ−HiWtk−1Γμ+2×hμn+1−kμ+1−n−kμn−j+1+μ,where
(40)W=X,Y,Z,h=Δt,W1=X,W2=Y,W3=Z,i=1,2,3,andH1,H2,H3 are defined in (20).By referring to Scenario3 and based on the work studied in [29], we can show that the numerical approximation of the system for the Caputo-Fabrizio fractional differential equation for t∈0,t1 is given as follows:
(41)Win+1=Win+1−μΓμHiWn−Wn−1+μhΓα512HiWn−2−43HiWn−1+2312HiWn,where
(42)W=X,Y,Z,h=Δt,W1=X,W2=Y,W3=Z,i=1,2,3,andH1,H2,H3are defined in (20).By referring to Scenario3 and based on the work studied in [29], we can show that the numerical approximation of the system for the Caputo fractional differential equation for t∈t1,t2 is given as follows:
(43)Win+1=Wit1+hμΓμ+1∑j=2nHitj−2,Wj−2n−j+1μ−n−jμ+hμΓμ+2∑j=2nHitj−1,Wj−1−Hitj−2,Wj−2×n−j+1μn−j+3+2μ−n−jμn−j+3+3μ+hμ2Γμ+3∑j=2nHitj,Wj−2Hitj−1,Wj−1+Hitj−2,Wj−2×n2n−j2+3μ+10n−j+2μ2+9μ+12−n2n−j2+5μ+10n−j+6μ2+18μ+12,whereW=X,Y,Z,h=Δt,W1=X,W2=Y,W3=Z,i=1,2,3,andH1,H2,H3 are defined in (20).


## 5. Numerical Simulations and Discussion

In this section, some of the numerical simulations of three different scenarios are depicted. Based on the dimensional parameter values given in Table1 and the relationship between parameters given in (5), the following nondimensional parameter values of the system (4) are used for the numerical simulations: c1=0.9,c2=0.65,c3=0.21,c4=0.55,c5=0.52,c6=0.39,c7=0.21,c8=2.48,K=0.2.By using the above parameter values, the interior equilibrium point given by expression (9) is calculated to be E4=X4,Y4,Z4=0.308,3.87,0.580 with the corresponding eigenvalues of the system given by −3.78,−0.22+0.88i,−0.22−0.88i. Thus, the equilibrium point E4 is stable.Let us consider the simulation of Scenario1 for different values of the parameters c5=r2k1/Λ=0.52,0.62,0.92 and c8=ρk1/Λ=0,1.6,2.It can be observed from Figures1–4 that when the values of c5 increase from 0, the hunting cells get activated and seem to increase and continue to oscillate with convergence for the deterministic part and cross-overs to randomness. The hunting cells can increase as indicated in the figures when the value of c5 is different from zero and gets larger and larger. In Figures 1–4, the tumor cells remain mitigated.Figure 1
Time series trajectories of the piecewise systems (10) and (11) for c1=0.9,c2=0.65,c3=0.21,c4=0.55,c5=0,c6=0.39,c7=0.21,c8=2.48,K=0.2.Figure 2
Time series trajectories of the piecewise systems (10) and (11) for c1=0.9,c2=0.65,c3=0.21,c4=0.55,c5=0.52,c6=0.39,c7=0.21,c8=2.48,K=0.2.Figure 3
Time series trajectories of the piecewise systems (10) and (11) for c1=0.9,c2=0.65,c3=0.21,c4=0.55,c5=0.62,c6=0.39,c7=0.21,c8=2.48,K=0.2.Figure 4
Time series trajectories of the piecewise systems (10) and (11) for c1=0.9,c2=0.65,c3=0.21,c4=0.55,c5=0.92,c6=0.39,c7=0.21,c8=2.48,K=0.2.It can be inferred from Figures5–7 that by increasing the parameter values of c8 from zero, the number of tumor cells reduces significantly as shown in both deterministic and stochastic parts of the figures. It means that the proliferation of the resting cells enhanced by the tumor cells characterized by the term c8XZ/X+K (see Equation (4)) has a positive effect on the mitigation of the tumor cells. It can be inferred from this result that proliferation of the hunting cells has a positive effect on mitigating the tumor cells which seems good news for treating the disease.Figure 5
Time series trajectories of the piecewise systems (10) and (11) for c1=0.9,c2=0.65,c3=0.21,c4=0.55,c5=0.52,c6=0.39,c7=0.21,c8=0,K=0.2.Figure 6
Time series trajectories of piecewise systems (10) and (11) for c1=0.9,c2=0.65,c3=0.21,c4=0.55,c5=0.52,c6=0.39,c7=0.21,c8=1.6,K=0.2.Figure 7
Time series trajectories of the piecewise systems (10) and (11) for c1=0.9,c2=0.65,c3=0.21,c4=0.55,c5=0.52,c6=0.39,c7=0.21,c8=2,K=0.2.Some of the phase portraits for Scenario2 which is the case for a cross-over from ABC fractional derivatives to the random process of systems in (15) and (16) are depicted in Figures 8–12. In these figures, the effect of different parameter values of c5 (Figure 8) and the effect of different values of the fractional orders μ (Figures 9–11) are shown.Figure 8
Phase portrait of the piecewise systems (15) and (16) projected onto the YZ plane for the parameter values of c1=0.9,c2=0.65,c3=0.21,c4=0.55,c6=0.39,c7=0.21,c8=2,K=0.2,forthederivativeorderofμ=1, and the different values of c5.Figure 9
Phase portrait of the piecewise systems (15) and (16) projected onto the YZ plane for the parameter values of c1=0.9,c2=0.65,c3=0.21,c4=0.55,c5=0.52,c6=0.39,c7=0.21,c8=2,K=0.2,anddifferentvaluesoffractionalorderμ.Figure 10
Phase portrait of the piecewise systems (15) and (16) projected onto the XY plane for the parameter values of c1=0.9,c2=0.65,c3=0.21,c4=0.55,c5=0.52,c6=0.39,c7=0.21,c8=2,K0.2,anddifferentvaluesoffractionalorderμ.Figure 11
Phase portrait of the piecewise system (15) & (16) projected onto the XYZ plane for the parameter values of c1=0.9,c2=0.65,c3=0.21,c4=0.55,c5=0.52,c6=0.39,c7=0.21,c8=2,K=0.2,and different luesoffractionalorderμ.Figure 12
Phase portrait of the piecewise systems (15) and (16) projected onto the XZ plane for the parameter values of c1=0.9,c2=0.65,c3=0.21,c4=0.55,c5=0.52,c6=0.39,c7=0.21,c8=2,K=0,anddifferentvaluesoffractionalorderμ.It can be observed from Figure8 that by increasing the parameter values of c5, the memory effect of the ABC fractional derivative increases. This effect led to a slow increment of the resting cells Zt, when the number of hunting cells Yt is less than one unit, and later on, as the number of hunting cells increased, the number of resting cells decreases until both of them begin decreasing at the same time as shown in Figure 8. This process repeats, and the trajectories spiral inwards and cross-over to randomness.Some of the simulation results of Scenario3 are described by piecewise systems (17)–(19); a cross-over from fading memory (Caputo-Fabrizio fractional derivative) to the power-law (Caputo fractional derivative) and then to the random process is depicted in Figures 13–15. These figures show the impact of different fractional orders (Figure 13 and zooming in shown in Figure 14 and different values of the parameter c5 (Figure 15 and zooming in shown in Figure 16) on the dynamics of the piecewise system described in Scenario 3.Figure 13
Phase portrait of the piecewise systems (17)–(19) projected onto the XYZ plane for the parameter values of c1=0.9,c2=0.65,c3=0.21,c4=0.55,c5=0.52,c6=0.39,c7=0.21,c8=2,K=0.2,anddifferentvaluesoffractionalorderμ.Figure 14
Zooming in on the inner part of Figure13.Figure 15
Phase portrait of the piecewise systems (15) and (16) projected onto the XZ plane for the parameter values of c1=0.9,c2=0.65,c3=0.21,c4=0.55,c6=0.39,c7=0.21,c8=2,K0.2,differentvaluesoftheparameterc5andfractionalorderμ=0.98.Figure 16
Zooming in on the inner part of Figure15.


## 6. Conclusions

The piecewise mathematical model representation of cancer-immune interaction used in this study has exposed a property that has never been considered or observed in earlier studies using mathematical models based on the classical or different fractional derivatives. The cancer-immune interaction, for instance, showed a cross-over behavior from deterministic to stochastic as shown in the Scenario1 of this study. We argue that the approach of piecewise mathematical model representation of different real-world dynamic systems is an eye-opener for researchers as the approach has the potential of uncovering hidden properties in the dynamics of a system. In this study, it would have been better to compare the model results with the actual data; it is what the researchers are supposed to consider in future work. It can be contested that the piecewise mathematical model approach is better closer to reality as compared to using only one classical or fractional derivative representation of a dynamic system. This is because of this fact that different dynamical systems have the property that cross-over behaviors cannot be cached without a piecewise mathematical model representation of the system. It is observed from the result of this study that (Scenario 1) the proliferation of the resting cells enhanced by the tumor cells characterized by the term c8XZ/X+K is a good target for treatment of the diseases as is shown in Figures 1–5.



---

*Source: 5075613-2022-10-15.xml*