Immunomodulatory role of black tea in the mitigation of cancer induced by inorganic arsenic
Ravi Kiran, Swati Tyagi, Syed Abbas, Madhumita Roy, A. Taraphder
aa r X i v : . [ q - b i o . T O ] M a y Immunomodulatory role of black tea in themitigation of cancer induced by inorganic arsenic
Ravi Kiran ∗ , Swati Tyagi , Syed Abbas † , Madhumita Roy , and A. Taraphder Centre for Theoretical Studies,Indian Institute of Technology Kharagpur, Kharagpur 721302, India Department of Applied Sciences, Punjab Engineering College, Chandigarh-160012, India School of Basic Sciences, IIT Mandi HP 175005 Department of Environmental Carcinogenesis and Toxicology, Chittaranjan National Cancer Institute, 37 S. P.Mukherjee Road, Kolkata 700026, India Department of Physics, Indian Institute of Technology Kharagpur, Kharagpur 721302, India
May 28, 2020
Abstract
We present a model analysis of the tumor and normal cell growth under the influence ofa carcinogenic agent, an immunomdulator (IM) and variable influx of immune cells includingrelevant interactions. The tumor growth is facilitated by carcinogens such as inorganic arsenicwhile the IM considered here is black tea (Camellia sinesnsis). The model with variable influxof immune cells is observed to have considerable advantage over the constant influx model,and while the tumor cell population is greatly mitigated, normal cell population remainsabove healthy levels. The evolutions of normal and tumor cells are computed from theproposed model and their local stabilities are investigated analytically. Numerical simulationsare performed to study the long term dynamics and an estimation of the effects of variousfactors is made. This helps in developing a balanced strategy for tumor mitigation withoutthe use of chemotherapeutic drugs that usually have strong side-effects.
Keywords : Tumor growth dynamics; Immune response; Mathematical model; Next gener-ation matrix; Basic reproduction number; Stability.
AMS subject classification : 93A30; 93D20; 37B25; 37N25.
The scourge of cancer is on the rise all around the globe. There have been a number of factorsresponsible for this increasing trend. Cancer is one of the leading causes of death and accordingto 2018 GLOBOCAN data, death due to cancer stands at 9.6 million in 2018 [1]. Canceris an uncontrolled and abnormal proliferation of cells, leading to formation of tumour whichcan infiltrate and destroy normal tissues. There have a number of etiologic factors responsiblefor this disease, including environmental exposures, lifestyle and others. Often there is someoccupational exposure to certain carcinogenic risk-factors like chemicals, radioactive materialsetc. ∗ [email protected], [email protected] † [email protected] In this section we describe the model in detail. As discussed above, we do not use chemotherapyto destroy the cancer cells and instead include the effects of BT in reducing tumor growth (albeitthrough the reduction of ROS) and the modulation of immunity.2 .1 The Model —Overview
The model presented has following components:1.
Arsenic as stimulator of tumor: iAs as an external source to stimulate the growth oftumor. We emphasize that while we consider iAs as driving the tumor growth, based onthe experiments discussed above, it can, in general, be viewed as any carcinogen presentin the environment to which exposure is common.2.
Immune Response:
The model includes immune cells whose growth may be stimulated bythe presence of tumor. These cells can inhibit the growth of tumor cells. Their effects aretaken in the model via the usual kinetic process. The important addition to usual modelsis a variable influx, s(t) in the evolution of immune cells. This reflects the fact that theimmune response is ideally never constant and the body can produce variable amounts ofimmune cells at different stages of response.3.
Competition Terms:
Normal cells and tumor cells compete for available resources, whileimmune cells and tumor cells compete in predator-prey fashion.4.
BT-Response:
BT has mitigating effects on cancer growth (via ROS reduction) and there-fore a term representing such effects is included [6]. The model also includes BT-tumorinteraction, as well as BT-immune cell interaction. The beneficial effect of BT is seenthrough the suppression of tumor cells and strengthening of the immune cells, and it hasno known adverse effects on the normal cells unlike chemotherapy.
Both the normal and tumor cells independently increase according to the usual logistic growthlaw. The interaction between normal and tumor cells is of predator-prey type, described by thefollowing system of equations, where normal cells are denoted by N and the tumor cells aredenoted by T. dNdt = r N (1 − b N ) − c T NdTdt = r T (1 − b T ) − c ′ T N
The interaction terms c and c ′ are competition terms and are both assumed to be positive. Anegative competition term would imply that instead of the normal cells destructively competingwith the tumor cells for resources and space, the presence of the normal cells would in factstimulate further growth of the tumor cell population. While some authors argue that c ′ couldbe negative [12, 13], we assume a destructive competition in this study.The inclusion of arsenic in the system has deleterious effects on normal cells, and a certainproportion of normal cells become tumorous. It is represented by the following terms dNdt = − α ′ AN and dTdt = αAN We let N ( t )( ≡ N ) denote the number of normal cells at time t , T ( t )( ≡ T ) denote the number of tumor cellsat time t , I ( t )( ≡ I ) denote the number of immune cells at time t , and s ( t )( ≡ s ) denote the variable influx ofimmune cells. A ( t )( ≡ A ) and D ( t )( ≡ D ) are the Arsenic and BT respectively.
3t is possible that not all of the normal cells turn into tumorous, as iAS can also result inthe death of normal cells. In earlier studies, a constant influx rate has been assumed for theimmune cells. The assumption has been relaxed here to accommodate real immune responses,which is variable, and along with a steady rate, production of immune cells at a rate of η isallowed. This is represented by the term ηb + s s To avoid immune cell proliferation and immune-upon-immune crowding, a saturation valueis assumed. The presence of tumor cells stimulates the immune response, represented by thepositive nonlinear growth term for the immune cells: ρITα + T where ρ and α are positive constants. This type of response terms is of the same form as usedin the models of Kuznetsov et.al. [14]. Moreover, the reaction of immune cells and tumor cellscan result in either the death of tumor cells or the deactivation of the immune cells, representedby two competition terms dTdt = − βIT and dIdt = − β ′ IT Arsenic and BT considered in the model have a simple form: there is a steady influx ofboth and a decay with certain rates. The presence of BT is assumed to stimulate the immuneresponse and hence a term similar to Michaelis-Menten equation is assumed. BT also has anadverse effect on the tumor, and both these effects are represented: dTdt = − γDT and dIdt = δIDα + D Combining all these terms, we propose and analyze the model described by the following systemof equations: dNdt = r N (1 − b N ) − c T N − α ′ ANdTdt = r T (1 − b T ) − c ′ T N + αAN − βIT − γDTdIdt = s ( t ) + ρITα + T − d I I − β ′ IT + δIDα + D (1) dAdt = a − d A AdDdt = b − d D Ddsdt = s + ηb + s s − µ s. In this section, we summarize the parameters of the mathematical model presented. This pa-rameter set may vary depending on the case one is analysing. However, the analyses of themodel are quite general and hence they still apply. The idea is to keep the tumor population aslow as possible with normal cell count above a healthy threshold. We start with a small amountof tumor cells and immune cells while initial values of arsenic and BT are zero.4
Per unit growth rates: r and r are growth rates for tumor cells and normal cells respec-tively. Here, we assume the tumor cell population grows more rapidly than the normalcell population, and let r > r . • Carrying capacities: b − ≤ b − = 1 . • Competition terms: c , c ′ , α, α ′ , β β ′ and γ are competition terms. • Death rates: d I , d A and d D are per capita death rates of immune cells, arsenic, BTrespectively; µ is the death rate of the stimulated immune cells. a and b are constants,assumed 0.4 here. • Immune source rate:
The immune source rate is considered to be variable here, denotedby s ( t ). The influx of immune cells is to be stimulated from outside; s is the constantinflux; η is the production rate, b is the saturation term and µ is decay constant. • Immune response rate: ρ, which is assumed to have a baseline value of 1. In the studybelow, we set ρ = 0 .
01 to simulate a compromised immune system. • Immune threshold rate: α , is inversely related to the steepness of the immune responsecurve. • BT-Immune term:
The BT-immune term is modelled similar to Michaelis-Menten term,with α as threshold rate and δ as response rate.We first analyse the equations analytically for positivity, boundedness, equilibrium pointsand the stability of the solutions. We also look at possible instabilities. Thereafter, we work outthe dynamics and relevant phase diagram numerically. The system (1) has initial conditions given by N (0) = N ≥ , T (0) = T ≥ , I (0) = I ≥ ,A (0) = A ≥ , D (0) = D ≥ , and s (0) = s ≥ . We consider all the variables and parametersof the model to be non-negative, since the model investigates cellular populations. Based on thebiological findings, we analyze the system (1) in the region Ω = (cid:8) ( N, T, I, A, D, s ) ∈ R (cid:9) . Wefirst assure that the system (1) is well posed such that the solutions with non-negative initialconditions remain non-negative for all 0 < t < ∞ and thus making the variables biologicallymeaningful. Theorem 3.1.
The region Ω ⊂ R given by Ω = (cid:26) ( N, T, I, A, D, s ) ∈ R : N ≤ b (cid:27) is posi-tively invariant with respect to system of equations (1) and non-negative solutions exist for alltime < t < ∞ . Proof.
Let Ω ⊂ R given by Ω = (cid:26) ( N, T, I, A, D, s ) ∈ R : N ≤ b (cid:27) . Then the solutions( N ( t ) , T ( t ) , I ( t ) , A ( t ) , D ( t ) , s ( t )) of (1) are positive for all t ≥ . It can be observed from thefirst compartment that dNdt = r N − b N − c T N − α ′ AN ≤ r N − b N . N (0) = N , we obtain N ( t ) ≤ b + ke − r t , (2)where k = 1 − N b N . Thus N = 1 b + k . Substituting the value of k in (2), we obtain N ( t ) ≤ b + − N b N e − r t ≤ b as t → ∞ . Since b > , we have N ( t ) > t > . Similarly we can show that T ( t ) > , I ( t ) > , A ( t ) > , D ( t ) > s ( t ) > . In this section, we discuss the existence of all possible equilibrium points of system (1). Themodel system (1) admits four equilibrium points, one tumor-free equilibrium point, two deadequilibrium points and one co-existing equilibrium point P ∗ = ( N ∗ , T ∗ , I ∗ , A ∗ , D ∗ , s ∗ ) respec-tively. We have N ∗ > , T ∗ > , I ∗ > , A ∗ > , D ∗ > , s ∗ > • Tumor-Free equilibrium point (TFE): P = ( 0 . α ′ − d A r . α ′ b , , s ∗ ( d D α + 0 . d ( d D α + 0 . − . δ , . d A , . d D , s ∗ ) , where s ∗ is the real positive so-lution of equation µ s + s ( µ b − η − s ) − bs = 0 . • Dead equilibrium point: An equilibrium point is said to be dead equilibrium, if the normalcell population is zero. For d α > ( δ − d ) 0 . d D , we have two dead equilibria in this case.1. Type 1 Dead equilibrium point P d = (0 , , s ∗ ( d D α + 0 . d ( d D α + 0 . − . δ , . d A , . d D , s ∗ ) .
2. Type 2 Dead equilibrium point P d = (0 , m ∗ , s ∗ ( d D α + 0 . d ( d D α + 0 . − . δ , . d A , . d D , s ∗ ) . Here m ∗ is the real positive solu-tion of equation T + AT + BT + C = 0 , where A = 1 β ′ r b (cid:18) − r b ρ + d r b + β ′ α r b − β ′ r − β ′ γD ∗ + δD ∗ r b α + D ∗ (cid:19) B = 1 β ′ r b ( s ∗ β + r ρ + ργD ∗ + d α r b − d r − d γD ∗ − β ′ α r − β ′ α γD ∗ − δαD ∗ r b α + D ∗ + δD ∗ ( r + γD ∗ ) α + D ∗ ) C = 1 β ′ r b (cid:18) s ∗ βα − d α r − d α γD ∗ + δD ∗ α ( r + γD ∗ ) α + D ∗ (cid:19) D ∗ = 0 . d D . • Co-existing equilibrium point P ∗ = ( N ∗ , T ∗ , I ∗ , A ∗ , D ∗ , s ∗ )6 The Reproduction Number and Tumor-Free Equilibrium Point
In this section, we find the basic reproduction number by following the next generation matrixmethods as described in [15, 16]. At TFE, the matrices F and V are given as follows: F = (cid:18) r αN ∗ (cid:19) and V = (cid:18) c ′ N ∗ + βI ∗ + γD ∗ d A (cid:19) . Then the basic reproduction number R is given by the largest eigenvalue of F V − . Thus,we have R = r r b d A [ d ( α d D + 0 . − . δ ]( d ( α d D + 0 . − . δ ) ( r d A c ′ − . α ′ c ′ + 0 . γb r ) + b r d A βs ∗ ( α d D + 0 .
4) (3)The Basic reproduction number R , which measures the rate of spread of tumor; R > , if eachcell produces on average more than one cell and thus the tumor grows over time. If R < , theneach cell produces on average less than one new cell and thus the therapy (drug administration)can eradicate the tumor. At each time step, a tumor cell either produces an offspring or dies. Theorem 6.1.
The Tumor-free equilibrium point P ∗ = ( N ∗ , , I ∗ , A ∗ , D ∗ , s ∗ ) is locally asymp-totically stable, if the following Routh-Hurwitz criterion is satisfied, B + B < B B ( B + B ) + B B B < . along with bη ( b + s ∗ ) < µ , otherwise unstable, where B ′ i s are as defined in (5).Proof. Linearizing system (1) around TFE P , we obtain the following Jacobian matrix J ( P ) J ( P ) = B B B B B B B B
10 0 0 B B
00 0 0 0 0 B where B = − c ( r d A − . α ′ ) b r d A , B = − α ′ ( r d A − . α ′ ) b r d A (5) B = 0 . αd A , B = r − c ′ N ∗ − βI ∗ − γD ∗ B = α ( r d A − . α ′ ) b r d A , B = ρα − βI ∗ B = 0 . δ . α d D − d , B = α δI ∗ ( α + D ∗ ) B = − d A , B = − d D , B = bη ( b + s ∗ ) − µ . (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − λ B B B B − λ B B B − λ B
10 0 0 B − λ B − λ
00 0 0 0 0 B − λ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = 0Solving this determinant, we obtain the following equation:( λ − B )( λ − B )( λ − B ) (cid:0) λ + A λ + A λ + A (cid:1) = 0 . (6)where A = − ( B + B ) , A = B B and A = B B B . Clearly the characteristic equation (6)has three roots given by λ = B = − d A < , λ = B = − d D < , λ = B = bη ( b + s ∗ ) − µ < . We are left with the cubic equation λ + A λ + A λ + A = 0 . Applying Routh-Hurwitzcriterion, the Tumor-free equilibrium is locally asymptotically stable, provided A > A A − A > . Theorem 6.2.
The Type 1 Dead equilibrium point P d = (0 , , I ∗ , A ∗ , D ∗ , s ∗ ) is locally asymp-totically stable provided the following hold; r < . α ′ d A , µ > bη ( b + s ∗ ) , r < βI ∗ + γD ∗ , d > . δ . α d D . (7) Proof.
The characteristic equation corresponding to P d is given by( r − α ′ A ∗ − λ )( r − βI ∗ − γD ∗ − λ )( − d + δD ∗ α + D ∗ − λ )( − d A − λ )( − d D − λ )( bη ( b + s ∗ ) − µ − λ ) = 0 . (8)Clearly the roots of characteristic equation (8) are r − α ′ A ∗ < , r − βI ∗ − γD ∗ < , − d + δD ∗ α + D ∗ < , − d A < , − d D < bη ( b + s ∗ ) − µ < . Thus the dead equilibrium P d is locally asymptot-ically stable. Theorem 6.3.
The co-existing equilibrium point P ∗ = ( N ∗ , T ∗ , I ∗ , A ∗ , D ∗ , s ∗ ) is locally asymp-totically stable, if the following Routh-Hurwitz criterion is satisfied, A + A + A < A + A )( A A + A A + A A + A − A A ) + A A < . otherwise unstable, where A ′ i s are as defined in (10)Proof. The characteristic equation of system (1) corresponding to P ∗ is given by (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) A − λ A A A A − λ A A A A A − λ A
10 0 0 A − λ A − λ
00 0 0 0 0 A − λ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = 0 , A = r (1 − b N ∗ ) − c T ∗ − α ′ A ∗ , A = − c ′ T ∗ + αA ∗ (10) A = − c N ∗ , A = r (1 − b T ∗ ) − c ′ N ∗ − βT ∗ A = α ρI ∗ ( α + T ∗ ) − β ′ I ∗ , A = − βT ∗ A = ρT ∗ α + T ∗ − d − β ′ T ∗ + δD ∗ ( α + D ∗ ) , A = − α ′ N ∗ , A = αN ∗ ,A = − d A , A = − γT ∗ A = δα I ∗ ( α + D ∗ ) A = − d D , A = bη ( b + s ∗ ) − µ . Thus we obtain the following equation( A − λ )( A − λ )( A − λ ) (cid:0) λ + B λ + B λ + B (cid:1) = 0 , (11)where B = − ( A + A + A ) , B = A A + A A + A A − A A and B = A A + A A A − A A A . Clearly the characteristic equation (11) has three roots given by λ = A = − d A < , λ = A = − d D < λ = A = bη ( b + s ∗ ) − µ < . We are left with the cubic equation λ + B λ + B λ + B = 0 . Applying Routh-Hurwitz criterion, the co-existing equilibrium P ∗ islocally asymptotically stable, provided B > B B − B > . This implies the co-existingequilibrium P ∗ is locally asymptotically stable if A + A + A < A + A )( A A + A A + A A + A − A A ) + A A < . Irrespective of the fact that tumor-free equilibrium point is stable or unstable, efforts of doctorshave been always oriented to reach to this point. As tumor cells can be completely destroyedat this point, thus resulting in complete therapy of the disease. Moreover, in the mean time, itbecomes necessary to find a therapeutic protocol to be able to incline the solution of equationstowards this stable point, irrespective of the initial conditions. For this, we need stimulationof immune cells and drug administration, which could guarantee the global stability of thisequilibrium point.In this section, we employ Lyapunov’s direct method [17] to design the desirable diseaseeradication protocol. This technique requires selecting a suitable Lyapunov function candidateand then finding a control law to make this candidate a real Lyapunov function.
Theorem 7.1.
The tumor-free steady state, P = ( 0 . α ′ − d A r . α ′ b , , s ∗ ( d D α + 0 . d ( d D α + 0 . − . δ , . d A , . d D , s ∗ )= ( N ∗ , , I ∗ , A ∗ , D ∗ , s ∗ ) is globally asymptotically stable if: ( i ) b < µ η , ( ii ) g = max n b N ∗ , β ′ I − β ′ II ∗ − α ρI αAN ∗ − αA ∗ N ∗ o ( iii ) δα > S − S ∗ D ∗ I + DI ∗ − D ∗ I ∗ roof. We define the following Lyapunov function, V ( t ) = a s − s ∗ ) + h A − A ∗ ) + c D − D ∗ ) + g N − N ∗ ) + f I − I ∗ ) + m T , where a, h, c, g, f, m are all positive constants. Computing time derivative of V along with (1),we obtain dVdt = a ( s − s ∗ )( s + ηb + s s − µ s ) + h ( A − A ∗ )(0 . − d A A ) + c ( D − D ∗ )(0 . − d D D )+ g ( N − N ∗ )( r N (1 − b N ) − c T N − α ′ AN )+ f ( I − I ∗ )( s ( t ) + ρITα + T − d I I − βIT + δIDα + D )+ mT ( r T (1 − b T ) − c ′ T N + αAN − βIT − γDT ) ≤ − bd A ( A − A ∗ ) − − cd D ( D − D ∗ ) − aµ ( s − s ∗ ) + aηb ( s − s ∗ ) ( b + s )( b + s ∗ ) + ( s − s ∗ )( I − I ∗ ) − s ∗ I ∗ ( I − I ∗ ) + α ρIT ( I − I ∗ ) α ( α + T ) − β ′ IT ( I − I ∗ ) + δα I ( I − I ∗ )( D − D ∗ )( α + D )( α + D ∗ ) + gr ( N − N ∗ ) − gr b ( N − N ∗ ) + αgr b N ∗ ( N − N ∗ ) − c gN T ( N − N ∗ ) − c gT ∗ ( N − N ∗ ) − gα ′ N ( A − A ∗ )( N − N ∗ ) − gα ′ A ∗ ( N − N ∗ ) − gr b N N ∗ ( N − N ∗ ) + gr T − gr bT − gc ′ T ( N − N ∗ ) − gc ′ N ∗ T + gαAT ( N − N ∗ ) + gαT N ∗ ( A − A ∗ ) − gβIT − gγT ( D − D ∗ ) − gγD ∗ T . Then using (i), (ii) and (iii), we can obtain dVdt < . Thus one can guarantee that solution ofequations goes to tumor free equilibrium point, if the parameters of model system (1) satisfy(i)-(iii).
Corollary 7.2.
The co-existing steady state, P ∗ = ( N ∗ , T ∗ , I ∗ , A ∗ , D ∗ , s ∗ ) is globally asymptot-ically stable, provided the following conditions hold: ( i ) b < µ η , ( ii ) g = max n b N ∗ , b T ∗ , βT ∗ ( s ∗ − s + β ′ I ∗ T ∗ − β ′ IT ∗ − α ρI ∗ T ∗ + α ρIT ∗ − δα DI − δα ID ∗ + δα DI ∗ − δα D ∗ I ∗ ) , β ′ I − β ′ II ∗ − α ρI + α ρIT ∗ αAN ∗ − αA ∗ N ∗ − βIT ∗ + βI ∗ T ∗ , − s ∗ I ∗ − sI ∗ αN ∗ A ∗ T ∗ − αN ∗ AT ∗ − βI ∗ T ∗ o for some positive constant g. In comparison to non-delayed models, delay differential equation (DDEs) systems can exhibitmuch richer dynamics since a time delay could cause the loss of stability of equilibrium and giverise to periodic solutions through the Hopf bifurcation. The instabilities and oscillatory behaviorcaused by delays are very common, however the delays may also have the opposite effect, namelythat they can suppress oscillations and stabilize equilibria which would be unstable in the absenceof delays. 10ue to chemical transportation of signals and the time needed for differentiation/divisionof cells, the production of tumor cells may not be instantaneous but, instead, it exhibits sometime lag. Moreover, due to the immunity, there may be time delay in competition of normalcells and tumor cells as well. To capture such phenomenon, we introduce two delays into thenon-delayed model (1). Thus we obtain the following system of DDEs corresponding to system(1) described by the following equations; dNdt = r N (1 − b N ) − c T N − α ′ A ( t − τ ) N ( t − τ ) dTdt = r T (1 − b T ) − c ′ T N − βIT − γDT + αA ( t − τ ) N ( t − τ ) dIdt = s ( t ) + ρITα + T − d I I − βIT + δIDα + D (12) dAdt = 0 . − d A AdDdt = 0 . − d D Ddsdt = s + ηb + s s − µ s. We denote by C the Banach space of continuous functions φ : [ − τ, → R equipped withsuitable norm, where τ = max { τ , τ } . Further let C + = { φ = ( φ , φ , φ , φ , φ , φ ) ∈ C : φ i ( θ ) ≥ ∀ θ ∈ [ − τ, , i = 1 , , · · · , . The initial conditions corresponding to delayed system(12) are N ( θ ) = φ ( θ ) , T ( θ ) = φ ( θ ) , I ( θ ) = φ ( θ ) , A ( θ ) = φ ( θ ) , D ( θ ) = φ ( θ ) , s ( θ ) = φ ( θ ) , (13)where φ = ( φ , φ , φ , φ , φ , φ ) ∈ C + . Qualitative Analysis: Preliminaries
In this subsection, we establish the non-negativity of the solutions of system (12) with initialconditions (13).
Proposition 8.1.
The solutions ( N ( t ) , T ( t ) , I ( t ) , A ( t ) , D ( t ) , s ( t )) with initial conditions (13) ofthe system (12) are non-negative.Proof. We can rewrite the system (12) in vector form by setting Z = ( N, T, I, A, D, s ) T ∈ R and F ( Z ) = F ( Z ) F ( Z ) F ( Z ) F ( Z ) F ( Z ) F ( Z ) = r N (1 − b N ) − c T N − α ′ A ( t − τ ) N ( t − τ ) r T (1 − b T ) − c ′ T N − βIT − γDT + αA ( t − τ ) N ( t − τ ) s ( t ) + ρITα + T − d I I − βIT + δIDα + D . − d A A . − d D Ds + ηb + s s − µ s. (14)where F : C + → R and F ∈ C ∞ ( R ) . Then delayed model (12) becomes˙ Z ( t ) = F ( Z t ) , (15)where ˙ ≡ ddt , Z t ( θ ) = Z ( t + θ ) , θ ∈ [ − τ, . It can be observed from (15) that whenever we choose Z ( θ ) ∈ C + such that Z i = 0 , we obtain F i ( z ) | Z i ( t )=0 , Z + ∈ C + ≥ , i = 1 , , · · · , . Then from1118], any solution (15) with Z t ( θ ) ∈ C + say Z ( t ) = Z ( t, Z (0)) is such that Z ( t ) ∈ R + ∀ t > . Further we define P ( t ) = N ( t ) + T ( t ) + I ( t ) + A ( t ) + D ( t ) + s ( t ) . Then using the similar analysis as done in Theorem 3.1, we have P ( t ) is bounded and hence are N ( t ) , T ( t ) , I ( t ) , A ( t ) , D ( t ) , s ( t ) . This completes the proof.
Stability and Bifurcation Analysis
To investigate the local stability of equilibria of system (12), we linearize the system and evaluatethe characteristic equation first at equilibrium P ∗ . The characteristic equation is (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) A − λ A A A A − λ A A A A A − λ A
10 0 0 A − λ A − λ
00 0 0 0 0 A − λ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = 0 , where A = r (1 − b N ∗ ) − c T ∗ − α ′ Ae − λτ = P − α ′ Ae − λτ ,A = − c ′ T ∗ + αAe − λτ = Q + αAe − λτ (16) A = − c N ∗ , A = r (1 − b T ∗ ) − c ′ N ∗ − βT ∗ A = α ρI ∗ ( α + T ∗ ) − β ′ I ∗ , A = − βT ∗ A = ρT ∗ α + T ∗ − d − β ′ T ∗ + δD ∗ ( α + D ∗ ) ,A = − α ′ N e − λτ , A = αN e − λτ ,A = − d A , A = − γT ∗ A = δα I ∗ ( α + D ∗ ) A = − d D , A = bη ( b + s ∗ ) − µ . (17)The characteristic equation is P ( λ ) + P ( λ ) e − λτ + P ( λ ) e − λτ = 0 , (18)where λ is an eigenvalue and P ( λ ) = λ + λ ( − P − A − A ) + λ ( A P + A P + A A − QA ) + ( A A + QA A − P A A )= λ + H λ + H λ + H P ( λ ) = − αAA λ + αA A = H λ + H P ( λ ) = α ′ Aλ + α ′ AA A = H λ + H . Case 1: τ = τ = 0This case is equivalent to (11) of non-delayed model. Thus P ∗ is locally asymptotically stableprovided (9) holds. Similarly, TFE is locally asymptotically stable provided (4) holds and Type1 dead equilibrium is locally asymptotically stable, if (7) holds.12 ase 2: τ > , τ = 0In this case, the characteristic equation becomes P ( λ ) + P ( λ ) e − λτ = 0 . We can rewrite thisequation as (cid:0) λ + H λ + H λ + H (cid:1) + ( H λ + H ) e − λτ = 0 . (19)If the time delay τ is able to destablize P ∗ and produces oscillations, this can occur only whencharacteristic roots cross the imaginary axis to the right. Let ω > λ = iω is a purelyimaginary root of (19). Separating real and imaginary parts, we have H cos ωτ + H ω sin ωτ = H ω − H (20) H ω cos ωτ − H sin ωτ = ω − H ω. (21)Eliminating τ in (20) and (21), we obtain the following sixth-degree polynomial equation F ( ω ) = ω + h ω + h ω + h = 0 , (22)where h = H − H , h = H − H H − H and h = H − H . Letting v = ω gives thefollowing simplified system G ( v ) = v + h v + h v + h = 0 . (23)Define ∆ = h − h . Next we follow the method described in [19] to investigate for the existenceof positive roots of equation (22). From (23), we have dG ( v ) dv = 3 v + 2 h v + h (24)If h < , since lim v →∞ G ( v ) = ∞ , then we obtain that (22) has atleast one positive root.When ∆ ≤ , dG ( v ) dv ≥ , that is G ( v ) is monotonically increasing. Thus, if h ≥ ≤ , then (22) has no positive root.When ∆ > , then graph of G ( v ) has two critical points v ∗ = − h + √ ∆3 , v ∗ = − h − √ ∆3 . Therefore, if h ≥ , then from Lemma 2.2 of [19], we can say that (24) has positive roots if andonly if ∆ > , v ∗ > G ( v ∗ ) ≤ . Assume that equation (24) has three positive roots, givenby v , v , v respectively. Then (24) has three positive roots ω k = √ v k , k = 1 , , . Solving (20)and (21) for τ , we obtain τ ( n )1 ,k = 1 ω k arccos H ( H ω ) − H ) + H ω ( ω − H )2( H + ω H ) + 2 nπω k , k = 1 , , , · · · , n = 0 , , , · · · (25)and ± iω k is pair of purely imaginary roots of (22) with τ ( n )1 ,k . We further define τ ∗ = τ (0)1 ,k = min k = { , , } n τ (0)1 ,k o , ω ∗ = ω k . (26)We now obtain the transversality condition for Hopf Bifurcation at τ = τ ∗ . Differentiating (19)with respect to τ and substituting expression for e − λτ from (19), we obtain dλdτ h λ + 2 λH + H − ( H − H λτ − H τ )( λ + H λ + H λ + H ) H λ + H i = − λ ( λ + H λ + H λ + H )Thus we have (cid:18) dλdτ (cid:19) − = H λ ( H λ + H − λ + 2 λH + H λ ( λ + H λ + H λ + H ) − τ λ . (cid:0) dλdτ (cid:1) − at τ = τ ∗ ( i.e.λ = iω ∗ ) and taking real part, we obtain Re "(cid:18) dλdτ (cid:19) − (cid:12)(cid:12)(cid:12) τ = τ ∗ = 3 ω + 2 ω ( H − H ) + ( H − H H − H ) H ω + H = G ′ ( ω ∗ ) H ω ∗ + H We have G ( v ) is non-increasing and positive. Thus sign (cid:26)(cid:18) dRe ( λ ) dτ (cid:19) (cid:12)(cid:12)(cid:12) τ = τ ∗ (cid:27) = sign ((cid:18) dλdτ (cid:19) − (cid:12)(cid:12)(cid:12) τ = τ ∗ ) = sign { G ′ ( ω ∗ ) } . Hence if G ′ ( ω ∗ ) = 0 , then transversality condition holds. Summarizing the above, we have thefollowing theorem Theorem 8.2.
For τ > , τ = 0 , assume that condition (9) holds. If either h < or h ≥ , ∆ > , v ∗ > and G ( v ∗ ) < , then the Endemic equilibrium P ∗ is locally asymptotically stablefor < τ < τ ∗ , where τ ∗ = 1 ω ∗ arccos H ( H ω ) − H ) + H ω ( ω − H )2( H + ω H ) + 2 nπω ∗ . (27) Furthermore, if G ′ ( ω ∗ ) = 0 , then system (12) undergoes Hopf bifurcation to periodic solutionsat P ∗ at τ = τ ∗ . Remark 8.3. If h ≥ and ∆ ≤ , then equation (24) has no positive real root, thus theequilibrium P ∗ is locally asymptotically stable for all τ > . Case 3: τ > , τ = 0In this case, the characteristic equation becomes P ( λ ) + P ( λ ) e − λτ = 0 . We can rewrite thisequation as (cid:0) λ + H λ + H λ + H (cid:1) + (cid:0) H λ + H (cid:1) e − λτ = 0 . (28)Using the similar analysis as done in case 2, we have the following theorem. Theorem 8.4.
For τ > , τ = 0 , assume that condition (9) holds. If either h < or h ≥ , ∆ > , v ∗ > and G ( v ∗ ) < , then the Endemic equilibrium P ∗ is locally asymptotically stablefor < τ < τ ∗ , where τ ∗ = 1 ω ∗ arccos ( H ω − H )2( H − ω H ) + 2 nπω ∗ (29) Furthermore, if G ′ ( ω ∗ ) = 0 , then system (12) undergoes Hopf bifurcation to periodic solutionsat P ∗ at τ = τ ∗ . Case 4: τ > , τ ∈ (0 , τ ∗ )In this case, we consider τ as a parameter and fix τ at a point in its stable interval. At P ∗ , thecharacteristic equation takes the following form; P ( λ, τ , τ ) = (cid:0) λ + H λ + H λ + H (cid:1) + ( H λ + H ) e − λτ + (cid:0) H λ + H (cid:1) e − λτ = 0 . (30)Assume that (30) has purely imaginary root given by λ = iω. Substituting in (30) and separatingreal and imaginary parts, we obtain H ω − H − H cos ωτ − H ω sin ωτ = ( H − H ω ) cos ωτ (31) − ω + H ω + H ω cos ωτ − H sin ωτ = ( H − H ω ) sin ωτ . (32)14liminating τ in (31) and (32), we obtain the following polynomial equation; ω + ω ( H − H − H ) + ω (2 H sin ωτ − H H sin ωτ )+ ω ( H − H H + H + 2 H H − H H cos ωτ ) = 0 , (33)which is a sixth-degree polynomial in ω. If (C1) holds, then let we denote the six positive rootsof (33) as ω k , k = 1 , , · · · , . Solving (31)-(32) for τ , we obtain τ n ,k = 1 ω k arccos H ω − H − H cos ωτ − H ω sin ωτ H − H ω + 2 nπω k , k = 1 , , · · · , n = 0 , , · · · (34)and ± ω k is pair of purely imaginary root of (33) with τ n ,k . For simplicity we define τ ∗ = τ ,k = min k =1 , , ··· , (cid:8) τ ,k (cid:9) , ω ∗ = ω k . (35)Now, differentiating (30) with respect to τ and further simplifying for transversality condition,we obtain Re "(cid:18) dλdτ (cid:19) − (cid:12)(cid:12)(cid:12) τ = τ ∗ = 2 ωH cos ωτ + ( H − ω ) sin ωτ ω ( H − H ω ) + H sin ω ( τ − τ ) ω ( H − H ω ) + 2 H ( H − H ω ) (cid:12)(cid:12)(cid:12) τ ∗ . Since (C2) holds, the transversality condition holds. Thus we have the following result.
Theorem 8.5.
Assume condition (C1)-(C2) hold and τ ∈ (0 , τ ∗ ) , then P ∗ is locally asymptot-ically stable for τ ∈ (0 , τ ∗ ) and system (12) undergoes Hopf bifurcation to periodic solutions at P ∗ for τ = τ ∗ , where(C1): Equation (33) has six positive roots,(C2): H > H ω (cid:12)(cid:12)(cid:12) τ ∗ , H > ω (cid:12)(cid:12)(cid:12) τ ∗ . Case 5: τ > , τ ∈ (0 , τ ∗ )In this case, τ is fixed at a point in its stable interval and τ is considered as a parameter. At P ∗ we have the same characteristic equation as (30). Following the similar procedure as in case4, we substitute λ = iω and separate the real and imaginary parts. Thus we obtain H cos ωτ + H ω sin ωτ = H ω − H − ( H − H ω ) cos ωτ (36) H ω cos ωτ − H sin ωτ = ω − H ω + ( H − H ω ) sin ωτ . Thus we obtain ω − (2 H sin ωτ ) ω + ( H + H + 2 H H cos ωτ ) ω + (2 H H sin ωτ + 2 H sin ωτ ) ω (37)+ ( − H H − H H − H H cos ωτ − H H cos ωτ + H − H ) ω − (2 H H sin ωτ ) ω + ( H + H − H H cos ωτ − H ) = 0 (38)and τ n ,k = 1 ω k arccos H ( H ω − H − ( H − H ω ) cos ωτ ) + H ω ( ω − H ω + ( H − H ω ) sin ωτ )2( H + H ω ) + 2 nπω k , (39) k = 1 , , · · · , n = 0 , , · · · . urthermore, we have sign h Re "(cid:18) dλdτ (cid:19) − (cid:12)(cid:12)(cid:12) τ = τ ∗ = sign h cos ωτ (3 H ω − H H ω + 2 H H ω ) − sin ωτ (3 ω H − H H − H H ω ) + H − H ω sin ω ( τ − τ ) ω ( H + H ω ) i . Theorem 8.6.
Assume conditions (C1’)-(C2’) hold and τ ∈ (0 , τ ∗ ) , then P ∗ is locally asymptoticallystable for τ ∈ (0 , τ ∗ ) and system (12) undergoes Hopf bifurcation to periodic solutions at P ∗ for τ = τ ∗ , where(C1’): Equation (37) has six positive roots,(C2’): H ω + 2 ωH H > H H ω, ω H < H H + 2 H H ω and H ω sin ω ( τ − τ ) < . Case 6: τ = τ = τ In this case, characteristic equation is P ( λ, τ ) = (cid:0) λ + H λ + H λ + H (cid:1) + (cid:0) H λ + H λ + H (cid:1) e − λτ = 0 , (40)where H = H + H . Then substituting λ = ı ω in (40) and following the similar procedure, we obtain τ ik = 1 ω k arccos H ( H ω − H ) + ( H − H ω )( ω − H )2 H ( H − H ω ) + 2 nπω k , k = 1 , , · · · , n = 0 , , · · · . (41)and ω + ( H − H − H ) ω + ( H − H H + 2 H H ) + ( H − H ) = 0 . (42)Furthermore, we have sign " Re "(cid:18) dλdτ (cid:19) − (cid:12)(cid:12)(cid:12) τ = τ ∗ = sign " ( H − H ω ) (cid:0) H ω cos ωτ + ( H − ω ) sin ωτ (cid:1) − H ω (cid:0) ( H − ω ) cos ωτ − H ω sin ωτ (cid:1) + 2 ωH ( H − H ω ) − H ωω (cid:0) ( H ω − H ) + H ω (cid:1) . Theorem 8.7.
Assume conditions (C1”)-(C2”) hold and τ ∈ (0 , τ ∗ ) , then P ∗ is locally asymptoticallystable for τ ∈ (0 , τ ∗ ) and system (12) undergoes Hopf bifurcation to periodic solutions at P ∗ for τ = τ ∗ , where(C1”): Equation (40) has six positive roots,(C2”): ( H − H ω ) (cid:0) H ω cos ωτ + ( H − ω ) sin ωτ (cid:1) +2 ωH ( H − H ω ) > H ω (cid:0) ( H − ω ) cos ωτ − H ω sin ωτ (cid:1) + H ω. Numerical Simulations
We begin by studying the original model (1) without the delay.
Time (in arbitrary units) C e ll P opu l a t i on ( × ) NTIADS (a) Cell Evolution with time
Time (in arbitrary units) C e ll P opu l a t i on β = 0 β = 0.02 β = 0.08 β =0.3 β =0.5 (b) Evolution of Tumor cells with various β Figure 1: Evolution of cells for (a) parameters given in (43) and (b) for various β Time (in arbitrary units) T u m o r C e ll P opu l a t i on α = 0 α = 0.75 α = 1 α =1.25 α =3.5 (a) Evolution of Tumor cells with various α δ α R (b) R for various α and δ Figure 2: Evolution of cells for (a) various α and (b) R for various α and δ For the numerical simulation of our model (1), we consider the parameter values as, r = 1 , r = 0 . , s = 0 . , b = 1 , b = 1 , c = 0 . , c ′ = 0 . , α = 0 . , α ′ = 0 . , β = 0 . , (43) β ′ = 0 . , γ = 0 . , d = 0 . , d D = 1 , d A = 0 . , η = 0 . , µ = 0 . , b = 0 . , ρ = 0 . , α = 0 . , (44) δ = 0 . , α = 0 . . (45)Corresponding to these parameter values, the equilibrium points are P = (0 , , . , , . , . ,P = (0 , . , . , , . , . P = (0 . , . , . , , . , . . For this set ofparameters, the reproduction number is R = 0 . . Here the four equilibrium points exist and tumor-free equilibrium does not exists. The simulation results for the model system (1) corresponding to theseparameter values and initial population (30000, 10000, 10000, 50000, 20000, 20000) are shown in Figure1(a). t can be observed from the time portrait that the system is asymptotically stable, with the solutionsconverging to the equilibrium point P . Furthermore, we observe the evolution of Tumor cells withvariations in parameters β and BT-immune threshold rate α . We can observe from Figure 1(b) thatas the competition term β increases, population of Tumor cells decreases. This is due to competitionbetween the Tumor cells and immune cells. Time(in arbitrary units) F r a c t i on o f P opu l a t i on Evolution of Tumor cells
Only ImmunityImmunity + Arsenic (I+A)Immunity + Arsenic + Drug (I+A+D)Constant InfluxI+A+D+Drug- Immune Term (a) Evolution of tumor cells with various interactionterms.
Time(in arbitrary units) F r a c t i on o f P opu l a t i on Evolution of Normal cells
Only ImmunityImmunity + Arsenic (I+A)Immunity + Arsenic + Drug (I+A+D)Constant InfluxI+A+D+Drug- Immune Term (b) Evolution of normal cells with various interactionterms.
Figure 3: Solution for variable-inlfux with BT-immune term converges to equilibrium state T*=0.3005 for tumors and N*= 0.6399 for normal cells
Time(in arbitrary units) F r a c t i on o f P opu l a t i on Evolution of Immunity cells
Only ImmunityImmunity + Arsenic (I+A)Immunity + Arsenic + Drug (I+A+D)Constant InfluxI+A+D+Drug- Immune Term (a) Evolution of immune cells with various interactionterms
Time(in arbitrary units) F r a c t i on o f P opu l a t i on Evolution of cells
Normal Cells(Variable flux)Tumor Cells(Variable flux)Normal Cells(Steady flux)Tumor Cells(Steady flux) (b) Evolution of cells under steady and variable influx.
Figure 4: Solution for steady immune converges to (N*, T*)= (0.0976,0.8431) and for variableinflux it converges to (N*, T*)= (0.6390, 0.3010)
The behaviours of tumor, normal and immune cells are studied using similar parameter set as before.In figure 3(a) and 3(b), the behaviours of the tumor and normal cells are illustrated in the presence ofvarious interaction components.As the system of tumor cells interacts with immune cells in steady influx model (black circles in fig.3(a)), the saturation value is high enough to consider the system dead. The population is simulataneously lowfor the normal cells as shown in fig. 3(b). This gives us the idea to supplement the immune cell growth tocounteract the rise of tumor cells. A variable influx is then assumed in system with immunity alone, andthe pronounced effect of variable influx is immediately observed as the tumor population comes down,while increasing the normal cells population (red dashed line in fig. 3(a) and 3(b)). Cancer inducing
50 100 150
Time(in arbitrary units) F r a c t i on o f P opu l a t i on N o r m a l c e ll s =0.1 =0.25 =0.5 Time(in arbitrary units) F r a c t i on o f P opu l a t i on T u m o r c e ll s =0.1 =0.25 =0.5 (a) Evolution of tumor cells with variation in µ Time(in arbitrary units) F r a c t i on o f P opu l a t i on N o r m a l c e ll s =0.1=0.25=0.5 Time(in arbitrary units) F r a c t i on o f P opu l a t i on T u m o r c e ll s =0.1=0.25=0.5 (b) Evolution of normal cells with variation in η Figure 5: Evolution of tumor and normal cells for various µ and η (a) Figure 6: phase portrait for model 1. Solution trajectories converge to its equilibrium state (N*,T*, A*)= (0.6399, 0.3005, 9.961) starting from different initial conditions. arsenic is then introduced in the system with the variable influx model and the tumor population israised, and decreasing the normal population. The is due to the fact that arsenic converts a fraction ofnormal cells into tumorous. Black tea is introduced in the system in form drug and while this does affectthe population of tumor and normal cells (green dash-dot in the figure), the introduction of BT-immuneinteraction has great impact in the system as shown in fig.3(a)(magenta dots).For variable influx model we observe that immune cells saturate at a higher value than for the steadyinflux and the effect is even pronounced if BT-immune term is incorporated. This indicates that it may bepossible to reach a healthy outcome with an immunomodulator and variable immune response, workingin tandem, without the intervention of chemotherapeutic drugs. Indeed, it depends on the values of theparameters and the severity of the disease, but what we are able to show here is that it may be possibleto dispense with chemotherapy or at least reduce its application to a large degree if suitable protocolswith BT and immune response could be achieved. hile the perfect cure for cancer is still not possible, we can aim to provide a good and prolonged life.One way to acheive that is by keeping the tumor cells always under control and lower than normal cellpopulation. This can be achieved by extensive chemotherapy but the price is paid in form of painful side-effects. Starting with a small tumor population, and supplementing the system with immunomodulatoryeffects of BT, we observe as shown in fig.4(b), that tumor population can be kept below a certain thresholdfor a longer amount of time.The variation in the immune cells in system(1)is brought about by s(t). These stimulated immunecells are assumed to have production rate of η and death rate of µ . Higher rate of production of immunecells has larger effect on tumor population as shown in fig.5(a). One has to be careful here as largeproduction rate can lead to immune cell proliferation and immune upon immune crowding. A similarscenerio is for death of stimulated immune cells (see fig.5(b)) and a smaller death rate can lead to similarproblems as that of larger production rate.Considering the parameteric values in sec. 9.1, the model system(1) has the equilibrium state solution(N*,T*,A*)= (0.6399, 0.3005, 9.961), and the solution curves stabilize to its equilibrium state(shown inphase portrait 6(a) ) In this section we present the results obtained for a delayed system (12). As discussed in Section 8, timedelay of τ and τ is incorporated in the normal and tumor cells respectively. This delay could be dueimmunity in the system .A delay τ = τ = τ has been added in Arsenic-Normal cell interaction. This kind of delay in thesystem can be analysed analytically as done in Section 8. Physically this delay implies that the arsenicin the system converts the normal cell into tumorous at a later time. The effect of the delay doesn’t havea large impact over the fraction of cells as shown in fig 7(a). For a delay of τ = 5 the normal cells farebetter very slightly before converging to its stable equilibrium point (see fig 7(b))A more physical delay in the form of delay in immune response is also studied. While the systemis complex, the analytic research about the characteristic equation can be carried out similarly, and wepresent the numerical analysis in this section to show the dynamical behaviour of the system with delayedimmune response. Time F r a c t i on o f c e ll s Normal Cells( =5)Tumor cells( =5)Normal Cells( =0)Tumor cells( =0) (a) Solution for model 12. Solution converges to its equi-librium state (N*, T*)= (0.6439, 0.3005)
Normal Cells T u m o r C e ll s =5=0 (b) Phase portrait for model 12. Solution converges toits equilibrium state (N*, T*)= (0.6439, 0.3005) Figure 7: Solution converges to its equilibrium state (N*, T*)= (0.6439, 0.3005)20
20 40 60 80 100 120 140
Time F r a c t i on o f c e ll s Normal Cells( =5)Tumor cells( =5)Normal Cells( =0)Tumor cells( =0) (a) Solution for model 46. Solution converges to its equi-librium state (N*, T*)= (0.6439, 0.3005)
Normal Cells T u m o r C e ll s =5=0 (b) Phase portrait for model 46. Solution converges toits equilibrium state (N*, T*)= (0.6439, 0.3005) Figure 8: Solution converges to its equilibrium state (N*, T*)= (0.6439, 0.3005) dNdt = r N (1 − b N ) − c T N − α ′ ANdTdt = r T (1 − b T ) − c ′ T N − βI ( t − τ ) T − γDT + αANdIdt = s ( t ) + ρI ( t − τ ) Tα + T − d I I − βI ( t − τ ) T + δI ( t − τ ) Dα + D (46) dAdt = a − d A AdDdt = b − d D Ddsdt = s + ηb + s s − µ s. A delay in tumor-immune and BT-immune response is studied in eqn.46. A delay of τ = 5 hasconsiderable effect on fraction of cells as shown in fig 8(a). The fraction of tumor cells rises above thenormal cells before converging towards equilibrium point (N*, T*)= (0.6439, 0.3005), as observed in Fig.8(b).
10 Conclusion
We have written down a model for the growth of cancer cells under the exposure of an environmentalcarcinogen and a mitigating agent. We deliberately exclude any chemotherapeutic treatment in the modelin order to test and establish a possible protocol without its use. The treatment of the model followstwo routes, one analytical looking for the local stability of tumor-free equilibrium, dead equilibrium andlocal and global stabilities of endemic equilibrium from the model.The Reproduction Number, which isan indicator of rate of spread has been computed by using next generation matrix method. The localstability of tumor-free equilibrium, dead equilibrium, endemic equilibrium and global stability of endemicequilibrium have been analyzed. Furthermore, stability and bifurcation analysis of delayed model havebeen discussed.A model system with compromised immune system has been studied with prime focus on usingimmunotherapy to obtain best results within certain general parameters. A comparison has been madebetween a constant influx model of immune cells studied earlier and the variable influx model of immunecells under consideration. Black tea, which has been studied extensively in literature is used as animmunomodulator to curb the side- effects of traditional chemotherapeutic treatments. Along with ariable influx model an important BT-immune interaction has been incorporated in the system. Constantinflux of immune cells saturates the tumor cells at unhealthy levels even under BT-immune interaction.The maximum effect on tumor cells comes from combining the two predominant factors, namely, variableinflux and BT-immune interaction. Starting from as high as 10 percent tumor cell count to begin with,the tumor cell population never exceeds the normal cells at any time during the evolution. The effect ofrate of production of immune cells and their death have been investigatyed and a higher production rate,and slower death rate of immune cells have greater effect on tumor cells; however this has to be exercisedwith caution as this may lead to immune proliferation and immune upon immune crowding.Finally, a physically relevant scenario with delay has been introduced. A delay in arsenic-normalcell interaction some impact on the general result, albeit small quantitative shifts in counts. On theother hand, a delay in immune cell response has considerable effect on the fractional count of cells. Adelay of τ = 5 results in higher peak value of tumor cells and even surpasses the fraction of normalcells highlighting the importance of rapid immune response in addition to the variable immune influx infighting cancer. References [1] Freddie Bray, Jacques Ferlay, Isabelle Soerjomataram, Rebecca L Siegel, Lindsey A Torre, andAhmedin Jemal. Global cancer statistics 2018: Globocan estimates of incidence and mortalityworldwide for 36 cancers in 185 countries.
CA: a cancer journal for clinicians , 68(6):394–424, 2018.[2] D Sinha, S Roy, and M Roy. Antioxidant potential of tea reduces arsenite induced oxidative stressin swiss albino mice.
Food and chemical toxicology , 48(4):1032–1039, 2010.[3] Sarkar Dibyendu and Rupali Datta. Biogeochemistry of arsenic in contaminated soils of superfundsites.
United States Environmental Protection Agency , 2007.[4] James Carelton. Final report: : Biogeochemistry of arsenic in contaminated soils of superfund sites.
United States Environmental Protection Agency , 2007.[5] DN Guha Mazumder. Diagnosis and treatment of chronic arsenic poisoning.
United Nations synthesisreport on arsenic in drinking water , 2000.[6] HM Srivastava, Urmimala Dey, Archismaan Ghosh, Jai Prakash Tripathi, Syed Abbas, A Taraphder,and Madhumita Roy. Growth of tumor due to arsenic and its mitigation by black tea in swiss albinomice.
Alexandria Engineering Journal , 2020.[7] Mohammad Rafiul Haque and Shahid Husain Ansari. Immunostimulatory effect of standardisedalcoholic extract of green tea (camellia sinensis l.) against cyclophosphamide-induced immunosup-pression in murine model.
International Journal of Green Pharmacy (IJGP) , 8(1), 2014.[8] Antony Gomes, Poulami Datta, Amrita Sarkar, Subir Chandra Dasgupta, and Aparna Gomes. Blacktea (camellia sinensis) extract as an immunomodulator against immunocompetent and immunodefi-cient experimental rodents.
Oriental Pharmacy and Experimental Medicine , 14(1):37–45, 2014.[9] Chandan Chattopadhyay, Nandini Chakrabarti, Mitali Chatterjee, Sonali Mukherjee, Kajari Sarkar,and A Roy Chaudhuri. Black tea (camellia sinensis) decoction shows immunomodulatory propertieson an experimental animal model and in human peripheral mononuclear cells.
Pharmacognosyresearch , 4(1):15, 2012.[10] Lisette G De Pillis and Ami Radunskaya. A mathematical tumor model with immune resistance anddrug therapy: an optimal control approach.
Computational and Mathematical Methods in Medicine ,3(2):79–100, 2001.[11] Lisette G De Pillis and Ami Radunskaya. The dynamics of an optimally controlled tumor model: Acase study.
Mathematical and computer modelling , 37(11):1221–1244, 2003.
12] John Carl Panetta. A mathematical model of periodically pulsed chemotherapy: tumor recurrenceand metastasis in a competitive environment.
Bulletin of mathematical Biology , 58(3):425–447, 1996.[13] S Michelson and JT Leith. Host response in tumor growth and progression.
Invasion & metastasis ,16(4-5):235–246, 1996.[14] Vladimir A Kuznetsov, Iliya A Makalkin, Mark A Taylor, and Alan S Perelson. Nonlinear dynamics ofimmunogenic tumors: parameter estimation and global bifurcation analysis.
Bulletin of mathematicalbiology , 56(2):295–321, 1994.[15] Carlos Castillo-Chavez, Zhilan Feng, and Wenzhang Huang. On the computation of ro and itsrole on.
Mathematical approaches for emerging and reemerging infectious diseases: an introduction ,1:229, 2002.[16] Pauline Van den Driessche and James Watmough. Reproduction numbers and sub-threshold endemicequilibria for compartmental models of disease transmission.
Mathematical biosciences , 180(1-2):29–48, 2002.[17] A Ghafari and N Naserifar. Mathematical modeling and lyapunov-based drug administration incancer chemotherapy.
Iranian Journal of Electrical and Electronic Engineering , 5(3):151–158, 2009.[18] Xia Yang, Lansun Chen, and Jufang Chen. Permanence and positive periodic solution for thesingle-species nonautonomous delay diffusive models.
Computers & Mathematics with Applications ,32(4):109–116, 1996.[19] Yongli Song and Sanling Yuan. Bifurcation analysis in a predator–prey system with time delay.
Nonlinear analysis: real world applications , 7(2):265–284, 2006., 7(2):265–284, 2006.