Exploring dynamical complexity in a time-

delayed tumor-immune model

Cite as: Chaos 30, 123118 (2020); https://doi.org/10.1063/5.0025510
Submitted: 17 August 2020 . Accepted: 14 November 2020 . Published Online: 04 December 2020

Parthasakha Das, Ranjit Kumar Upadhyay, Pritha Das, and Dibakar Ghosh

ARTICLES YOU MAY BE INTERESTED IN

Analysis of tumor-immune dynamics in a delayed dendritic cell therapy model

Chaos: An Interdisciplinary Journal of Nonlinear Science 30, 113108 (2020); https://
doi.org/10.1063/5.0006567

Loewner driving force of the interface in the 2-dimensional Ising system as a chaotic
dynamical system
Chaos: An Interdisciplinary Journal of Nonlinear Science 30, 113130 (2020); https://
doi.org/10.1063/5.0023261

Explosive synchronization in temporal networks: A comparative study

Chaos: An Interdisciplinary Journal of Nonlinear Science 30, 113135 (2020); https://
doi.org/10.1063/5.0023329

Chaos 30, 123118 (2020); https://doi.org/10.1063/5.0025510 30, 123118

© 2020 Author(s).
 Chaos ARTICLE scitation.org/journal/cha

Exploring dynamical complexity in a

time-delayed tumor-immune model
Cite as: Chaos 30, 123118 (2020); doi: 10.1063/5.0025510
Submitted: 17 August 2020 · Accepted: 14 November 2020 ·
Published Online: 4 December 2020 View Online Export Citation CrossMark

Parthasakha Das,1,a) Ranjit Kumar Upadhyay,2 Pritha Das,1 and Dibakar Ghosh3

AFFILIATIONS
1
Department of Mathematics, Indian Institute of Engineering Science and Technology, Shibpur, Howrah 711103, India
2
Department of Mathematics and Computing, Indian Institute of Technology (Indian School of Mines), Dhanbad 826004, India
3
Physics and Applied Mathematics Unit, Indian Statistical Institute, Kolkata 700108, India

a)
Author to whom correspondence should be addressed: parthasakha87das@gmail.com

ABSTRACT
The analysis of dynamical complexity in nonlinear phenomena is an effective tool to quantify the level of their structural disorder. In particu-
lar, a mathematical model of tumor-immune interactions can provide insight into cancer biology. Here, we present and explore the aspects of
dynamical complexity, exhibited by a time-delayed tumor-immune model that describes the proliferation and survival of tumor cells under
immune surveillance, governed by activated immune-effector cells, host cells, and concentrated interleukin-2. We show that by employing
bifurcation analyses in different parametric regimes and the 0–1 test for chaoticity, the onset of chaos in the system can be predicted and
also manifested by the emergence of multi-periodicity. This is further verified by studying one- and two-parameter bifurcation diagrams for
different dynamical regimes of the system. Furthermore, we quantify the asymptotic behavior of the system by means of weighted recurrence
entropy. This helps us to identify a resemblance between its dynamics and emergence of complexity. We find that the complexity in the model
might indicate the phenomena of long-term cancer relapse, which provides evidence that incorporating time-delay in the effect of interleukin
in the tumor model enhances remarkably the dynamical complexity of the tumor-immune interplay.
Published under license by AIP Publishing. https://doi.org/10.1063/5.0025510

Exploring the mechanisms in tumor-immune dynamics is of in the tumor model enriches the dynamical complexity of the
utmost importance in oncology. Typically, anti- and pro-tumor tumor-immune interplay.
factors cause complex variations in the tumor-immune interplay.
Here, we study a time-delay, tumor-immune model that con-
sists of tumor cells, activated immune-effector cells, host cells,
and concentrated interleukin-2. The time-delay in the model is I. INTRODUCTION
introduced due to the delay in the stimulation and transporta- Complexity science has significantly contributed to improving
tion of immune cells. We discuss the qualitative properties of the and saving human lives over the last few decades.1,2 Indeed, under-
model and perform one- and two-parameter bifurcation analyses standing the complex mechanism of tumor-immune (TI) dynamics
to identify regimes of regular and chaotic behavior. We verify the is a challenging task for bio-mathematicians, bio-physicists, and bio-
appearance of chaos by employing the 0–1 test and augment the medical engineers3 alike. Modeling of TI dynamics is an interesting
study of the asymptotic behavior with the analysis of its behavior research topic on its own.4 The proliferation of tumor cells depends
in appropriate parametric spaces. This is further corroborated by on immune-effector cells, i.e., on natural killer cells, macrophages,
computing Poincaré maps and by means of weighted recurrence activated cytotoxic-T lymphocytes, tumor and host cells, cytokines
entropy calculations, both allowing for the study of dynamical (IL-2, IL-6, and IL-10), and the micro-environment of the tumor.5–10
complexity in the system. Our analysis shows the appearance of Part of the proliferation process is the slow-fast stage. During that
regular and chaotic dynamics and the emergence of complex- stage, routing imaging is essential for monitoring the recurrence
ity. We find that the complexity in the model might indicate of cancer, the so-called “tumor dormancy.”11 The immune system
the phenomena of long-term cancer relapse, which provides evi- of a patient with tumor cells shows rather irregular and unpre-
dence that incorporating time-delay in the effect of interleukin dictable behavior due to complex inter-cellular interactions in the

Chaos 30, 123118 (2020); doi: 10.1063/5.0025510 30, 123118-1

Published under license by AIP Publishing.
 Chaos ARTICLE scitation.org/journal/cha

tumor.5,12–14 As time-delays affect the inter-cellular interactions and that larger delay leads to stronger tumor proliferation and to long-
signaling, it is necessary to consider the delayed response when term cancer relapse. They also suggested that the growth rate of
modeling tumor-immune interactions.15–17 Due to the diversity and immune-effector cells and the strength of host cells are necessary to
non-linearity of cancer-immune interactions, different attractors develop in their competition against tumor cells. However, the effec-
(e.g., point attractors, strange attractors, limit cycles, etc.) have been tor cells are among the most significant cells in the human immune
observed.16,18 In fact, the existence of strange attractors was found to system as they fight viruses, foreign micro-organisms, and tumor
govern the chaotic behavior in dynamical systems at large.19–21 It is cells. Natural killer cells, cytotoxic T-lymphocytes and cytokines,
worth noting that, in the context of cancer dynamics, chaos has been such as IL-2 and IL-6, are also involved in the tumor-immune inter-
studied in the time-delay TI system.22,23 play, that is in the interplay between tumor and immune cells.
Chaos is attributed to sensitive dependence of the system on Effector cells annihilate tumor cells by delivering biochemical sig-
initial conditions. In the framework of TI dynamics, this sensitiv- nals. Moreover, the identification and annihilation of tumor cells
ity gives rise to local instability and invokes the proliferation of by effector cells are not instantaneous processes. Cytokines play
cancer-specific patterns.24–26 This translates to different develop- an important role in the activation of immune cells. Cytokine-
ment of tumors in patients. Consequently, understanding better TI mediated cell communication and orchestrated immune response
dynamics is a challenging task for oncologists and clinicians but are time dependent and, thus, necessitate the use of time-delay when
attractive to researchers working on the mathematical modeling of modeling cancer dynamics at the cellular level.
cancer dynamics. The study of the emergence of chaos in cancer Here, we consider the time-delay TI model proposed by
dynamics can potentially illustrate advanced knowledge in cellular Das et al.,38
interactive processes. In this regard, Lyapunov exponents quantify γ1 ET
the chaotic behavior in a dynamical system by providing the rate Ṫ = α1 T(1 − β1 T) − − r1 HT,
of divergence of trajectories that emanate from infinitesimally close k1 + T2
initial conditions. The authors in Refs. 27 and 28 proposed the 0–1 Ḣ = α2 H(1 − β2 H) − r2 TH,
test to characterize the dynamics in a system by utilizing the asymp- (1)
γ2 ET
totic growth of state-space trajectories. Complementary to these I˙L = − δ2 IL ,
approaches, Poincaré maps prove useful in visualizing the qualita- k2 + T2
tive properties of the dynamics in a system.29 Interestingly, as we γ3 E(t − τ )IL (t − τ )
study here for the time-delay TI system (1), the presence of chaos in Ė = aT − δ3 E + ,
k3 + I2L (t − τ )
the system gives rise to complex dynamics and to entropy increase.
Dynamical complexity is associated with structural disor- with initial conditions
der in a system and can be quantified by means of entropy η(ν) = (η1 (ν), η2 (ν), η3 (ν), η4 (ν)), (2)
computations.30,31 For example, the weighted entropy is a measure
of uncertainty provided by probabilistic experiments on qualitative and η(ν) ∈ C([−τ , 0], R4+ ) such a way that T(ν)
weights of possible events.32 Weights are computed by the inverse = η1 (ν), H(ν) = η2 (ν), IL (ν) = η3 (ν), and E(ν) = η4 (ν), where
exponential of a distance metric, indicating the dispersion of points ν = [−τ , 0]. Here, C is a set of continuous maps on Banach
in the phase space of the system. Consequently, weighted entropy space in [−τ , 0] → R4+ with norm defined by kηk = sup[−τ ,0] {kηi k;
measures the complexity in a system,33 attributed to chaotic behav- i = 1, 2, . . . , 4} and τ ∈ [0, ∞). For biologically, practility, t is
ior. Moreover, entropic quantities based on recurrence plots have assumed that ηi (0) > 0; i = 1, 2, . . . , 4.
also been employed to quantify the recurrence properties of chaotic In system (1), T denotes the number of tumor cells in time, H
dynamics, which is an effective approach to determine the diver- denotes the number of host or normal cells in time, IL denotes the
gence of trajectories.34 Recurrence plots have been used to study change in cytokines in time, and E denotes the number of cancer-
heart beats, neural signals, biophysical systems, and ecosystems,35,36 specific effector cells in time. Also, τ is the time-delay introduced
and weighted recurrence plots and entropic measures have been to the model to account for delayed interactions in the activa-
proposed, for example, in Ref. 37. tion of immune cells, cytokine-mediated cell communication, and
This paper is organized as follows: In Sec. II, we discuss the orchestrated immune response.
time-delay TI model used in our study. The qualitative dynamics The first equation in system (1) describes the growth dynam-
is investigated in Sec. III by performing bifurcation analyses. In ics of the number of tumor cells. In particular, the number of
Sec. IV, we focus on the dynamical complexity in the TI system and tumor cells, T, increases logistically with proliferation rate α1 > 0
study the long-term behavior of its dynamics and weighted entropy and reciprocal carrying capacity β1 > 0. Immune-effector cells anni-
in recurrence plots to measure the disorder produced in regimes hilate tumor cells with rate γ1 > 0, governed by Monod–Haldane
related to cancer dynamics. Finally, in Sec. V, we discuss our work kinetics. The second equation describes the dynamics of the number
and compare it with the work in the field. of host or normal cells, H, which follows a logistic growth with rate
α2 > 0 and reciprocal carrying capacity β2 > 0. The competition
between tumor and host cells is described by the last two terms in the
II. THE TIME-DELAY TI MODEL first two equations related to rates r1 > 0 and r2 > 0, respectively.
In the last few decades, various mathematical models have The third equation describes the rate of change in cytokines IL ,
been proposed to model cancer dynamics.3,4 Ghosh et al.22 and and γ2 > 0 indicates the recruitment rate in the presence of tumor
Khajanchi et al.23 studied a time-delay tumor model and observed cells by self-limiting production of IL . Parameter δ2 > 0 denotes the

Chaos 30, 123118 (2020); doi: 10.1063/5.0025510 30, 123118-2

Published under license by AIP Publishing.
 Chaos ARTICLE scitation.org/journal/cha

Next, we perform a detailed numerical analysis of system (1)

for various combinations of the system parameters and initial con-
ditions, which will show that the solution remains positive and
bounded.

III. QUALITATIVE ANALYSIS

We start our qualitative analysis by studying the dynamics of
system (1) and investigating its dynamical properties. In particular,
we begin by analyzing the dynamical behavior in the neighborhood
of its fixed points.

A. Linear stability analysis of fixed points

System (1) admits the following four fixed points:
• The trivial or no-living cell fixed point E(0) = (0, 0, 0,0).
FIG. 1. Schematic illustration of the cellular interactions among tumor cells T, 
host cells, H, concentrated interleukin-2, IL , and effector cells, E, in the time-delay • The axial or only host-cell-present fixed point E(1) = 0, β12 , 0, 0 .
TI system (1). Here, the arrows show the direction of influence among the different  
types of cells, and δ2 > 0 and δ3 > 0 denote the decay rate of interleukin-2, IL , • The host-cell-free fixed point E(2) = T̂, 0, ÎL , Ê , where
and the decay rate of cancer-specific effector cells, E, respectively. Also, τ is the    
α T̂ 1−β1 T̂ k1 +T̂2

time-delay that accounts for delayed interactions in the activation of immune cells, 1 γ3 ÊÎL γ2 ÊÎL
T̂ = a δ2 E − k +Î2 , ÎL = δ k +Î2 , and Ê = .
cytokine-mediated cell communication, and orchestrated immune response. 3 L ( 2 L)  δ1 T̂

(∗)
• The cancer-present fixed point E = Ť, Ȟ, ǏL , Ě , where
 
Ť = a1 δ2 E − kγ3+ĚǏǏL2 , Ȟ = α2α−r 2 Ť
, ǏL = δ γk2 Ě+ǏǏL2 , and Ě
3 L 2 β2 ( 2 L)
 k1 +Ť2 

decay rate of IL . The last equation corresponds to the number of
  
α2 −r2 Ť
= α Ť 1 − β1 Ť − Ť α2 β2 .
cancer-specific effector cells, E, changing in time, where tumor anti- 1δ Ť

genicity, a > 0, increases the activation rate of effector cells at a rate To study the stability of the fixed points, we employ the Jacobian
of γ3 > 0. Here, δ3 > 0 indicates the decay rate of the effector cells. matrix J of system (1)
Moreover, IL activates effector cells by a paracrine and autocrine
process with a delayed interaction, τ . Figure 1 shows schematically g11 −r1 T 0 g14
 
the cellular interactions in the TI system (1), and Table I summarizes −r H g22 0 0
J= 2 ,

the values and parameter ranges used. g31 0 −δ2 g34
a 0 g43 e−λτ −δ3 + g44 e−λτ
γ1 E(k1 −T2 ) 1T
where g11 = α1 (1 − 2β1 T) − − r1 H, g14 = − k γ+T 2 , g22
TABLE I. The values and ranges of the parameters used in the time-delay (k1 +T2 )2 1

TI model (1). γ2 E(k2 −T2 ) 2T

= α2 (1 − 2β2 H) − r2 T, g31 = , g34 = k γ+T 2,
(k2 +T2 )2 2
γ3 E(k3 −I2L )
Parameter Description Unit Value g43 = , and g44 = kγ33+IILL .
(k3 +I2L )2
α1 Growth rate of T day−1 1 Next, we compute the eigenvalues of J at four fixed points in
β1 (Carrying capacity of T)−1 cell−1 0.1 the absence of time-delay, i.e., for τ = 0, considering the instanta-
γ1 Decay rate of T day−1 5.56 neous interaction between IL and E, as it is an easier case to treat
k1 Half-saturation constant volume 10 semi-analytically. We will get back to the study of the dynamics of
r1 Fractional killing rate of T by H cell−1 day−1 [0.001, 0.2] system (1) for τ > 0 in Subsection. III B. To calculate the eigenval-
α2 Growth rate of H day−1 1 ues for τ = 0, we used the parameter values in Table I, keeping r1
β2 (Carrying capacity of H)−1 cell−1 day−1 1 and δ2 fixed at 0.05 and 0.65, respectively. These values are within
r2 Fractional killing rate of H by T cell−1 day−1 0.55 the ranges for r1 and δ2 reported in the table.
γ2 Recruitment rate of IL day−1 27.77 The eigenvalues at the fixed point E(0) are λ(0) 1 = α1 > 0, λ2
(0)

k2 Half-saturation constant volume 10 (0) (0)

= α2 > 0, λ3 = −δ2 < 0, and λ4 = −δ3 < 0, as α1 , α2 , δ2 , and δ3
δ2 Decay rate of IL day−1 [0.2, 0.7] are all positive. Therefore, there are stable and unstable manifolds
a Antigenicity of T day−1 0.2 emanating from E(0) , and thus, it is a saddle fixed point.
δ3 Decay rate of E day−1 0.167 The eigenvalues at the fixed point E(1) are λ(1) r1
1 = α1 − β 2 ,
γ3 Recruitment rate of E cell−1 day−1 1.1
k3 Half-saturation constant volume 10 λ(1) (1) (1)
2 = 0, λ3 = −δ2 < 0, and λ4 = −δ3 < 0. Thus, system (1)
τ Time delay day (0, 25] undergoes a transcritical bifurcation at the critical point α1 = βr12 as
r1 ∈ [0.001, 0.2] (see Table I).

Chaos 30, 123118 (2020); doi: 10.1063/5.0025510 30, 123118-3

Published under license by AIP Publishing.
 Chaos ARTICLE scitation.org/journal/cha

The eigenvalues at the fixed point E(2) are λ(2) 1 ≈ 0.894 648,
λ(2)
2 ≈ −0.719 801, λ(2) 3 ≈ −0.008 328 68 − 0.267 752i, and λ4
(2)

≈ −0.008 328 68 + 0.267 752i. Consequently, E(2) is a saddle-focus

fixed point and, thus, unstable.
The linear stability analysis at the cancer present fixed point
E(∗) in the absence of delay shows that it is asymptotically stable
under certain restrictions, as discussed in Das et al.38
So far, we have investigated the dynamics of system (1) in the
FIG. 3. Two-parameter bifurcation diagrams by varying δ2 and r1 in (a), τ and r1
absence of time-delay, i.e., for τ = 0. Since our aim is to explore the
in (b), and τ and δ2 in (c). The color bar denotes the regions of periodic behavior
time-delay dynamics of system (1), in the following, we will study its of increasing orders. Here, we have used the parameter values in Table I, keeping
properties for τ > 0. τ = 23 in panel (a), δ2 = 0.65 in panel (b), and r1 = 0.05 in panel (c).

B. Bifurcation analysis
periodic dynamics of different periods can be seen as r1 increases
To do so, we compute one- and two-parameter bifurcation dia- in (0, 0.2]. On the other hand, the increase in the decay rate δ2 in
grams for different settings. We start by keeping all parameters as in [0.2, 0.7] of interleukin-2 promotes instability and periodic behav-
Table I and vary the interaction rate between tumor and host cells, ior in the system as seen in Fig. 2(b). Similarly, a period-doubling
r1 in (0, 0.2] in Fig. 2(a), the decay rate of IL cells, δ2 in [0.2, 0.7] in route to chaos can be seen in Fig. 2(c) with the increase in the time-
panel (b), and the time-delay of IL activation of the effector cells, τ delay τ ∈ (0, 25]. This shows that large time-delay is not enough to
in (0, 25] in panel (c). In Fig. 3, we plot two-parameter bifurcation instigate effector cells against tumor cells. In all cases, system (1)
diagrams for r1 and δ2 in panel (a), for τ and r1 in panel (b), and for exhibits multi-periodic states and chaotic behavior as can be seen
δ2 and τ in panel (c), keeping τ = 23 in panel (a), δ2 = 0.65 in panel in the bifurcation diagrams in Fig. 2, depending on the values of the
(b), and r1 = 0.05 in panel (c). bifurcation parameters.
We compute the local maxima of the tumor variable T, Tmax , for Furthermore, to study the effects of r1 , δ2 , and τ on the qual-
r1 ∈ (0, 0.2], δ2 ∈ [0.2, 0.7], and τ ∈ (0, 25] in panels (a), (b), and (c) itative behavior of system (1), we plot in Fig. 3 two-parameter
in Fig. 2, respectively. Figure 2(a) shows the one-parameter bifurca- bifurcation diagrams of Tmax by varying δ2 and r1 in [0.2, 0.7] and
tion diagram varying r1 for fixed δ2 = 0.65 and τ = 23. Chaotic and (0, 0.5] [panel (a)], τ and r1 in (0, 25] and (0, 0.5] [panel (b)], and τ
and δ2 in [0, 25] and (0.2, 0.7] [panel (c)]. The color bar represents
different orders of periodic behavior. The results in Fig. 3 signify the
existence of higher-order periodicity in system (1). Consequently,
higher-order periodic behavior leads to period doubling cascade to
chaos as these parameters vary. These results show rich dynamical
behavior that we explore next.

C. Dynamical behavior analysis employing the 0–1 test

Here, we complement the qualitative analysis in Sec. III by
studying the chaotic behavior of the time-delay TI model (1) and
employing the 0–1 test,27,28 a method that can discriminate effec-
tively between regular and chaotic dynamics.
According to the 0–1 test,27,28 one considers the solution
y(j), j = 1, . . . , N, where N is the length of the trajectory y, and
constructs the two variables,
n
X
pα (n) = y(i) cos(iα),
i=1
n
X
qα (n) = y(i) sin(iα),
i=1
 N
where α ∈ (0, π) and n ∈ 1, 10 . In the numerical computations,
we have considered α = π/5. Regular geometric structures in the
FIG. 2. One-parameter bifurcation diagrams by varying r1 in (a), δ2 in (b), and (pα , qα ) plane correspond to regular dynamics and Brownian-like
τ in (c). Each bifurcation diagram is computed by considering the local maxima, motion to chaotic dynamics. From now on, we will drop the depen-
Tmax , of T of system (1). Here, we have used the parameter values in Table I,
keeping δ2 = 0.65 and τ = 23 in panel (a), r1 = 0.05 and τ = 23 in panel (b),
dence of p and q on α and will refer to the (p, q) plane for simplicity.
and r1 = 0.05 and δ2 = 0.65 in panel (c). By applying the 0–1 test in our case, we consider as y the solu-
tion T to system (1) that represents the evolution of the number of

Chaos 30, 123118 (2020); doi: 10.1063/5.0025510 30, 123118-4

Published under license by AIP Publishing.
 Chaos ARTICLE scitation.org/journal/cha

FIG. 5. Results from the application of the 0–1 test for the dynamics of the
time-delay TI system (1): Panels (a) and (b) show Kα as a function of δ2 for τ = 23
and Kα as a function of τ for δ2 = 0.65, respectively. Panel (c) shows the param-
eter space of Kα in the plane (τ , δ2 ). The color bar shows the values of Kα , where
Kα ≈ 0 corresponds to regular behavior (blue) and Kα ≈ 1 to chaotic behavior
N
(red). Here, Kα was computed for n ≤ 10 , where N = 1000 is the length of the
solution T of system (1).

To quantify the dynamics of the time-delay TI system (1), we

compute Kα in Fig. 5 as a function of δ2 in [0.2, 0.7] for fixed τ = 23
[panel (a)], as a function of τ in (0, 25] for fixed δ2 = 0.65 [panel
(b)], and in the parameter space (τ , δ2 ) in (0, 25] × [0.2, 0.7] in panel
(c). Panel (a) shows that depending on δ2 , Kα fluctuates around 0.5
FIG. 4. Characteristic (p, q) plots from the application of the 0–1 test for chaoticity
or increases to 1, indicating regular or chaotic dynamics, respec-
for the dynamics of system (1): Panels (a) and (c) show (p, q) plots for δ2 = 0.23
and 0.65, respectively, for fixed time-delay τ = 6, what correspond to regular tively. A similar situation appears for Kα in Fig. 5(b) where again,
dynamics. Panels (b) and (d) show (p, q) plots for τ = 8 and 23, respectively, for depending on τ , the system either lies in a regime of regular behavior
fixed δ2 = 0.65, what correspond to chaotic behavior. Note that a regular geomet- depicted by Kα < 0.5 or Kα around 1 after about τ ≈ 10, signaling
ric structure corresponds to regular behavior [panels (a) and (c)] and Brownian-like chaotic dynamics.
motion to chaotic dynamics [panels (b) and (d)]. Finally, Fig. 5(c) shows how Kα varies in the parameter space
(τ , δ2 ), similarly to our study in Fig. 3(c) for the two-parameter
bifurcation diagram. Comparing the results in Fig. 5(c) with those of
tumor cells in time. In Fig. 4, we plot (p, q) for δ2 = 0.23 [panel (a)] the two-parameter bifurcation diagram in Fig. 3(c), one can appre-
and δ2 = 0.65 [panel (c)], both for fixed τ = 6, and (p, q) for τ = 8 ciate high resemblance and similar conclusions with respect to the
[panel (b)] and τ = 23 [panel (d)], both for fixed δ2 = 0.65. Regu- dynamical properties of system (1). Concluding, these results show
lar dynamics are depicted by regular geometric structures, such as that the 0–1 test can successfully discriminate between regular and
those seen in panels (a) and (c). On the other hand, chaotic dynam- chaotic dynamics in the time-delay TI system (1).
ics is depicted by Brownian-like, irregular motion as can be seen in
panels (b) and (d). These results show the potential of the 0–1 test
to discriminate effectively between regular and chaotic dynamics in IV. DYNAMICAL COMPLEXITY ANALYSIS
the case of the time-delay TI system (1). Here, we investigate the dynamical complexity in system (1) by
Subsequently, the diffusive behavior on the (p, q) plane can be means of weighted recurrence plots and weighted entropy. Before
measured by means of the square dispersion, doing so, however, we will start by discussing our results on the long-
N term behavior of representative trajectories of system (1) and their
1 X
2
connection to cancer dynamics.
Mα (n) = lim pα (i + n) − pα (i)
N→∞ N
i=1

2  A. Phase space analysis and connection to cancer
+ qα (i + n) − qα (i) ,
dynamics
where n is very small compared to the length of the trajectory N (i.e, We start by studying the asymptotic behavior of system
n  N). The limit is assured by computing Mα (n) only for n ≤ ncut , (1). In particular, panels (a), (b), and (c) in Fig. 6 show three-
where ncut  N. In particular, as the test for chaos depends on the dimensional projections of regular and chaotic trajectories. Panels
N
growth rate of Mα as a function of n, we consider ncut = 10 . (a) and (b) show, for δ2 = 0.65 and r1 = 0.05, period 1 and 2
Following Refs. 27 and 28, regular and chaotic dynamics can be trajectories for τ = 5 and τ = 11, respectively, which correspond to
quantified by the asymptotic growth rate of tumor-immunoediting that goes through an equilibrium process of
log (Mα (n)) dual host protection and tumor-promoting actions in the immune
Kα = lim , system.12,14 Panel (c) shows a chaotic structure for τ = 23. These
n→∞ log(n)
results provide evidence of a nonlinear relation among immune,
where n → ∞ implies that all possible time lags are considered in host, and tumor cells.18 In Fig. 6(c), it can be observed that host cells,
the computation of Mα . Consequently, Kα ≈ 0 denotes regular and H, and tumor cells, T, nearly attain their maximum cell numbers,
Kα ≈ 1 denotes chaotic behavior. i.e., their carrying capacity. This represents a long-term, chaotic, and

Chaos 30, 123118 (2020); doi: 10.1063/5.0025510 30, 123118-5

Published under license by AIP Publishing.
 Chaos ARTICLE scitation.org/journal/cha

(a) 1000 1
(b) 1000 1

800 0.8 800 0.8

600 0.6 600 0.6

Time

Time
400 0.4 400 0.4

200 0.2 200 0.2

0 0
0 200 400 600 800 1000 0 200 400 600 800 1000
Time Time

FIG. 7. Recurrence plots for the regular and chaotic dynamics in Subsection. III C.
Panels (a) and (b) show the recurrence plots for the regular dynamics for
τ = 6 and chaotic dynamics for τ = 23 with δ2 = 0.65. Note the difference in
the patterns in the two panels. The color bars show the values of the weights wij
of the weighted N × N matrix W, where N = 1000.

For periodic trajectories, the occurrence of identical states

results in wij = 1, and for chaotic trajectories, wij tends to zero due
FIG. 6. Three-dimensional projections of periodic and Poincaré map of a chaotic to exponential divergence. We note that definition (3) results in a
trajectory of the time-delay TI system (1): Panels (a)–(c) show three-dimensional symmetric matrix W. Thus, different dynamical behaviors can be
projections of trajectories of periods 1 and 2 and a chaotic trajectory for τ = 5, encoded in the weighted matrix W with entries in [0, 1], and dynam-
11, and 23, respectively, where δ2 = 0.65 and r1 = 0.05. Panel (d) shows the
Poincaré map of the chaotic trajectory in panel (c). Note that all other parameters
ical complexity can be observed based on the entries of W, as we
used are those in Table I. discuss below.
In the following, we consider T-solution component of sys-
tem (1) as the trajectory y and compute W using Eq. (3). We then
thus unpredictable growth in the number of tumor cells, pinpointing use the function imagesc in Matlab to scale the data in W as a col-
to cancer dynamics.20,39 orful image. Each entry of W specifies the color for one pixel of
We complement these results by studying the Poincaré map of the image. The resulting image is an N × N grid of pixels, where
the local maxima, Tmax , of the chaotic solution T in panel (c) that N = 1000. The results of the weighted-recurrence computations are
represents the evolution of the number of tumor cells in time. In presented in Fig. 7, where panel (a) shows a periodically repeated
particular, the Poincaré map of the chaotic structure in panel (c) for pattern for δ2 = 0.23 and τ = 6, for the regular dynamics discussed
τ = 23, δ2 = 0.65, and r1 = 0.05 shows scattered points in panels in Subsection. III C. This pattern is due to the periodic trajectory
(d), providing further evidence of its chaotic nature. These results visiting the same region in the state-space over and over in time,
are in line with those from the 0–1 test in Subsection. III C. what is a recurrence. On the other hand, panel (b) shows different
morphological patterns for δ2 = 0.23 and τ = 23, which amount to
chaotic dynamics and indicate a disordered structure (see Subsec-
B. Weighted recurrence entropy analysis tion. III C). It is worth noting that the results from the recurrence
Finally, we turn to the study of dynamical complexity in sys- plots are in agreement with the regular and chaotic nature of the tra-
tem (1) by employing the computation of weighted entropy based jectories depicted by the 0–1 test in Subsection. III C and Poincaré
on weighted recurrence plots.33,37 maps in Subsection. IV A.
In particular, in an m-dimensional space, for a given trajectory Next, we proceed to the study of the dynamical complexity in
y ∈ Rm , the weighted recurrence is defined as the matrix W, where system (1) and introduce the weighted entropy, Ewr , based on the
its entries, wij , are defined by weighted recurrence matrix W in Eq. (3). The weighted entropy Ewr 33
is defined by
wij = e−kyi −yj k , i, j = 1, . . . , N, (3) X
Ewr = − p(ek ) log p(ek ), (4)
where k · k denotes the Euclidean distance and N denotes the length ek ∈E
of the trajectory y. Here, wij is the inverse exponential of the
Euclidean distance between pairs of points yi and yj at times i and where p(ek ) is the probability of ek to occur in E, where
j, respectively. ( N
)
The weights wij encode the information whether points i and j 1 X
E = ek : e k = wki ; k = 1, . . . , N .
are close in the m-dimensional space at times i and j by measuring N i=1
the distances kyi − yj k, computed as the inverse of the exponen-
tial function, normalized in [0, 1]. In this context, 0 corresponds For the computation of Ewr , the probabilities in Eq. (4) were
to divergent and 1 to identical states. Consequently, the weights wij computed by constructing the histogram of the values of the set
describe the disorder in the system. E using the bin-counting method and 20 bins. The histogram was

Chaos 30, 123118 (2020); doi: 10.1063/5.0025510 30, 123118-6

Published under license by AIP Publishing.
 Chaos ARTICLE scitation.org/journal/cha

(a) 3 (b) 3 (c) 0.7 2.8

performing 0–1 tests for trajectories in the characteristic of regu-
2.8 2.6

2.6
2.8
0.6
2.4
lar and chaotic regimes. In this context, regular dynamics indicates
0.5 2.2
an equilibrium state of cancer immunoediting on the simultane-
Ewr

Ewr
2.4

2
2.6
2
2.2 0.4
1.8
ous stage in host-protective and cancer-promoting actions of the
2.4
2 0.3 1.6 immune system. Further clinical investigation is thus necessary.12
1.8 1.4
0.2 0.3 0.4 0.5 0.6 0.7
2.2
0 5 10 15 20 25
0.2
0 5 10 15 20 25
Interestingly, chaotic dynamics is also found here and corresponds
2
to long-term phenomena of cancer relapse.14,39
From the viewpoint of dynamical systems, regular patterns cor-
FIG. 8. Recurrence entropy plots: Panel (a) shows Ewr as a function of δ2 in
[0.2, 0.7] for fixed τ = 23 and panel (b), Ewr as a function of τ in (0, 25] for fixed respond to self-organization. This means that cellular interactions
δ2 = 0.65. Panel (c) shows the plot of Ewr for different combinations of τ and δ2 . are more robust. On the other hand, complex behavior is asso-
Note that we have used 20 bins in the computations of the probabilities in Ewr . ciated with weak interactions generating irregular patterns.40 The
entropy computed from recurrence plots provides further evidence
of the structural disorder in the delayed system studied here. Our
further normalized so that the probabilities sum up to 1. The depen- work provides evidence that high entropy production, driven by
dence of Ewr on δ2 and τ is shown in panels (a) and (b) in Fig. 8. chaoticity, corresponds to higher complexity, which leads to abrupt,
Panel (a) shows the dependence of Ewr on δ2 in [0.2, 0.7] for fixed rapid proliferation of cancer cells. This is the result of the loss of
τ = 23, and panel (b) shows the dependence on τ in (0, 25] for autonomic and synchronized balance in the tumorigenic (low ener-
fixed δ2 = 0.65. A similar trend in fluctuations can be seen in pan- getic) and tumoricidal (high energetic) properties of immunity in
els (a) and (c) in Fig. 5, showing that chaotic behavior corresponds defending cancerous cells.41,42
to higher entropy, depicted by Ewr . We also plot the weighted recur- Here, we have shown that our model exhibits regular and
rence entropy Ewr in Fig. 8(c) to study its dependence on τ and δ2 on chaotic dynamics, depending on the choice of parameters and time-
(0, 25] × [0.2, 0.7]. Comparing Figs. 5(c) and 8(c), we observe sim- delay. Chaos is attributed to sensitive dependence of the trajectories
ilar pattern variations between Kα and Ewr . These results provide of the system on initial conditions. This sensitivity gives rise to local
additional evidence of the dynamical complexity in the system for instability. Chaotic behavior gives rise to increased entropy produc-
different parameter regimes and settings. tion (as opposed to entropy production in regular dynamics) and to
Our study has shown that the chaotic behavior in system (1) complexity in the system. Thus, higher entropy production is driven
is the result of disorder and chaotic oscillations, which are asso- by chaotic dynamics, which produces complexity in the system. In
ciated with the dynamics of cancerous cells.24 In contrast, regular the case of chaotic trajectory in our model, the ability to predict
dynamics, as those discussed for system (1), are indicative of dual the proliferation rate of tumor cell is lost, and thus, the behavior of
host protection and tumor-promoting action that lead to cancer the trajectory becomes unpredictable. However, the trajectories do
immunoediting.5 Disorderly rapid oscillations of cancerous cells are not escape to infinity due to bounded solutions for all finite times.
related to increased rates of entropy production and complexity Chaoticity, thus, means that the evolution of the number of tumor
in the tumor-immune interaction.36 This is corroborated by our cells, T, becomes unpredictable in time. Typically, the innate and
finding that the rate of entropy production of the cancerous cells adaptive immune cells suppress the progression of tumor prolifer-
is higher than that of the immune and healthy cells. Thus, the ation by either eliminating tumor cells or by attempting to regulate
increase in entropy production, driven by chaoticity, is related to their outgrowth (tumor-immune interplay). When the competition
the increased level of proliferation of cancerous cells. Consequently, between tumor cells and their micro-environment (attributed to
higher entropies reflect rapid and highly uncontrolled proliferation host cells, effector cells, and cytokines) starts to lose its efficiency,
of cancer cells.31,40 the proliferation of tumor cells is no longer strongly maintained,
and thus, fast tumor growth is possible. Thus, in the framework
of tumor-immune dynamics, sensitivity to initial conditions gives
V. DISCUSSION rise to local instability, which leads to unpredictability, increased
Cellular interactions in the cancer-immune system give rise to production of entropy, complexity, and invokes the uncontrollable,
different cancer-growth patterns and contribute to highly compli- unpredictable, proliferation of cancer-specific patterns24–26 and also
cated and not well-understood processes. Here, we sought to explore to long-term cancer relapse.39
the dynamics and complexity of a time-delay TI system. To this end, On the other hand, regular dynamics indicates an equilib-
we have performed stability and bifurcation analyses for different rium state of cancer immunoediting on the simultaneous process
parameter settings. The bifurcation diagrams show a stable, peri- in host-protective and cancer-promoting actions.12 This might lead
odic, and chaotic behavior. Furthermore, two-dimensional bifurca- to tumor dormancy, i.e., to a non-cancerous stage. Our study has
tion diagrams in the form of parametric plots revealed regions of shown that the chaotic behavior in our system is the result of dis-
periodic and chaotic behavior that were further explored. order and chaotic oscillations that are associated with the dynamics
The appearance of chaotic behavior pinpoints the exis- of cancerous cells. In contrast, regular dynamics are indicative of
tence of metastatic cancer in the viewpoint of oncology.39 Two- dual host protection and tumor-promoting action that lead to can-
dimensional bifurcation diagrams have identified regions of cer immunoediting.5 Disorderly rapid oscillations of cancerous cells
low- and high-order periodic behavior leading to period-doubling are related to increased rates of entropy production and complex-
cascades to chaos. In our work, the regular and chaotic behavior ity in the tumor-immune interaction.36 This is corroborated by our
revealed in the bifurcation diagrams was further corroborated by finding that the rate of entropy production of the cancerous cells

Chaos 30, 123118 (2020); doi: 10.1063/5.0025510 30, 123118-7

Published under license by AIP Publishing.
 Chaos ARTICLE scitation.org/journal/cha

12
is higher than that of the immune and healthy cells. Thus, the C. M. Koebel, W. Vermi, J. B. Swann, N. Zerafa, S. J. Rodig, L. J. Old, M. J.
increase in entropy production, driven by chaoticity, is related to Smyth, and R. D. Schreiber, Nature 450, 903–907 (2007).
13
the increased level of proliferation of cancerous cells. Consequently, M. D. Vesely, M. H. Kershaw, R. D. Schreiber, and M. J. Smyth, Annu. Rev.
Immunol. 29(1), 235–271 (2011).
higher entropies reflect rapid and highly uncontrolled proliferation 14
M. Moghtadaei, M. R. H. Golpayegani, and R. Malekzadeh, J. Theor. Biol. 334,
of cancer cells.31,40 130–140 (2013).
Finally, our study has revealed that the complexity of the 15
H. M. Byrne, Math. Biosci. 144(2), 83–117 (1997).
16
tumor-immune interaction builds up with the increase in entropy P. Das, P. Das, and S. Das, Appl. Math. Comput. 361, 536–551 (2019).
17
production of cancerous cells attributed to unpredictable growth in P. Das, P. Das, and S. Das, Math. Model Nat. Phenom. 15, 45 (2020).
18
the number of tumor cells, that is to chaotic dynamics. The pro- S. Abernethy and R. J. Gooding, Physica A 507(C), 268–277 (2018).
19
M. Itik and S. Banks, Int. J. Bifurcat. Chaos 20, 71–79 (2010).
duction of entropy, chaotic behavior, and complex dynamics are 20
C. Letellier, F. Denis, and L. A. Aguirre, J. Theor. Biol. 322, 7–16 (2013).
concepts that can guide computational biologists elucidate associ- 21
P. Das, S. Mukherjee, and P. Das, Chaos Soliton Fractals 128, 297–305 (2019).
ations between oscillation patterns and variations in parameters to 22
D. Ghosh, S. Khajanchi, S. Mangiarotti, F. Denis, S. K. Dana, and C. Letellier,
develop cancer treatments and novel drug therapies. Thus, it would BioSystems 158, 17–30 (2017).
23
be useful in a future publication to compare model solutions studied S. Khajanchi, M. Parc, and D. Ghosh, Chaos 28(10), 103101 (2018).
24
here with oncology and immunological data to see how closely they J. Post, R. J. Sklarew, and J. Hoffma, Cancer 39(4), 1500–1507 (1977).
25
P. Das, S. Mukherjee, P. Das, and S. Banerjee, Nonlinear Dyn. 101(1), 675–685
match. (2020).
26
P. Das, P. Das, and S. Mukherjee, Physica A 541, 123603 (2019).
27
DATA AVAILABILITY G. A. Gottwald and I. Melbourne, Proc. R. Soc. Lond. A 460, 603–611 (2004).
28
G. A. Gottwald and I. Melbourne, Phys. Rev. E 77, 028201 (2008).
The data that support the findings of this study are available 29
W. Just and H. Kantz, J. Phys. A Math. Theor. 33(1), 163–170 (1999).
within the article. 30
D. Daems and G. Nicolis, Phys. Rev. E 59, 4000–4006 (1999).
31
C. K. Volos, S. Jafari, J. Kengne, J. M. Munoz-Pacheco, and K. Rajagopal,
Entropy 21(4), 370 (2019).
REFERENCES 32
S. Guiasu, Rep. Math. Phys. 2(3), 165–179 (1971).
1 33
D. Helbing, D. Brockmann, T. Chadefaux, K. Donnay, U. Blanke, O. Woolley- B. Yan, S. Mukherjee, and S. He, Eur. Phys. J. Spec. Top. 228(12), 2769–2777
Meza, M. Moussaid, A. Johansson, J. Krause, S. Schutte, and M. Perc, J. Stat. Phys. (2019).
34
158(3), 735–781 (2015). N. Marwan, M. C. Romano, M. Thiel, and J. Kurths, Phys. Rep. 438(5), 237–329
2 (2007).
T. Bountis, J. Johnson, A. Provata, and G. Tsironis, Eur. Phys. J. Spec. Top.
35
225(6), 883–890 (2016). S. Banerjee, S. K. Palit, S. Mukherjee, M. R. K. Ariffin, and L. Rondoni, Chaos
3 26(3), 033105 (2016).
R. P. Araujo and D. L. S. McElwain, Bull. Math. Biol. 66(5), 1039 (2004).
4 36
R. A. Gatenby and P. K. Maini, Nature 421(6921), 321 (2003). A. Kleidon, Y. Kleidon, and P. M. Cox, Philos. Trans. R. Soc. B 365(1545),
5 1297–1302 (2010).
G. P. Dunn, A. T. Bruce, H. Ikeda, L. J. Old, and R. D. Schreiber, Nat. Immunol.
37
3(11), 991–998 (2002). D. Eroglu, T. K. DM. Peron, N. Marwan, F. A. Rodrigues, L. da F. Costa,
6 M. Sebek, I. Z. Kiss, and J. Kurths, Phys. Rev. E 90, 042919 (2014).
R. Weinberg, Biology of Cancer (Garland Science, London, 2013).
7 38
S. Rakoff-Nahoum, Yale J. Biol. Med. 79, 123–130 (2006). P. Das, S. Das, R. K. Upadhyay, and P. Das, Chaos Solitons Fractals 136, 109806
8 (2020).
K. P. Wilkie and P. Hahnfeldt, Bull. Math. Biol. 79, 1426–1448 (2017).
9 39
C. Letellier, S. K. Sasmal, C. Draghi, F. Denis, and D. Ghosh, Chaos Solitons A. Cucuianu, Nat. Med. 4(12), 1342 (1998).
40
Fractals 99, 297–311 (2017). D. S. Coffey, Nat. Med. 4(8), 882–885 (1998).
10 41
S. Khajanchi and D. Ghosh, Appl. Math. Comput. 271, 375–388 (2015). W. Ritchie, S. Granjeaud, D. Puthier, and D. Gautheret, PLoS Comput. Biol.
11 4(3), e1000011 (2008).
K. P. Wilkie, A Review of Mathematical Models of Cancer-Immune Interactions
42
in the Context of Tumor Dormancy (Springer, New York, NY, 2013). M. Khatami, Clin. Transsl. Med. 7(1), 20 (2018).

Chaos 30, 123118 (2020); doi: 10.1063/5.0025510 30, 123118-8

Published under license by AIP Publishing.