SQU Journal for Science, 2018, 23(1), 19-31 DOI: http://dx.doi.org/10.24200/squjs.vol23iss1pp19-31 Sultan Qaboos University 19 Numerical Simulations of a Delay Model for Immune System-Tumor Interaction Mohamed A. Hajji* and Qasem Al-Mdallal Department of Mathematical Sciences, College of Science, UAE University, P.O. Box 15551, Al Ain, UAE. *Email: mahajji@uaeu.ac.ae ABSTRACT: In this paper we consider a system of delay differential equations as a model for the dynamics of tumor - immune system interaction. We carry out a stability analysis of the proposed model. In particular, we show that the system can have up to two steady states: the tumor free steady state, which always exist, and the tumor persistent steady state, which exists only when the relative rate of increase of the tumor cells exceeds the ratio between the natural proliferation rate and the relative death rate of the effector cells. We also determine an upper bound for the delay, such that stability is preserved. Numerical simulations of the system under different parameter values are performed. Keywords: Delay differential equations, Asymptotic stability, Numerical simulations. باستخدام المعادالت التفاضلية التأخيريةلسرطاني ل نظام المناعة والورم االمحاكاة العددية لنموذج تفاع قاسم المدللو *حاجيعلي محمد ي. ونقوم نفترض في هذا البحث جملة من المعادالت التفاضلية التأخيرية كنموذج لديناميكيات نظام التفاعل بين األورام السرطانية والجهاز المناع :صالملخ نبين بشكل خاص امكانية احتواء النظام على حالتين ثابتتين: الحالة المستقرة الخالية من الورم، وهي موجودة دائًما، بتحليل استقرار النموذج المقترح. كما ي ومعدل الوفيات والحالة المستقرة للورم المستمر، والموجودة فقط عند تجاوز المعدل النسبي لزيادة الخاليا الورمية النسبية بين معدل االنتشار الطبيع .نعين أيضا الحد األعلى للتأخير مع المحافظة على االستقرار. وتم إجراء محاكاة عددية للنظام بتحديد قيم متنوعة للمعاملو للخاليا المستجيبه. النسبي .، محاكاة عدديةاالستقرار قاربمعادالت التفاضلية التأخيرية، ت :مفتاحيةالكلمات ال 1. Introduction he immune system (IS) is a very complex one. It is composed of a complex network of different cells which operate collectively by communicating information through signalling. The duty of these cells is to keep our system natural, i.e., free of foreign entities. When a foreign entity enters our system, particular immune system cells raise an alarm and send signals through the network to other immune system cells calling for an attack on the foreign entity. Once the foreign substance is located and identified as non-self, the immune system plans an attack strategy. Among the most deadly foreign entities are cancer cells. Cancer cells are difficult to deal with because of their ability to multiply so fast that the IS cannot keep up. Another property of cancer cells is that they are able to camouflage the antigens (the substance that triggers the IS) such that the IS cannot recognise them as foreign cells. Once cancer cells are identified as foreign, the IS starts an attack to destroy them. The interactions between tumour cells and the immune system are very complex and require sophisticated models to describe them. Mathematical models, based on ordinary differential equations, delay differential equations or partial differential equations, have proven to be useful tools in analysing and understanding the IS-tumor interactions. Several mathematical models have been suggested to describe the interactions between tumour and immune system [1-9]. T M.A. HAJJI and Q. AL-MDALLAL 20 Many of these papers consider the tumor-IS interaction with other factors such as treatment. In this paper, we consider a mathematical model of IS-tumor interactions based on a system of delay differential equations. The delay is introduced to reflect the non-instantaneous outcome of the interaction between the immune system cells and the tumor cells. We carry out a stability analysis of the system and we simulate the system under different parameter values to explore the different asymptotic behaviors. Many models that have appeared in the literature are based of the well-known Kuznetsov and Taylor model [1] which we describe in the next section. 2. Kuznetsov and Taylor model The Kuznetsov and Taylor mathematical model to describe the tumor-immune system interaction is based on the following: 1. The growth of tumor cells population T (in the absence of the immune system cells) follows a logistic model: (1)= (1 ), tot dT aT bT dt  where a is the maximum growth rate, 1 b is the carrying capacity of the biological environment for the tumor cells, and totT is the total population of unhit T cells, CTTtot = . 2. The rate of production of the immune system cells E (the effector cells), has two sources: (i) a constant normal production rate s (in the absence of the tumor cells) and (ii) a production rate caused by the presence of tumor cells. The magnitude of the second production rate is a function, ),( TCF , of the current concentrations of C and T , where C is the concentration of the complex TE  . The immune system cells E die out at a rate 1 d . 3. The interaction between the effector cells E and the tumor cells T forming the complexes TEC  is described by the following kinetic Figure 1. Interaction kinetics between E and T cells. where the C complex is formed at rate of 1 k and breaks at three different rates: 1 k giving E and T unaltered, 2 k given E (unaltered) and T D , and 3k giving ED and T (unaltered), where TD and ED are altered T and E cells which are bound to die. Kuznetsov and Taylor’s complete model is given by the following set of differential equations 1 1 1 2 1 1 3 1 1 2 3 3 2 2 3 (2) = ( , ) ( ) , = (1 ) ( ) , = ( ) , = , = , tot E E T T dE s F C T d E k ET k k C dt dT aT bT k ET k k C dt dC k ET k k k C dt dD k C d D dt dD k C d D dt                  where, again, E , T , C , E D , T D are the concentrations at the tumor site of the immune system cells (effector cells), the tumor cells, the effector-tumor cells complex, the inactived effector cells, and the inactived tumor cells, respectively. The above model was simplified based on the fact that the TE  bond, which is formed at a rate of 1 k , lasts a relatively short period of time [10] and disassociates giving rise to E and T (no win situation) with a rate 1 k , NUMERICAL SIMULATIONS OF A DELAY MODEL 21 or E D and T (tumor wins) with a rate 3k , or TD amd E ( E wins) with a rate 2k , see Figure 1. The inactivated effector cells E D and inactivated tumor cells T D are bound to die with rates 2 d and 3d , respectively. In [1], it was claimed that C is approximately constant, that is 0 dt dC , hence KETC  , with )/(= 3211 kkkkK   . In addition, totT was approximated by TTtot  , and ),(),( TEFTCF  . These assumptions reduce the above system (2) to the following IS-tumor interaction model, in non-dimensional form, (3) = ( , ) , = (1 ) , dE s F E T mET dE dt dT aT bT nET dt      where the parameters 3= Kkm , 2= Kkn , and 1= dd . The function ),( TEF describes the rate at which cytotoxic effector cells accumulate around the region where the T cells are localized. Based on Kuznetsov and Taylor’s model, many researchers have derived other models by considering different forms for the function ),( TEF and/or by adding additional considerations (see [3]-[8], and references therein). In [1], the authors analyzed model (3) with a Holling Type II form for F . In non-dimensionl variables, their model takes the form (4) = , = (1 ) . dx xy xy x dt y dy y y xy dt              Later, in [6], the authors considered Kuznetsov and Taylor’s model (3) with Holling Type I response function kETTEF =),( , with and without delay, where the delay was introduced to reflect the delay in the response of the immune system before proliferating effector cells. In this work, we consider (3) with a Holling Type I response function kETTEF =),( with delay introduced in both equations as described in the next section. 3. A delay model As in [6], we consider ),( TEF of Holling Type I, but we consider the delay  according to the following biologically possible scenario. We assume that there is delay in both the proliferation of effector cells and in the disposal of tumor cells. In dimensionless form, the model we consider is (5) ( ) = ( ) ( ) , ( ) = (1 ) ( ) ( ), x t wx t y t x y t y y x t y t                   where )(tx and )(ty are the nondimensional density at the tumor site of the effector cells and of the tumor cells, respectively. All parameters in (5) are nonnegative real numbers and are approximated (see [1] and [6]) by: 3 (6) = 0.1181, = 0.04, = 0.3743, = 1.636, = 2.0 10 . w      The delay factor in model (5) models the natural delays in both the proliferation of effector cells and in the disposal of tumor cells. This is biologically meaningful as the IS takes time to react and migrate effector cells to the tumor site. Also, in the second equation of model (5), it models the delay in completely eradicating the tumor cells as they will not die immediately upon interacting with effector cells. Our aim in studying the above model (5) is to investigate the effect, if any, of the time delay,  , in the disposal of the tumor cells. We carry out the stability analysis of the equilibrium states of the model (5) and determine necessary conditions for the local asymptotic stability. Also, we determine an upper bound for the delay parameter  such that stability is preserved. Numerical simulations of the model are performed to investigate the asymptotic behavior of the system. 4. Steady states The steady states, ),( yx , of (5) are the solutions of the nonlinear homogeneous system M.A. HAJJI and Q. AL-MDALLAL 22 (7) = 0, (1 ) = 0. wxy x y y xy         From the second equation of (7) we have either 0=y or )(1= yx   . If 0=y , then from the first equation of (7) we have /=x , giving the tumor free steady state ,0)/(=0 E . If 0y , then )(1= yx   and from the first equation of (7) we have the quadratic equation for y : 0=)( 2   ywyw whose discriminant  is 0.>4)(= 22 ww   Then depending on the parameter values, there may exist up to two more steady states ),( 111 yxE and ),( 222 yxE , where 1 1 (8) ( ) ( ) = , = 2 2 w w x y w w             2 2 (9) ( ) ( ) = , = 2 2 w w x y w w             The existence of 1 E and 2 E as biologically meaningful (i.e., with positive coordinates) depends on the parameters. If, as assumed, all parameters are positive, then it is clear that 0< 2 x because |)(>| w  . So 2 E does not exist as biologically meaningful. As for 1 E , since |)(>| w  , we have 0> 1 x . It remains to check when 0> 1 y . Calculations reveal that 0> 1 y if and only if  )/( , this biologically means that the effector cells will take over the tumor cells, which as a result, are bound to die in the long run. On the other hand, if  <)/( , the tumor cells will always exist. Now, we turn to the stability analysis of the steady states 0E and 1E . 5. Stability analysis We know that a steady state is locally asymptotically stable if all the eigenvalues of the Jacobian matrix of the system at the steady state have negative real parts. Before we continue, we recall an important stability criterion known the Mikhailov criterion [9]: Mikhailov criterion: Let ( )P s and ( )Q s be two polynomials with deg( P )= N and deg(Q ) < N . If the quasi-polynomial ( ) = ( ) ( ) s R s P s Q s e   has no roots on the imaginary axis, then all roots of ( )R s have negative real parts if and only if the argument of ( )R js , = 1j  , increases by / 2N  as s . For the system (5) under consideration, the Jacobian matrix at a steady state ),( yx is given by (10)= . 2 wye wxe J ye y xe                     The characteristic polynomial is (11)( ) = ( ) ( ) ,W P Q e       where 2 (11)( ) = ( ) , = 2 ,P A A A y          (12)( ) = ( ) ( ).Q x wyA x wy     NUMERICAL SIMULATIONS OF A DELAY MODEL 23 Stability of the tumor-free steady state 0E For the steady state ,0)/(=0 E , the characteristic polynomial reduces to (13)( ) = ( )( / ).W e             We have the following proposition. Proposition 1 (i) If  /> , then 0E is an unstable saddle point. (ii) If  /< , then there exists a 0>0 such that 0E is locally asymptotically stable for )[0, 0  and unstable for 0> where 1 0 2 2 2 (14) ( / )cos = . /         Proof. From (14), we see that one eigenvalue is 0<=   . Next, we need to analyze the roots of the quasi- polynomial (15)( ) = / .R e          It is clear that if  /> , )(R in (5.7) has a positive root, since the functions  =)(f and    eg /=)( satisfy  /=(0)>=(0) gf and )(f and 0)( g , as  . This proves (i). For (ii), assume that  /< . According the Mikhailov criterion, if )(R has no imaginary roots, then all its roots have negative real parts if and only if   /2.=))(( = 0=  s s jsRArg Let js= ,   Rs , 1= j . Then (16) becomes (16)( ) = [ / cos( )] [ / sin( )].R js s j s s          Let ))((= jsRArg . Then as s , we have and (17)cos( ) 0 sin( ) 1.   This implies that /2=))((lim jsRArg s  . We have 0>/=(0)  R (for  /< ) and 0=(0))(RArg . Thus   /2,=))(( = 0=  s s jsRArg which implies, according the Mikhailov criterion, that all roots of )(R have negative real parts, provided that )(R has no imaginary roots. Now let us see under what condition )(R has imaginary roots. If 0=)( jsR for some   Rs , then (18)/ cos( ) = 0,s    (19)/ sin( ) = 0.s s   Squaring and adding (19) and (20) gives ./=/= 2222222   ss Conversely, if 222 0 /==  ss , the first value of  , say 0 , that satisfies both (19) and (20) is 1 1 0 2 2 2 0 (20) 1 ( / )cos = ( / ) = .cos /s            This implies that for 0< , )(R has no imaginary roots, and 0E is locally asymptotically stable. This completes the proof of Proposition 1. □ Stability of the steady state ),(= 111 yxE For ),(= 111 yxE , from equations (11)-(13), the characteristic polynomial takes the form (21)( ) = ( ) ( ) ,W P Q e       M.A. HAJJI and Q. AL-MDALLAL 24 where 2 0 1 1 (22)( ) = ( ) , ( ) = ( ) ( ),P A A Q C A x wy             where .=,2= 1101 AwyxACyA   Before we proceed, we state and prove the following Lemma about the positivity of 0C , which will be used later. Lemma 2 0>== 1110  yAwyxAC  . Proof. From the expression of   1 2= yA and )(1= 11 yx   , we can write A as 11 = xyA  . Substituting this into 0C , we get )].([= 1110 xywyC   Now, from the expressions of 1 x and 1 y (Eq. (8)), we have .= 2 )( 2 )( = 11 ww w w w xy       It follows that 0.>=)]([= 11110  yxywyC  We have the following proposition for the stability of 1 E . □ Proposition 3. Assume that AC 2>0 . The steady state ),(= 111 yxE is locally asymptotically stable, for )[0, 0   , where 0 is given by (24). Proof. Let   Rsjs,= . We have (23)( ) = ( ) ( ) = ( ) ( ) js W js P js Q js e X s jY s    where 2 1 1 1 1 (24)( ) = ( ) ( ) cos( ) ( ) sin( ),X s s A x wy A s s x wy s         1 1 1 1 (25)( ) = ( ) ( ) cos( ) ( ) sin( ).Y s s A s x wy s x wy A s        If ))((= jsWArg , it can be easily verified that (26) ( ) ( ) cos( ) = 1, sin( ) = 0. | ( ) | | ( ) |w w X s Y s W j W j         This implies that =))((lim jsWArg s  . For 0=s , we have 1 1 0 (27)(0) = = > 0,W x wy A A C   which implies that 0=(0))(WArg . Hence,   =))(( =0= s s jsWArg , and 1 E would be locally asymptotically stable, provided that )(W has no roots on the imaginary axis. Now, we check for pure imaginary roots of )(W . If there exists   Rs such that 0=)( jsW , then 0=)(=)( sYsX , or 2 11 12 (28)cos( ) sin( ) = ( ),a s sa s s A    12 11 (29)cos( ) sin( ) = ( ),sa s a s s A     where .=,= 11121111 wyxaAwyxa  Squaring both sides of (29) and (30) and adding, we get 4 2 1 2 (30)= 0,s a s a  where )2(=)(=,= 00 2 11 2 2 2 12 22 1 ACCaAaaAa   The solutions of (31) are formally written as NUMERICAL SIMULATIONS OF A DELAY MODEL 25 2 1 1 22 (31) 4 = . 2 a a a s    It is clear that a sufficient condition for the right hand side of (32) to be positive is 0< 2 a , i.e., AC 2>0 , since 0> 0 C . In this case, the positive solution of (32) is 2 1 1 2 0 (32) 4 = 2 a a a s    Therefore, if AC 2>0 , then )(W has only two pure imaginary roots 0js , where 0s is given by (33). With 0s as in (33), the first positive value of  such that 0=)(sX and 0=)(sY is 2 1 11 12 0 11 0 2 2 2 0 11 0 12 (33) ( ( ))1 = .cos a a A s a A s a s a             Therefore, if AC 2>0 , then for 0< , )(W has no imaginary roots, and 1E is locally asymptotically stable. □ 6. Numerical simulations In this section, we perform a number of numerical simulations to reveal the dynamics of system (5). We consider the two cases (i)  /< where only 0E exists, and (ii)  /> where both 0E and 1E exist. Case 1: For the case  /< , we simulate the system with the values of 0.1= and 0.3= . The other parameters are fixed as in (6). In this case, the only equilibrium state is 0 = ( / , 0) = (0.315522, 0)E   which was proved to be asymptotically stable up to an upper bound for the delay  (see Proposition 1 and Eq. (21)). The calculated upper bounds according to (21) are 0 = 4.17134 and 0 = 3.22237 for = 0.1 and 0.3 , respectively. We remark that these upper bounds are valid only when the solutions ( )x t and ( )y t remain non-negative. In fact, simulations reveal that the solutions become negative for much smaller values of  , 0.05  . The outcomes of the simulations for this case are displayed in Figures 2, 4, and 5. In Figure 2, we display the phase portraits of system (5) for = 0.1 and = 0.3 and = 0 and = 0.05 , each with 5 differential initial conditions. They all confirm the stability of 0E , as the solutions )(tx and )(ty asymptotically tend to 0E . In Figures 4 and 5, the solutions of system (5) are displayed with = 0.1 and = 0.3 , respectively, each with = 0.05 . Both figures confirm the convergence of the solutions to the steady state 0 E . It is important to mention that for > 0.05 , under certain initial conditions, the solutions become negative at some time * t at which the system becomes biologically invalid. Case 2: For the case  /> , we fix all parameters as in (6) and simulate the system with different delay values  . In this case, we have two steady states, 0 = ( / , 0) = (0.315522, 0)E   which is unstable and .52522)(1.61138,7=),(= 111 yxE which was proven to be locally asymptotically stable for up to 0< . From Eq. (34), 0 is calculated to be 0 = 0.0880838 . To confirm these results, we simulated the system for = 0, 0.01, 0.05, 0.08, 0.09 and various initial conditions. The results are displayed in Figures 6 and 7, where Figure 6 displays the phase portraits for = 0, 0.01, 0.05, 0.08, 0.09 and Figure 7 displays the solutions of the system for = 0.05, 0.08, 0.09 and initial condition (0) = 1, (0) = 30x y . It is to be noted that for 0 = 0, 0.01, 0.05, 0.08 <  , the solutions converge to 1 E , confirming the local asymptotic stability of 1 E . However, in the case 0 = 0.9 >  , we notice the formation of the closed orbit around 1 E . This also can be seen in the right- most plot in Figure 7, where the effector and tumor cells oscillate around 1 E . This confirms the existence of a bifurcation with respect to the parameter  as it crosses 0 = 0.0880838 . It should be also mentioned that for certain initial conditions the solutions of the system become negative in finite time. M.A. HAJJI and Q. AL-MDALLAL 26 Figure 2. Phase portraits of system (5) with different initial conditions ( (0), (0)) = (0,1), (1, 5), (2,10), (4,15)x y and (5, 20) : (a) = 0.1 and 0= , (b) = 0.1 and = 0.05. The equilibrium point is 0 = (0.315522, 0)E . Figure 3. Phase portraits of system (5) with different initial conditions ( (0), (0)) = (0,1), (1, 5), (2,10), (4,15)x y and (5, 20) : (a) 0.3= , 0= , (b) 0.3= and 0.05.= The equilibrium point is 0 = (0.315522, 0)E . a b a b NUMERICAL SIMULATIONS OF A DELAY MODEL 27 Figure 4. Solutions of system (5) with = 0.1 and = 0.05 and different initial conditions: (a) (0) = 0, (0) 1x y  ; (b) (0) = 1, (0) 5x y  ; (c) (0) = 2, (0) 10x y  ; (d) (0) = 4, (0) 15.x y  a b c d M.A. HAJJI and Q. AL-MDALLAL 28 Figure 5. Solutions of system (5) with = 0.3 and = 0.05 and different initial conditions: (a) (0) = 0, (0) 1x y  ; (b) (0) = 1, (0) 5x y  ; (c) (0) = 2, (0) 10x y  ; (d) (0) = 4, (0) 15.x y  a b c d NUMERICAL SIMULATIONS OF A DELAY MODEL 29 Figure 6. Solutions of system (5) with different initial conditions: (a) = / 0.1636    and = 0 ; (b) = / 0.1636    and = 0.05; (c) = / 0.1636    and = 0.08; (d) = / 0.1636    and = 0.09. a b c d M.A. HAJJI and Q. AL-MDALLAL 30 Figure 7. Solutions of system (5) with initial conditions: (0) = 1, (0) = 30x y : (a) = / 0.1636    and = 0.05; (b) = / 0.1636    and = 0.08; (c) = / 0.1636    and = 0.09. 7. Conclusion A mathematical model for the interaction between the immune system cells and tumor cells was proposed and studied. The model was derived from Kuznetsov and Taylor’s model [1] by introducing a delay term which reflects the slow response of the immune system cells to the presence of tumor cells. The importance in studying this delay model is in determining the effect of the delay term on the dynamics of the system. A stability analysis of the equilibrium states was carried out. Numerical simulations of the system were done for a number of different initial conditions. The simulation results confirm the theoretical stability results. Numerically, it was also observed that as the delay term  crosses the upper bound 0 , the steady state 1E loses its stability, resulting in the formation of closed orbits around 1 .E This is an important result in the sense that as the IS gets weaker (slower response, i.e., larger delay), the size of the tumor cells persists in an oscillating fashion (closed orbits). References 1. Kuznetsov, V.A., Makalkin, I.A., Taylor, M. and Perelson, A.S. Nonlinear dynamics of immunogenic tumors: parameter estimation and global bifurcation analysis. Bulletin of Matematical Biology, 1994, 56, 295-321. a b c NUMERICAL SIMULATIONS OF A DELAY MODEL 31 2. Kuznetsov, V.A. and Volkenshtein, M.V. Dynamics of cellular immunological anti-tumor reactions II. Qualitative analysis of the model (in Russian), Mathematical Methods of Systems Theory, Frunze: Kirghiz State University, 9191,1, 72-100. 3. Rihan, F.A., Abdel Rahman, D.H., Lakshmanan, S. and Alkhajeh, A.S. A time delay model of tumour-immune system interactions: Global dynamics, parameter estimation, sensitivity analysis. Applied Mathematics and Computation, 2014, 232, 606-623. 4. Rihan, F.A., Hashish, A., Al-Maskari, F., Hussein, M.S., Ahmed, E., Riaz, M.B. and Yafia, R. Dynamics of tumor- immune system with fractional-order. Journal of Tumor Research, 2016, 2(1),109-115. 5. Rihan, F.A. and Rihan, N.F. Dynamics of cancer immune system with external treatment and optimal control. Journal of Cancer Science and Therapy, 2016, 8(10), 257-261 . 6. Galach, M. Dynamics of the tumor-immune system competition-the effect of time delay. International Journal of Applied Mathematics and Computer Science, 2003, 13(3), 395-406. 7. Arciero, J.C., Jackson, T.L. and Kirschner, D.E. Mathematical model of tumor-immune evasion and siRNA treatment. Discrete and Continuous Dynamical Systems - Series B, 2004, 4(1),39-58. 8. Kirschner, D. and Panetta, J.C. Modeling immunotherapy of the tumor-immune interaction. Journal of Mathematical Biology, 1998, 37, 235-252. 9. Y. Kuang, 1993, Delay Differential Equations with Applications in Population Dynamics, London, Academic Press. 10. Dionysiou, D.D. and Stamatakos, G.S. Applying a 4D multiscale in vivo tumor growth model to the exploration of radiotherapy scheduling: the effects of weekend treatment gaps and p53 gene status on the response of fast growing solid tumors. Cancer Inform, 2006, 2,113-121. Received 15 June 2017 Accepted 17 December 2017