[1]Enrique C. Gabrick These authors contributed equally to this work.

These authors contributed equally to this work. These authors contributed
equally to this work. These authors contributed equally to this work. These
authors contributed equally to this work. These authors contributed equally to
this work. These authors contributed equally to this work. These authors
contributed equally to this work. These authors contributed equally to this
work. These authors contributed equally to this work. These authors
contributed equally to this work.

1]Graduate Program in Science, State University of Ponta Grossa, 84030-900,
Ponta Grossa, PR, Brazil 2]Department of Mathematics and Statistics, State
University of Ponta Grossa, 84030-900, Ponta Grossa, PR, Brazil 3]Department
of Physics, State University of Ponta Grossa, 84030-900, Ponta Grossa, PR,
Brazil 4]Department of Physiology and Pharmacology, State University of New
York Downstate Health Sciences University, 11203, Brooklyn, NY, USA
5]University Center UNIFATEB, 84266-010, Telêmaco Borba, PR, Brazil
6]Institute of Physics, University of São Paulo, 05508-090, São Paulo, SP,
Brazil 7]Department of Physics, Federal University of Paraná, 82590-300,
Curitiba, PR, Brazil

# Fractional dynamics and recurrence analysis in cancer model

<EMAIL_ADDRESS>Matheus R. Sales  Elaheh Sayari  José Trobia  Ervin K.
Lenzi  Fernando da S. Borges  José D. Szezech Jr  Kelly C. Iarosz  Ricardo L.
Viana  Iberê L. Caldas  Antonio M. Batista  [  [  [  [  [  [  [

###### Abstract

In this work, we analyze the effects of fractional derivatives in the chaotic
dynamics of a cancer model. We begin by studying the dynamics of a standard
model, i.e., with integer derivatives. We study the dynamical behavior by
means of the bifurcation diagram, Lyapunov exponents, and recurrence
quantification analysis (RQA), such as the recurrence rate (RR), the
determinism (DET), and the recurrence time entropy (RTE). We find a high
correlation coefficient between the Lyapunov exponents and RTE. Our
simulations suggest that the tumor growth parameter ($\rho_{1}$) is associated
with a chaotic regime. Our results suggest a high correlation between the
largest Lyapunov exponents and RTE. After understanding the dynamics of the
model in the standard formulation, we extend our results by considering
fractional operators. We fix the parameters in the chaotic regime and
investigate the effects of the fractional order. We demonstrate how fractional
dynamics can be properly characterized using RQA measures, which offer the
advantage of not requiring knowledge of the fractional Jacobian matrix. We
find that the chaotic motion is suppressed as $\alpha$ decreases, and the
system becomes periodic for $\alpha\lessapprox 0.9966$. We observe limit
cycles for $\alpha\in(0.9966,0.899)$ and fixed points for $\alpha<0.899$. The
fixed point is determined analytically for the considered parameters. Finally,
we discover that these dynamics are separated by an exponential relationship
between $\alpha$ and $\rho_{1}$. Also, the transition depends on a supper
transient which obeys the same relationship.

###### keywords:

Cancer model, Fractional Calculus, Recurrence analysis

##  1 Introduction

Cancer is a set of diseases that arises from the abnormal growth and
uncontrolled division of body cells. It can spread to the body’s cells,
causing many deaths [1]. Each year the American Cancer Society estimates the
number of new cancer cases. For 2020 there was an estimated number of 19.3
million new cancer cases and almost 10.0 million deaths from cancer [2]. Only
in the United States of America were projected 1,918,030 new cancer cases and
609,360 cancer deaths for 2022 [3]. In this way, cancer is a crucial public
health problem worldwide [4] that requires many efforts to understand the
mechanism behind the illness and improve the treatment methods [5, 6, 7].
There are several methods to study the dynamics of cancer cell proliferation
and one of the most successful method is through mathematical models.

Mathematical model is a powerful tool for understanding cancer dynamics [8, 9]
and simulate treatment measures, such as effects of drug resistance [10, 11],
immunotherapy [12], chemotherapy [13, 14], radiotherapy [15], biochemotherapy
[16] and predict fluctuations associated with the growth rate of cancer cells
[17]. A fundamental question that arises is how to understand the dynamics
among healthy cells, the immune system and the tumor cells. In light of this
question, several models have been proposed and studied [18, 19, 20, 21, 22,
23, 24].

Dehingia et al. [25] studied the effect of time delay in tumor-immune
interaction and stimulation process. They obtained conditions for the
existence of equilibrium points and Hopf bifurcation. Also, they derived
conditions for periodic solutions. Díaz-Marín et al. [26] proposed a model to
describe cell population dynamics in a tumor with periodic radiation as
treatment. As global attractors, they found almost periodic solutions or the
vanishing equilibrium. López et al. [27] formulated a model for tumor growth
in the presence of cytotoxic chemotherapeutic agents. Their model allows an
investigation of the Norton-Simon hypothesis in the context of dose-dense
chemotherapy. In another work, López et al. [28] validated a model by means of
experimental results. In this model, they considered tumor growth, including
tumor-healthy cell interactions, immune response, and chemotherapy. Most of
these models deal with three or more dimensional systems, where chaotic
behavior is a possibility [29].

A three-dimensional cancer model with chaotic dynamics was proposed by Itik
and Banks [30]. In their model, they considered the interactions of tumor
cells with healthy host cells and immune systems. Through the calculation of
Lyapunov exponents, they showed the existence of chaotic attractors for some
parameters. Letellier et al. [31], made a topological and observability
analysis of this model. They also investigated the equilibrium points and the
bifurcation diagrams. Khajanchi et al. [32] studied a similar model with time
delay. They analyzed the existence and stability of biologically feasible
points and the emergence of Hopf bifurcations. Chaos in the 3-cell cancer
model was also studied by Abernethy and Gooding [33], Gallas et al. [34],
Khajanchi [35] and others [36, 37, 38, 39].

Although the literature about the model proposed by Itik and Banks [30] is
extensive, the works are restricted to integer-order differential equations.
However, this formulation does not incorporate non-local effects. To do that,
it is necessary to consider fractional differential equations [40]. In
general, non-local operators are more efficient in describing some situations
since they capture non-localities and have memory effects [41]. Fractional
operators have been used to model many real problems, such as photo acoustic
[42, 43], viscoelastic properties [44], Quantum Mechanics [45, 46, 47],
Epidemiology [48, 49], Ecology [50, 51], Duffing oscillator [52] and many
others [53, 54, 55].

The dynamical behavior of the cancer model proposed by Itik and Banks [30] was
studied in light of fractional operators in Ref. [56]. It was considered three
different fractional derivatives formulations: the power-law [57] (singular
kernel), the exponential [58], and the Mittag-Leffler [59]. The results show
that the dynamics is changed as a consequence of the extension to fractional
differential operators. Ghanbari [60] also explored the same model and
investigated the influence of power-law, exponential decay-law, and Mittag-
Leffler in the fractal-fractional approach. These researchers obtained
conditions for the existence and uniqueness of the solutions and expanded a
numerical method to study the dynamical behavior. The results showed that the
dynamics change from chaotic to a limit cycle depending on the factional
order. Naik et al. [61], considering an extension given by the Caputo
derivative, analyzed the stability of the model. In a numerical scheme, they
reported that the system goes to a limit cycle in a chaotic regime for the
standard model. Xuan et al. [62], using the Caputo fractal-fractional
derivative in the cancer model, demonstrated that hidden attractors emerge
under fractional operators. To study the system’s complexity, they studied
bifurcation diagrams and stability. Their findings also showed that the system
converges to a limit cycle attractor when the fractional order is decreased.
These works extend the model of Itik and Banks [30]. Extensions of other
cancer models can be found in Refs. [63, 64, 65, 66, 67, 68, 69, 70, 71, 72].

In this work, we analyze the influence of fractional operators in a cancer
model [30]. As a definition of the fractional operator, we follow the Caputo
scheme [41], with an order equal to $\alpha$. We consider the parameters in
which the trajectories are chaotic [31]. As a novelty, we characterize the
model using recurrence quantification, such as recurrence rate (RR),
recurrence time entropy (RTE), and determinism (DET). Furthermore, we localize
the $\alpha$ value in which the system transits from chaos to limit cycle
utilizing RTE. Also, considering the mean square deviation of the time series,
we construct the plane parameters composed of $\rho_{1}\times\alpha$, where
$\rho_{1}$ is the tumor growth. Our results suggest an exponential relation
between $\alpha$ and tumor growth rate $\rho_{1}$ in which the dynamics
transit from the limit cycle to a fixed point. This transition depends on a
supper transient time that follows the same relationship between $\rho_{1}$
and $\alpha$.

Our work is organized as follows. In Section 2, we revisit the integer-order
model analyzing the fixed points and their stabilities. Section 3 discusses
the fractional influences in the cancer model. Our conclusions are drawn in
Section 4.

##  2 Standard model

We consider a cancer model that describes the interactions among host cells
($H$), effector immune cells ($E$) and tumor cells ($T$), governed by the
following equations:

$\displaystyle\frac{dH}{dt}$ $\displaystyle=$
$\displaystyle\rho_{1}H\left(1-\frac{H}{\kappa_{1}}\right)-\alpha_{13}TH,$ (1)
$\displaystyle\frac{dE}{dt}$ $\displaystyle=$
$\displaystyle\rho_{2}\frac{TE}{T+\kappa_{2}}-\alpha_{23}TE-\delta_{2}E,$ (2)
$\displaystyle\frac{dT}{dt}$ $\displaystyle=$
$\displaystyle\rho_{3}T\left(1-\frac{T}{\kappa_{3}}\right)-\alpha_{31}TH-\alpha_{32}TE.$
(3)

Equation 1 describes the growth of the host cells by a logistic function with
a rate equal to $\rho_{1}$ and biotic capacity equal to $\kappa_{1}$. The host
cells are killed by tumor cells. This term is represented by the interaction
$-\alpha_{13}TH$, where $\alpha_{13}$ is the host cell killing rate by tumor
cells. Equation 2 gives the rate at which immune cells appear. These cells are
the tumor antigens whose growth depends on the quantity of $T$, the growth
rate of effector immune cells ($\rho_{2}$), and a positive constant
$\kappa_{2}$. The tumor cell not only contributes to the growth of immune
cells but also to its death. In the term $-\alpha_{23}TE$, $\alpha_{23}$ is
the rate at which tumor cells inhibit the immune cells. In addition, the $E$
cells die according to $\delta_{2}$. The last equation, Eq. 3, gives the
dynamics of the tumor cells. Similarly, with the $H$ cells, the tumor growth
is given by a logistic function with a growth rate equal to $\rho_{3}$ and
biotic capacity equal to $\kappa_{3}$. The $T$ cells are killed by $H$ cells
at a rate equal to $\alpha_{31}$ and by the $E$ cells with a rate equal to
$\alpha_{32}$.

Considering the transformation $(H,E,T)\rightarrow(x,y,z)$ [30], the
normalized equations are

$\displaystyle\frac{dx}{dt}$ $\displaystyle=$
$\displaystyle\rho_{1}x(1-x)-\alpha_{13}xz,$ (4) $\displaystyle\frac{dy}{dt}$
$\displaystyle=$
$\displaystyle\frac{\rho_{2}yz}{1+z}-\alpha_{23}yz-\delta_{2}y,$ (5)
$\displaystyle\frac{dz}{dt}$ $\displaystyle=$ $\displaystyle
z(1-z)-xz-\alpha_{32}yz.$ (6)

A schematic representation of the model is displayed in Fig. 1. The red,
green, and blue circles represent the $x$, $y$, and $z$ cells, respectively.
The arrows indicate the interactions. If the interaction is constructive
(growth cells), we denote with a signal $+$. If the interaction is destructive
(kill cells), we denote it by a signal $-$. The $x$ cells growth with a rate
$\rho_{1}$, which is indicated by $+\rho_{1}$ in Fig. 1. The $x$ cells are
killed by the tumor cells at a rate equal to $\alpha_{13}$. The tumor cells
grow according to a rate $\rho_{3}$ and are killed by $x$ cells at a rate
$\alpha_{31}$ and by $y$ at a rate equal to $\alpha_{32}$. The $y$ cells grow
when interacting with tumor cells by a rate $\rho_{2}$ and are killed by $z$
cells according to a rate equal to $\alpha_{23}$. Also, the $y$ cells have a
natural death $\delta_{2}$. The immune cells do not interact with the host
cells. In this way, the tumor cells are the generalist competitors and the
other two are specialist competitors [31]. In this type of interaction, there
is no particular prey.

Figure 1: Schematic representation of the cancer model. $x$, $y$ and $z$
represent the host, immune, and tumor cells. $\rho_{1}$ is the growth rate of
$x$, $\alpha_{13}$ is the death rate of $x$ due to $z$, $\alpha_{31}$ is the
death rate of $z$ due to $x$, $\rho_{3}$ is the growth rate of $z$,
$\alpha_{32}$ is the death rate of $z$ due to $y$, $\rho_{2}$ is the growth
rate of $y$ due to $z$, $\alpha_{32}$ is the death rate of $y$ due to $z$ and
$\delta_{2}$ is the natural death of $y$.

The fixed point solutions of this model are found by solving:

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\rho_{1}x(1-x)-\alpha_{13}xz,$ (7) $\displaystyle 0$
$\displaystyle=$
$\displaystyle\frac{\rho_{2}yz}{1+z}-\alpha_{23}yz-\delta_{2}y,$ (8)
$\displaystyle 0$ $\displaystyle=$ $\displaystyle z(1-z)-xz-\alpha_{32}yz,$
(9)

admitting seven solutions, that are $F_{i}\equiv(x^{i},y^{i},z^{i})$, where
$i=1,...,7$. Solving Eqs. 7, 8 and 9, we obtain

$\displaystyle F_{1}$ $\displaystyle=$ $\displaystyle(0,0,0),\ F_{2}=(0,0,1),\
F_{3}=(1,0,0)$ (10) $\displaystyle F_{4}$ $\displaystyle=$
$\displaystyle\Bigg{(}0,\frac{\sqrt{\left(\alpha_{23}+\delta_{2}-\rho_{2}\right){}^{2}-4\alpha_{23}\delta_{2}}+3\alpha_{23}+\delta_{2}-\rho_{2}}{2\alpha_{23}\alpha_{32}},$
(11) $\displaystyle-$
$\displaystyle\frac{\sqrt{\left(\alpha_{23}+\delta_{2}-\rho_{2}\right){}^{2}-4\alpha_{23}\delta_{2}}+\alpha_{23}+\delta_{2}-\rho_{2}}{2\alpha_{23}}\Bigg{)}$
$\displaystyle F_{5}$ $\displaystyle=$
$\displaystyle\Bigg{(}\frac{\alpha_{13}\left(\sqrt{-2\rho_{2}\left(\alpha_{23}+\delta_{2}\right)+\left(\alpha_{23}-\delta_{2}\right){}^{2}+\rho_{2}^{2}}+\alpha_{23}+\delta_{2}-\rho_{2}\right)}{2\alpha_{23}\rho_{1}}+1,$
(12) $\displaystyle-$
$\displaystyle\frac{\left(\alpha_{13}-\rho_{1}\right)\left(\sqrt{\left(\alpha_{23}+\delta_{2}-\rho_{2}\right){}^{2}-4\alpha_{23}\delta_{2}}+\alpha_{23}+\delta_{2}-\rho_{2}\right)}{2\alpha_{23}\alpha_{32}\rho_{1}},$
$\displaystyle-$
$\displaystyle\frac{\sqrt{\left(\alpha_{23}+\delta_{2}-\rho_{2}\right){}^{2}-4\alpha_{23}\delta_{2}}+\alpha_{23}+\delta_{2}-\rho_{2}}{2\alpha_{23}}\Bigg{)},$
$\displaystyle F_{6}$ $\displaystyle=$
$\displaystyle\Bigg{(}0,-\frac{\sqrt{\left(\alpha_{23}+\delta_{2}-\rho_{2}\right){}^{2}-4\alpha_{23}\delta_{2}}-3\alpha_{23}-\delta_{2}+\rho_{2}}{2\alpha_{23}\alpha_{32}},$
(13)
$\displaystyle\frac{\sqrt{\left(\alpha_{23}+\delta_{2}-\rho_{2}\right){}^{2}-4\alpha_{23}\delta_{2}}-\alpha_{23}-\delta_{2}+\rho_{2}}{2\alpha_{23}}\Bigg{)},$
$\displaystyle F_{7}$ $\displaystyle=$
$\displaystyle\Bigg{(}\frac{\alpha_{13}\left(-\sqrt{-2\rho_{2}\left(\alpha_{23}+\delta_{2}\right)+\left(\alpha_{23}-\delta_{2}\right){}^{2}+\rho_{2}^{2}}+\alpha_{23}+\delta_{2}-\rho_{2}\right)}{2\alpha_{23}\rho_{1}}+1,$
(14)
$\displaystyle\frac{\left(\alpha_{13}-\rho_{1}\right)\left(\sqrt{\left(\alpha_{23}+\delta_{2}-\rho_{2}\right){}^{2}-4\alpha_{23}\delta_{2}}-\alpha_{23}-\delta_{2}+\rho_{2}\right)}{2\alpha_{23}\alpha_{32}\rho_{1}},$
$\displaystyle\frac{\sqrt{\left(\alpha_{23}+\delta_{2}-\rho_{2}\right){}^{2}-4\alpha_{23}\delta_{2}}-\alpha_{23}-\delta_{2}+\rho_{2}}{2\alpha_{23}}\Bigg{)}.$

To simplify the analysis of the fixed points, the parameters are equal to
$\alpha_{13}=1.5$, $\rho_{2}=4.5$, $\alpha_{23}=0.2$, $\delta_{2}=0.5$ and
$\alpha_{32}=2.5$. $\rho_{1}$ is maintained as our control parameter in the
entire work. For these values, the fixed points are $F_{1}=(0,0,0)$,
$F_{2}=(0,0,1)$, $F_{3}=(1,0,0)$, $F_{4}=(0,0.3469,0.1325)$,
$F_{5}=\frac{1}{\rho_{1}}(\rho_{1}-0.1987,0.0530(1.5-\rho_{1}),0.1325\rho_{1})$,
$F_{6}=(0,-7.1470,18.8675)$ and
$F_{7}=\frac{1}{\rho_{1}}(\rho_{1}-28.3012,7.5470(1.5-\rho_{1}),18.8675\rho_{1})$.
$F_{1}$ is a trivial solution and exhibits a population without cells. The
eigenvalues associated with the Jacobian calculated in $F_{1}$ are
$(1,-0.5,\rho_{1})$. The signs depend on $\rho_{1}$, which is always positive.
Then, we find two unstable and one stable point, which characterize a saddle.
$F_{2}$ is the point that represents a situation where the cancer cells occupy
the whole host. The Jacobian eigenvalues in $F_{2}$ are
$(1.55,\rho_{1}-1.5,-1)$. Considering $\rho_{1}=0.6$, we obtain
$(1.55,-0.9,-1)$, which characterizes a saddle point. $F_{3}$ corresponds to
the cancer-free solution. The Jacobian eigenvalues calculated at this point
are $\left(0,-0.5,-\rho_{1}\right)$, as $\rho_{1}\geq 0$, we have two negative
eigenvalues and one equal zero. Thus, this point is a saddle. In $F_{4}$, all
coordinates are positives; all solutions are in the positive octant. The
Jacobian eigenvalues in this point are
$(\rho_{1}-0.198754,\,\,-0.0662515+0.6131239\,i,\,\,-0.0662515-0.6131239\,i)$.
The $F_{5}$ point has a biological significance for $\rho_{1}\in[0.1987,1.5]$.
Selecting $\rho_{1}=0.6$, we obtain $(0.6688,0.0795,0.1325)$ which is an
unstable point. The $F_{6}$ point has a negative term associated with the
effector immune cells. Therefore, this solution does not have biological
relevance. $F_{7}$ has a biological relevance when $\rho_{1}\geq 28.3012$.
However, considering this range, the second element of $F_{7}$ is negative.
This way, the point $F_{7}$ is irrelevant.

The numerical solutions for Eqs. 4 (blue line), 5 (red line) and 6 (green
line) are displayed in Fig. 2(a) for $\rho_{1}=0.4$, $\alpha_{13}=1.5$,
$\rho_{2}=4.5$, $\alpha_{23}=0.2$, $\alpha_{32}=2.5$ and $\delta_{2}=0.5$. For
these parameters, the solution is a limit cycle. The solutions show that the
number of host cells decreases while the immune and cancer increase. After
that, the solution reaches an oscillatory behavior between the three
population cells. For different $\rho_{1}$ values, it is possible to observe
chaotic solutions. As shown in previous works [31], the dynamical system has a
chaotic attractor for $\rho_{1}=0.6$, as exhibited in Fig. 2(b). Fig. 2(c)
shows the 2-dimensional projection together with the fixed points.

Figure 2: (a) Time series for $x$, $y$ and $z$ in blue, red and green lines,
respectively. In panel (a), we consider $\rho_{1}=0.4$. (b) Phase space in
3-dimensional space and projection in the plane $y-z$ in the panel (c) for
$\rho_{1}=0.6$. We consider $\alpha_{13}=1.5$, $\rho_{2}=4.5$,
$\alpha_{23}=0.2$, $\delta_{2}=0.5$, $\alpha_{32}=2.5$, $x_{0}=0.80$,
$y_{0}=0.15$, and $z_{0}=0.05$.

In order to distinguish periodic dynamics from chaotic ones, we use the
recurrence plots (RPs), namely, recurrence quantification analysis (RQA),
which remain valid under fractional operators. The RPs, introduced by Eckmann
et al. [73], provide a visual representation of recurrences of the states of a
system in a $d$-dimensional phase space within a small deviation $\epsilon$,
as well as a quantitative analysis of the dynamical behavior exhibited by the
system. For this purpose, given a trajectory $\vec{x}(t)\in\mathbb{R}^{d}$, we
define a recurrence matrix R as

$R_{ij}=H(\epsilon-\|\vec{x}(t_{i})-\vec{x}(t_{j})\|),$ (15)

where $i,j=1,2,\ldots,N$, $N$ is the length of the time series, $H$ is the
Heaviside unit step function, $\epsilon$ is a small threshold and
$\|\vec{x}(t_{i})-\vec{x}(t_{j})\|$ is the spatial distance between two
states, $\vec{x}(t_{i})$ and $\vec{x}(t_{j})$, in phase space in terms of a
suitable norm.

The recurrence matrix is symmetric and binary with the elements $0$
representing the non-recurrent states and the elements $1$ representing the
recurrent ones. Two states are recurrent when the state at $t=t_{i}$ is close
(up to a distance $\epsilon$) to a different state at $t=t_{j}$. For the
threshold $\epsilon$, we choose it to be $\epsilon=0.01$. For a discussion on
the choice of $\epsilon$, see Appendix.

Several measures based on the RPs have been proposed [74, 75, 76, 77, 78]. The
most simple of them is the recurrence rate, RR, defined as

$\mathrm{RR}=\frac{1}{N^{2}}\sum_{i,j=1}^{N}R_{ij},$ (16)

which is a measure of the density of recurrent points. There are also measures
based on the diagonal lines formed in an RP, such as determinism (or
predictability)

$\mathrm{DET}=\frac{\sum_{\ell=\ell_{min}}^{N}\ell
P(\ell)}{\sum_{\ell=1}^{N}\ell P(\ell)},$ (17)

where $P(\ell)$ is the total number of diagonal lines with length $\ell$ and
minimal length $\ell_{min}$. Determinism is also a density as the recurrence
rate. It is the density of recurrent points that lie on a diagonal line and
refers to the degree to which the dynamics of a system is predictable and
repeatable over time. Systems with stochastic or chaotic behaviors cause no or
very short diagonal lines, whereas deterministic behavior exhibits longer
diagonal lines. Therefore, the determinism is expected to be small in the
first case and large in the second case. There are also measures based on the
vertical lines of length $v$, such as the laminarity (LAM) and the trapping
time (TT), for example. For a detailed discussion about these and other
measures, see Refs. [74, 75, 76, 77, 78] and references therein.

Entropy-based measures of RPs have been employed to characterize the dynamical
behavior of nonlinear systems [79, 80, 81, 82, 83, 84, 89]. The entropy of the
distribution of recurrence times (recurrence time entropy) has been utilized
as a tool for the detection of chaotic orbits. The vertical distances between
the diagonal lines, i.e., the gaps between them, are an estimate of the
recurrence times of the trajectory [85, 86, 87, 88]. The Shannon entropy of
the distribution of white vertical lines (the estimate of the recurrence
times) is defined as [79, 89, 84]

$\mathrm{RTE}=-\sum_{v=v_{\mathrm{min}}}^{v=v_{\mathrm{max}}}p_{w}(v)\ln{p_{w}(v)},$
(18)

where $v_{\mathrm{min}}$ ($v_{\mathrm{max}}$) is the length of the shortest
(longest) white vertical line. $p_{w}(v)=P_{w}(v)/N_{w}$ and $P_{w}(v)$ are
the relative distribution and the total number of white vertical lines of
length $v$, respectively, and $N_{w}$ is the total number of them. Due to the
finite size of a RP, the distribution of white vertical lines might be biased
by the white border lines, which are cut short by the border of the RP, thus
affecting the RQA measures, such as the RTE [90]. In order to avoid these
border effects, it is necessary to exclude the distribution of the lines that
begin and end at the border of the RP.

Figure 3: (a) Bifurcation diagram of $x_{\rm max}$ (black points) and $x_{\rm
min}$ (red points), (b) largest (blue line) and second largest Lyapunov (red
line) exponents, (c) recurrence rate, RR, (d) determinism, DET, and (e)
recurrence time entropy, RTE, as a function of $\rho_{1}$. We consider
$\alpha_{13}=1.5$, $\rho_{2}=4.5$, $\alpha_{23}=0.2$, $\delta_{2}=0.5$,
$\alpha_{32}=2.5$, and $\epsilon=0.01$.

To investigate the dynamic behavior of the cancer model, we consider
$\rho_{1}$ as the control parameter and compute the bifurcation diagram by
recording the local maxima, $x_{\rm max}$ (black dots) and local minima,
$x_{\rm min}$ (red dots), of $x(t)$, as shown in Fig. 3(a). For
$\rho_{1}<0.5$, the dynamics is periodic and becomes chaotic via period
doubling. The bifurcation diagram exhibits some periodic windows. We consider
the Lyapunov exponents ($\lambda$) as a classical method to compute chaotic
solutions. Figure 3(b) shows the largest, $\lambda_{1}$ and second largest,
$\lambda_{2}$, Lyapunov exponents by the blue and red lines, respectively. The
chaotic regimes are marked by at least one Lyapunov exponent greater than
zero. In the periodic windows, we observe $\lambda_{1}\approx 0$ and
$\lambda_{2}<0$. On the other hand, when the dynamics is chaotic,
$\lambda_{1}>0$. We also calculate the RQA measures RR, DET, and RTE in panels
(c), (d), and (e), respectively. To construct the recurrence matrix, we
consider the time series of $x_{\rm max}$, $\\{x_{\rm
max}^{(i)}\\}_{i=1,2,\ldots,N}$. We verify that periodic solutions exhibit
large values of RR and DET, with $\mathrm{DET}\rightarrow 1$, indicating that
all recurrent points lie on a diagonal line, whereas the RTE is close to zero
in these cases. RR and DET show smaller values during the chaotic windows,
while the RTE is larger. Therefore, there is a correlation between
$\lambda_{1}$ and these RQAs measures. To quantify this correlation, we use
the Pearson correlation coefficient,

$\rho_{x,y}=\frac{{\rm cov}(x,y)}{\sigma_{x}\sigma_{y}},$ (19)

where ${\rm cov}(x,y)$ is the covariance between the two data series and
$\sigma$ their respective standard deviation. We obtain for the correlation
coefficient between $\lambda_{1}$ and RR, DET, and RTE, the values
$\rho_{\lambda_{1},\mathrm{RR}}=-0.54$,
$\rho_{\lambda_{1},\mathrm{DET}}=-0.41$,
$\rho_{\lambda_{1},\mathrm{RTE}}=0.88$, respectively. Thus, the RTE is a great
alternative for the characterization of the dynamics of this system.

##  3 Fractional approach

The fractional extension of the model described by Eqs. (4), (5) and (6) is
obtained by making the following substitution: $Df\rightarrow D^{\alpha}f$,
where $D$ is the integer differential operator, i.e., $df/dt$ and
$\;{}_{0}D_{t}^{\alpha}f$ is the fractional differential operator and is
defined, in the Caputo sense, by

$\;{}_{0}D_{t}^{\alpha}f\equiv\frac{1}{\Gamma\left(1-\alpha\right)}\int_{0}^{t}dt^{\prime}\frac{1}{(t-t^{\prime})^{\alpha}}\frac{\partial}{\partial
t}f(\vec{r},t),$ (20)

where $\Gamma(\cdot)$ is the gamma function and $0<\alpha<1$ [41]. Therefore,
the extended cancer model is given by

$\;{}_{0}D_{t}^{\alpha}x$ $\displaystyle=$
$\displaystyle\rho_{1}x(1-x)-\alpha_{13}xz\;,$ (21) $\;{}_{0}D_{t}^{\alpha}y$
$\displaystyle=$
$\displaystyle\frac{\rho_{2}yz}{1+z}-\alpha_{23}yz-\delta_{2}y\;,$ (22)
$\;{}_{0}D_{t}^{\alpha}z$ $\displaystyle=$ $\displaystyle
z(1-z)-xz-\alpha_{32}yz\;.$ (23)

The time series for $\alpha=0.995$ (blue line), $\alpha=0.98$ (red line) and
$\alpha=0.95$ (green line) are displayed in Fig. 4. The numerical integration
is made by the algorithm described in Ref. [95]. The black dotted lines are
for the standard case, as exhibited in Fig. 2. These results show that the
fractional order attenuates the oscillations until a fixed point as $\alpha$
decreases. For instance, for $\alpha=0.95$, the fixed point is equal to
$(x,y,z)=(0.5032,\ 0.1456,\ 0.1325)$, which corresponds to $F_{5}$ calculated
for $\rho_{1}=0.4$.

Figure 4: Time series for $x$ in the panels (a) and (d), for $y$ in the panels
(b) and (e), and for $z$ in panels (c) and (f). The panels (a-c) are for
$\rho_{1}=0.4$ and (d-f) for $\rho_{1}=0.6$. The black dotted lines are for
$\alpha=1.0$, the blue line is for $\alpha=0.995$, the red line is for
$\alpha=0.98$ and the green line is for $\alpha=0.95$. We consider
$\alpha_{13}=1.5$, $\rho_{2}=4.5$, $\alpha_{23}=0.2$, $\delta_{2}=0.5$,
$\alpha_{32}=2.5$, $x_{0}=0.80$, $y_{0}=0.15$, and $z_{0}=0.05$.

The oscillatory behavior is attenuated in the chaotic regime for
$\rho_{1}=0.6$ and the dynamics is transited to periodic behavior. By
computing the RTE as a function of $\alpha$, we see that there are some
periodic windows. For $\alpha\lessapprox 0.9966$, the dynamics become periodic
(Fig. 5(a)). The RTE information agrees with the bifurcation diagram, as shown
in Fig. 5(b). The red and black points correspond to the $x$ minimum and
maximum values. The attractors in Fig. 6 display that the periodic behavior is
a limit cycle that converges to a fixed point as $\alpha$ approaches to $0.9$.
The green line is for $\alpha=0.995$, the red line for $\alpha=0.975$, the
blue line for $\alpha=0.936$, the orange line for $\alpha=0.915$ and the black
dotted point for $\alpha=0.885$.

Figure 5: (a) RTE and (b) bifurcation diagram as a function of $\alpha$. We
consider $\alpha_{13}=1.5$, $\rho_{1}=0.6$, $\rho_{2}=4.5$, $\alpha_{23}=0.2$,
$\delta_{2}=0.5$, $\alpha_{32}=2.5$, $x_{0}=0.80$, $y_{0}=0.15$, $z_{0}=0.05$,
and $\epsilon=0.01$. Figure 6: Phase portrait for different $\alpha$ values.
The green curve is for $\alpha=0.995$, the red curve for $\alpha=0.975$, the
blue curve for $\alpha=0.936$, the orange curve for $\alpha=0.915$ and the
black dot for $\alpha=0.885$. We consider $\rho_{1}=0.6$, $\alpha_{13}=1.5$,
$\rho_{2}=4.5$, $\alpha_{23}=0.2$, $\delta_{2}=0.5$, and $\alpha_{32}=2.5$.

One way to differentiate fixed points and limit cycles is by computing the
mean standard deviation ($\sigma$). To do that, we discard a transient time
($\tau$). Considering this time series, we compute $\sigma$. Figure 7(a)
displays the parameter plane $\rho_{1}\times\alpha$ as a function of $\sigma$
in color scale. The black points indicate when the solution is a fixed point
($\sigma<10^{-3}$) and the colorful points indicate limit cycle solutions.
Although some black points are mixed with colorful ones, a relation between
$\alpha$ and $\rho_{1}$ separates most parts of the limit cycle from the fixed
point solutions. This curve is indicated by the white dotted line in the panel
(a) and is given by
$\alpha=\beta_{1}e^{\beta_{2}\rho_{1}}+\beta_{3}e^{\beta_{4}\rho_{1}}$, where
$\beta_{1}=3.150$, $\beta_{2}=-7.800$, $\beta_{3}=0.808$ and
$\beta_{4}=0.110$. Figure 7(b) exhibits the same plane parameter with the time
to the solution for one initial condition reaches the steady solution in the
color scale. The black dotted line indicates a region of super transient,
which is defined by
$\alpha=\xi_{1}e^{\xi_{2}\rho_{1}}+\xi_{3}e^{\xi_{4}\rho_{1}}$, where
$\xi_{1}=3.25$, $\xi_{2}=-7.7$, $\xi_{3}=0.802$, and $\xi_{4}=0.108$. Note
that $\beta_{i}\approx\xi_{i}$. In this way, the transition from the limit
cycle to fixed points crosses a supper transient region.

Figure 7: $\rho_{1}\times\alpha$ as a function of $\sigma$ (mean standard
deviation) in the panel (a) and $\tau$ (average transient time) in the panel
(b). We consider $\alpha_{13}=1.5$, $\rho_{2}=4.5$, $\alpha_{23}=0.2$,
$\delta_{2}=0.5$, and $\alpha_{32}=2.5$.

As our results suggest, the transition from a limit cycle to fixed points is
associated with a supper transient time. In this way, we select $\rho_{1}=0.6$
and compute the mean transient time $\langle\tau\rangle$ for 100 different
initial conditions close to the critical point. The critical point defines the
transition and equals $\alpha_{c}=0.899$. Figure 8 shows the mean transient
time ($\langle\tau\rangle$) as a function of $|\alpha-\alpha_{c}|$ by the
black points. The average transient lifetime scales with $\alpha$ as a power-
law according to

$\langle\tau\rangle\propto|\alpha-\alpha_{c}|^{\gamma},$ (24)

and has two slopes. By fitting the data (red line in Fig. 8), we obtain an
exponent of $\gamma\approx-3/2$, with a correlation coefficient of 0.98 and
$\gamma\approx-1/3$ (green line) with a correlation coefficient of 0.99. All
the points evolve to $F_{5}$.

Figure 8: $\langle\tau\rangle$ versus $|\alpha-\alpha_{c}|$, where
$\alpha_{c}=0.899$. The curve has two slopes, in the red curve
$\gamma\approx-3/2$ and in the green $\gamma\approx-1/3$. The respective
correlation coefficients are 0.98 and 0.99. We consider $\rho_{1}=0.6$,
$\alpha_{13}=1.5$, $\rho_{2}=4.5$, $\alpha_{23}=0.2$, $\delta_{2}=0.5$, and
$\alpha_{32}=2.5$.

##  4 Conclusions

In this work, we consider a fractional extension of a cancer model. This model
describes the interactions among host, effector immune, and tumor cells. Due
to the non-linearity, the dynamical behavior exhibits chaotic solutions for
some parameter ranges. For the equations that govern the model, we obtain
analytical solutions for the fixed points and discuss their stability.

In order to study a more broad case, we consider numerical solutions. As
previously reported works, we verify chaotic solutions for $\rho_{1}=0.6$.
Furthermore, we construct bifurcation diagrams considering the maxima and
minimum points of the host cells. The dynamical behavior exhibited by the
bifurcation diagram also was verified by the Lyapunov exponents and recurrence
quantification analysis (RQA) measures. With regard to the RQA measures, we
compute the recurrence rate (RR), determinism (DET), and recurrence time
entropy (RTE). All three of them agree with the results from the Lyapunov
exponents. By computing the Pearson correlation coefficient between them and
$\lambda_{1}$, we observe a higher coefficient for the RTE. Thus, we consider
only the RTE to analyze the effects of fractional order.

We study the effects of fractional operators in the cancer model by the Caputo
definition. Our main goal is to understand what happens in the chaotic
behavior. By computing the RTE, our results suggest that the dynamical system
becomes periodic for $\alpha\lessapprox 0.9966$. When we look for the phase
portrait, the results show that the solutions transit from a limit cycle to a
fixed point, which is equal to the fixed points calculated analytically. The
limit cycle regime is in the range $\alpha\in(0.9966,0.899)$ for
$\rho_{1}=0.6$. For different $\rho_{1}$ values, we construct the parameter
space $\rho_{1}\times\alpha$ as a function of the mean standard deviation
$\sigma$. Values of $\sigma$ close to zero show fixed point solutions. Our
results show the existence of an exponential curve that delimits fixed point
solutions from limit cycles in this parameter space. Also, we analyze the
transient time related to the transition from the limit cycle to a fixed
point. For $\rho_{1}=0.6$, we find that for $\alpha<0.899$ the dynamics go to
the fixed point $F_{5}$ by a super transient with order $10^{4}$. The super
transient curve follows a linear decay according to two slopes, equal to
$-3/2$ and $-1/3$.

## Acknowledgments

This work was possible with partial financial support from the following
Brazilian government agencies: CNPq, CAPES, Fundação Araucária and São Paulo
Research Foundation (FAPESP 2018/03211-6, 2022/13761-9). R. L. V. received
partial financial support from the following Brazilian government agencies:
CNPq (403120/2021-7, 301019/2019-3), CAPES (88881.143103/2017-01), FAPESP
(2022/04251-7). E. K. L. received partial financial support from the following
Brazilian government agency: CNPq (301715/2022-0). E.C.G. received partial
financial support from Coordenação de Aperfeiçoamento de Pessoal de Nível
Superior - Brasil (CAPES) \- Finance Code 88881.846051/2023-01. We would like
to thank the 105 Group Science (www.105groupscience.com).

## Declarations

Conflict of Interest:
The author declares that there exists no competing financial interest or
personal relationships that could have appeared to influence the work reported
in this paper.

## Appendix

In Section 2, we introduced the cancer model and analyzed its dynamics for
different parameter values by means of a bifurcation diagram, the Lyapunov
exponents, and three RQA measures: RR, DET, and RTE. We chose the threshold
$\epsilon$ to be $\epsilon=0.01$. Many researchers have addressed the problem
of finding the best value of $\epsilon$ [91, 92, 93, 94]. Here, we used the
correlation coefficient, Eq. (19), to determine $\epsilon$. We computed the
RTE as a function of $\rho_{1}$ for different threshold values, ranging from
$10^{-5}$ to $10^{-1}$, and calculated the correlation coefficient between
$\lambda_{1}$ and RTE (Fig. 9). Even if we choose a very small $\epsilon$
($10^{-5}$) or a relatively large $\epsilon$ ($10^{-1}$), we still obtain a
rather high correlation coefficient, with the highest value being around
$\epsilon\approx 10^{-2}$, and hence our choice of $\epsilon$. Furthermore,
the fact that $\rho_{\lambda_{1},\mathrm{RTE}}$ does not change too much
within this range of $\epsilon$ indicates that, in our case, the choice of
$\epsilon$ is not as sensible as it would be in other cases. In fact, there is
a wide range of values for $\epsilon$ in which the results are good.

Figure 9: The correlation coefficient between $\lambda_{1}$ and RTE as a
function of the threshold $\epsilon$.

## References

  * [1] Hausman DM. What is cancer? Perspectives in Biology and Medicine 2019;62(4):778-784.
  * [2] Ferlay J, Colombet M, Soerjomataram I, Parkin DM, Piñeros M, Znaor A, Bray F. Cancer statistics for the year 2020: An overview. Cancer Epidemiology 2021;149:778-789.
  * [3] Siegel RL, Miller KD, Fuchs HE, Jemal A. Cancer statistics, 2022. CA Cancer J Clin 2022;72:7-33.
  * [4] Siegel RL, Miller KD, Jemal A. Cancer statistics, 2018. CA Cancer J Clin 2018;68:7–30.
  * [5] Jones PA, Baylin SB. The Epigenomics of Cancer. Cell 2007;128:683-692.
  * [6] Mufudza C, Sorofa W, Chiyaka ET. Assessing the Effects of Estrogen on the Dynamics of Breast Cancer. Computational and Mathematical Methods in Medicine 2012;473572.
  * [7] Waldman AD, Fritz JM, Lenardo MJ. A guide to cancer immunotherapy: from T cell basic science to clinical practice. Nature Reviews Immunology 2020;20:651-668.
  * [8] López AG, Seoane JM, Sanjuán MAF. Modelling Cancer Dynamics Using Cellular Automata. In: Berezovskaya, F., Toni, B. (eds) Advanced Mathematical Methods in Biosciences and Applications 2019;159-205.
  * [9] Iarosz KC, Borges FS, Batista AM, Baptista MS, Siqueira RAN, Viana RL, Lopes SR. Mathematical model of brain tumour with glia–neuron interactions and chemotherapy treatment. Journal of Theoretical Biology 2015;368;113-121.
  * [10] Trobia J, Gabrick EC, Seifert EG, Borges FS, Protachevicz PR, Szezech Jr. JD, Iarosz KC, Santos MS, Caldas IL, Tian K, Ren HP, Grebogi C, Batista AM. Effects of drug resistance in the tumour-immune system with chemotherapy treatment. Indian Academy of Sciences Conference Series 2020;3(1):39-44.
  * [11] Trobia J, Tian K, Batista AM, Grebogi C, Pen HP, Santos MS, Protachevicz PR, Borges FS, Szezech Jr JD, Viana RL, Caldas IL, Iarosz KC. Mathematical model of brain tumour growth with drug resistance. Communications in Nonlinear Science and Numerical Simulation 2021;103:106013.
  * [12] Castiglione F, Piccoli B. Cancer immunotherapy, mathematical modeling and optimal control. Journal of Theoretical Biology 2007;247(4):723-732.
  * [13] López AG, Iarosz KC, Batista AM, Seoane JM, Viana RL, Sanjuán MAF. The role of dose density in combination cancer chemotherapy. Communications in Nonlinear Science and Numerical Simulation 2019;79:104918.
  * [14] Borges FS, Iarosz KC, Ren HP, Batista AM, Baptista MS, Viana RL, Lopes SR, Grebogi C. Model for tumour growth with treatment by continuous and pulsed chemotherapy. BioSystems 2014;116:43-48.
  * [15] Liu Z, Yang C. A Mathematical Model of Cancer Treatment by Radiotherapy. Computational and Mathematical Methods in Medicine 2014;2014:172923.
  * [16] Mamat M, Subiyanto, Kartono A. Mathematical Model of Cancer Treatments Using Immunotherapy, Chemotherapy and Biochemotherapy. Applied Mathematical Sciences 2013;7(5):247-261.
  * [17] Sayari E, da Silva ST, Iarosz KC, Viana RL, Szezech Jr JD, Batista AM. Prediction of fluctuations in a chaotic cancer model using machine learning. Chaos, Solitons and Fractals 2022;164:112616.
  * [18] de Pillis LG, Radunskaya A. The Dynamics of an Optimally Controlled Tumor Model: A Case Study. Mathematical and Computer Modelling 2003;37:1221-1244.
  * [19] Iarosz KC, Martins CC, Batista AM, Viana RL, Lopes SR, Caldas IL, Penna TJP. On a cellular automaton with time delay for modelling cancer tumors. Journal of Physics: Conference Series 2011;285:012015.
  * [20] Kuznetsov VA, Makalkin IA, Taylor MA, Perelson AS. Nonlinear dynamics of immunogenic tumors: Parameter estimation and global bifurcation analysis. Bulletin of Mathematical Biology 1994;56(2):295-321.
  * [21] Pinho STR, Freedman HI, Nani F. A chemotherapy model for the treatment of cancer with metastasis. Mathematical and Computer Modelling: An International Journal 2002;36(7-8):773–803.
  * [22] Amatruda JF, Shepard JL, Stern HW, Zon LI. Zebrafish as a cancer model system. Cancer Cell 2002;1(3):229-231.
  * [23] Altrock PM, Liu LL, Michor F. The mathematics of cancer: integrating quantitative models. Nature Reviews Cancer 2015;15:730-745.
  * [24] Tuveson D, Clevers H. Cancer modeling meets human organoid technology. Science 2019;364(6444):952-955.
  * [25] Dehingia K, Sarmah HK, Alharbi Y, Hosseini K. Mathematical analysis of a cancer model with time-delay in tumor-immune interaction and stimulation processes. Advances in Difference Equations 2021;2021:473.
  * [26] Díaz-Marín H, López-Hernández JF, Osuna O. Stability and eradication of tumor in a model with almost periodically radiated cells. Journal of Applied Mathematics and Computing 2022;68:3781–3797.
  * [27] López AG, Iarosz KC, Batista AM, Seoane JM, Viana RL, Sanjuán MAF. Nonlinear cancer chemotherapy: Modelling the Norton-Simon hypothesis. Communications in Nonlinear Science and Numerical Simulation 2019;70:307-317.
  * [28] López AG, Seoane JM, Sanjuán MAF. A Validated Mathematical Model of Tumor Growth Including Tumor–Host Interaction, Cell-Mediated Immune Response and Chemotherapy. Bulletin of Mathematical Biology 2014;76:2884-2906.
  * [29] Tél T, Gruiz M. Chaotic dynamics: An introduction based on Classical Mechanics. Cambridge University Press 2006.
  * [30] Itik M, Banks SP. Chaos in a three-dimensional cancer model. International Journal of Bifurcation and Chaos 2010;20(1):71-79.
  * [31] Letellier C, Denis F, Aguirre LA. What can be learned from a chaotic cancer model? Journal of Theoretical Biology 2013;322:7–16.
  * [32] Khajanchi S, Perc M, Gosh D. The influence of time delay in a chaotic cancer model. Chaos: An Interdisciplinary Journal of Nonlinear Science 2018;28:103101.
  * [33] Abernethy S, Gooding RJ. The importance of chaotic attractors in modelling tumour growth. Physica A: Statistical Mechanics and its Applications 2018;507:268-277.
  * [34] Gallas MR, Gallas MR, Gallas JAC. Distribution of chaos and periodic spikes in a three-cell population model of cancer. The European Physical Journal Special Topics 2014;223:2131-2144.
  * [35] Khajanchi S. Bifurcation analysis of a delayed mathematical model for tumor growth. Chaos, Solitons and Fractals 2015;77:264-276.
  * [36] Kemwoue FF, Dongo JM, Mballa RN, Gninzanlong CL, Kemayou MW, Mokhtari B, Biya-Motto F, Atangana J. Bifurcation, multistability in the dynamics of tumor growth and electronic simulations by the use of Pspice. Chaos, Solitons and Fractals 2020;134:109689.
  * [37] Li J, Xie X, Chen Y, Zhang D. Complex dynamics of a tumor-immune system with antigenicity. Applied Mathematics and Computation 2021;400:126052.
  * [38] Valle PA, Coria LN, Gamboa D, Plata C. Bounding the Dynamics of a Chaotic-Cancer Mathematical Model. Mathematical Problems in Engineering 2018:9787015.
  * [39] Uthamacumaran A. A review of dynamical systems approaches for the detection of chaotic attractors in cancer networks. Patterns 2021;2(4):100226.
  * [40] Evangelista LR, Lenzi EK. An Introduction to Anomalous Diffusion and Relaxation, Springer Nature, 2023.
  * [41] Evangelista LR, Lenzi EK. Fractional Diffusion Equations and Anomalous Diffusion, Cambridge University Press, 2018.
  * [42] Somer A, Galovic S, Lenzi EK, Novatski A, Djordjevic K. Temperature profile and thermal piston component of photoacoustic response calculated by the fractional dual-phase-lag heat conduction theory. International Journal of Heat and Mass Transfer 2023;203:123801.
  * [43] Somer A, Popovic MN, da Cruz GK, Novatski A, Lenzi EK, Galovic SP. Anomalous thermal diffusion in two-layer system: The temperature profile and photoacoustic signal for rear light incidence. International Journal of Thermal Sciences 2022;179:107661.
  * [44] do Carmo WP, Santos AF, Lenzi MK, Fortuny M, Lenzi EK. A new fractional model applied to description of the viscoelastic creep behavior of two Brazilian oils and their w/o emulsions. Digital Chemical Engineering 2023;6:100069.
  * [45] Cius D, Menon Jr L, dos Santos MAF, de Castro ASM, Andrade FM. Unitary evolution for a two-level quantum system in fractional-time scenario. Physical Review E 2022:106;054126.
  * [46] Lenzi EK, Ribeiro HV, dos Santos MAF, Rossato R, Mendes RS. Time dependent solutions for a fractional Schrödinger equation with delta potentials. Journal of Mathematical Physics 2013:54(8);082107.
  * [47] Gabrick EC, Sayari E, de Castro ASM, Trobia J, Batista AM, Lenzi EK. Fractional Schrödinger equation and time dependent potentials. Communications in Nonlinear Science and Numerical Simulation 2023;107275.
  * [48] Srivastava HM, Saad KM, Khader MM. An efficient spectral collocation method for the dynamic simulation of the fractional epidemiological model of the Ebola virus. Chaos, Solitons and Fractals 2020:140;110174.
  * [49] Dong NP, Long HV, Khastan A. Optimal control of a fractional order model for granular SEIR epidemic with uncertainty. Communications in Nonlinear Science and Numerical Simulation 2020:88;105312.
  * [50] Kumar A, Kumar S. A study on eco-epidemiological model with fractional operators. Chaos, Solitons and Fractals 2022:156;111697.
  * [51] Mahmoud EE, Trikha P, Jahanzaib LS, Almaghrabi OA. Dynamical analysis and chaos control of the fractional chaotic ecological model. Chaos, Solitons and Fractals 2020:141;110348.
  * [52] Coccolo M, Seoane JM, Lenci S, Sanjuán MAF. Fractional damping effects on the transient dynamics of the Duffing oscillator. Communications in Nonlinear Science and Numerical Simulation 2023:117;106959.
  * [53] Nadeem M, Jafari H, Akgül A, De la Sen M. A Computational Scheme for the Numerical Results of Time-Fractional Degasperis–Procesi and Camassa–Holm Models. Symmetry 2022:14(12);2532.
  * [54] Vivekanandhan G, Abdolmohammadi HR, Natiq H, Rajagopal K, Jafari S, Namazi H. Dynamic analysis of the discrete fractional-order Rulkov neuron map. Mathematical Biosciences and Engineering 2023:20(3);4760-4781.
  * [55] Sun HG, Zhang Y, Baleanu D, Chen W, Chen Y. A new collection of real world applications of fractional calculus in science and engineering. Communications in Nonlinear Science and Numerical Simulation 2018:64;213-231.
  * [56] Atangana A, Jain S. The role of power decay, exponential decay and Mittag-Leffler function’s waiting time distribution: Application of cancer spread. Physica A 2018;512:330-351.
  * [57] Caputo M, Fabrizio M. On the Singular Kernels for Fractional Derivatives. Some Applications to Partial Differential Equations. Progress in Fractional Differentiation and Application 2021;7(2):79-82.
  * [58] Caputo M, Fabrizio M. A new Definition of Fractional Derivative without Singular Kernel. Progress in Fractional Differentiation and Applications 2015;1(2):73-85.
  * [59] Atangana A, Baleanu D. New fractional derivatives with non-local and non-singular kernel: Theory and application to heat transfer model. Thermal Science 2016;20(2):763-769.
  * [60] Ghanbari B. On the modeling of the interaction between tumor growth and the immune system using some new fractional and fractional-fractal operators. Advances in Difference Equations 2020;2020:585.
  * [61] Naik PA, Zu J, Naik M. Stability analysis of a fractional-order cancer model with chaotic dynamics. International Journal of Biomathematics 2021;14(6):2150046.
  * [62] Xuan L, Ahmad S, Ullah A, Saifullah S, Akgül A, Qu H. Bifurcations, stability analysis and complex dynamics of Caputo fractal-fractional cancer model. Chaos, Solitons and Fractals 2022;159:112113.
  * [63] Ahmad S, Ullah A, Abdeljawad T, Akgül A, Mlaiki N. Analysis of fractal-fractional model of tumor-immune interaction. Results in Physics 2021;25:104178.
  * [64] Ahmed E, Hashis AH, Rihan FA. On fractional order cancer model. Journal of Fractional Calculus and Applications 2012;3(2):1-6.
  * [65] Arfan M, Shah K, Ullah A, Shutaywi M, Kumam P, Shah Z. On fractional order model of tumor dynamics with drug interventions under nonlocal fractional derivative. Results in Physics 2021;21:103783.
  * [66] Baleanu D, Jajarmi A, Sajjadi SS, Mozyrska D. A new fractional model and optimal control of a tumor-immune surveillance with non-singular derivative operator. Chaos 2019;29:083127.
  * [67] Debbouche N, Ouannas A, Grassi G, Al-Hussein ABA, Tahir FR, Saad KM, Jahanshahi H, Aly AA. Chaos in Cancer Tumor Growth Model with Commensurate and Incommensurate Fractional-Order Derivatives. Computational and Mathematical Methods in Medicine 2022;5227503.
  * [68] Hassani H, Machado JAT, Avazzadeh Z, Safari E, Mehrabi S. Optimal solution of the fractional order breast cancer competition model. Scientific Reports 2021;11:15622.
  * [69] Iyiola OS, Zaman FD. A fractional difusion equation model for cancer tumor. AIP Advances 2014;4:107121.
  * [70] Kumar S, Ajay Kumar, Samet B, Gómez-Aguilar JF, Osman MS. A chaos study of tumor and effector cells in fractional tumor-immune model for cancer treatment. Chaos, Solitons and Fractals 2020;141:110321.
  * [71] Rehaman MAU, Ahmad J, Hassan A, Awrejcewicz J, Pawlowski W, Karamti H, Alharbi FM. The Dynamics of a Fractional-Order Mathematical Model of Cancer Tumor Disease. Symmetry 2022;14:1694.
  * [72] Solís-Pérez JE, Gómez-Aguilar JF, Atangana A. A fractional mathematical model of breast cancer competition model. Chaos, Solitons and Fractals 2019;127:38-54.
  * [73] Eckmann JP, Kamphorst SO, Ruelle D. Recurrence Plots of Dynamical Systems. Europhysics Letters 1987;4(9):973.
  * [74] Marwan N, Wessel N, Meyerfeldt U, Schirdewan A, Kurths J. Recurrence-plot-based measures of complexity and their application to heart-rate-variability data. Phys. Rev. E 2002;66(2):026702.
  * [75] Marwan N, Kurths. Nonlinear analysis of bivariate data with cross recurrence plots. Physics Letters A 2002;302(5):299-307.
  * [76] Marwan N. A historical review of recurrence plots. The European Physical Journal Special Topics 2008;164(1):3-12.
  * [77] Bedartha G. A Brief Introduction to Nonlinear Time Series Analysis and Recurrence Plots. Vibration 2019;2(4):332-368.
  * [78] Marwan N, Webber CL. Mathematical and Computational Foundations of Recurrence Quantifications. Springer International Publishing 2015.
  * [79] Little MA, McSharry PE, Roberts SJ, Costello DAE, Moroz IM. Exploiting Nonlinear Recurrence and Fractal Scaling Properties for Voice Disorder Detection. BioMedical Engineering OnLine 2007;6(1):23
  * [80] Webber CL, Zbilut JP. Dynamical assessment of physiological systems and states using recurrence plot strategies. Journal of Applied Physiology 1994;76(2):965-973.
  * [81] Trulla LL, Giuliani A, Zbilut JP, Webber CL. Recurrence quantification analysis of the logistic equation with transients. Physics Letters A 1996;223(4):255-260.
  * [82] Ngamga EJ, Senthilkumar DV, Kurths J. Dynamics between order and chaos revisited. The European Physical Journal Special Topics 2010;191(1):15-27.
  * [83] Eroglu D, Peron TKDM, Marwan N, Rodrigues FA, Costa LF, Sebek M, Kiss IZ, Kurths J. Entropy of weighted recurrence plots. Physical Review E 2014;90(4):042919.
  * [84] Sales MR, Mugnaine M, Szezech Jr. D, Viana RL, Caldas IL, Marwan N, Kurths J. Stickiness and recurrence plots: An entropy-based approach. Chaos: An Interdisciplinary Journal of Nonlinear Science 2023;33(3):033140.
  * [85] Zou Y, Thiel M, Romano MC, Kurths J. Characterization of stickiness by means of recurrence. Chaos: An Interdisciplinary Journal of Nonlinear Science 2007;17(4):
  * [86] Zou Y, Pazó D, Romano MC, Thiel M, Kurths J. Distinguishing quasiperiodic dynamics from chaos in short-time series. Phys. Rev. E 2007;76(1):016210.
  * [87] Baptista MS, Ngama EJ, Pinto PRF, Brito M, Kurths J. Kolmogorov–Sinai entropy from recurrence times. Physics Letters A 2010;374(9):1135-1140.
  * [88] Ngamga EJ, Senthilkumar DV, Prasad A, Parmananda P, Marwan N, Kurths J. Distinguishing dynamics using recurrence-time statistics. Phys. Rev. E 2012;85 (2):026217.
  * [89] Kraemer KH, Donner RV, Heitzig J, Marwan N. Recurrence threshold selection for obtaining robust recurrence characteristics in different embedding dimensions. Chaos: An Interdisciplinary Journal of Nonlinear Science 2018;28(8):085720.
  * [90] Kraemer KH, Marwan N. Border effect corrections for diagonal line based recurrence quantification analysis measures. Physics Letters A 2019;383(34):125977.
  * [91] Zbilut JP, Zaldivar-Comenges JM, Strozzi F. Recurrence quantification based Liapunov exponents for monitoring divergence in experimental data. Physics Letters A 2002;297(3):173-181.
  * [92] Thiel M, Romano MC, Kurths J, Meucci R, Allaria E, Arecchi FT. Influence of observational noise on the recurrence quantification analysis. Physica D: Nonlinear Phenomena 2002;171(3):138-152.
  * [93] Schinkel S, Dimigen O, Marwan N. Selection of recurrence threshold for signal detection. The European Physical Journal Special Topics 2008;164(1):45-53.
  * [94] Medrano J, Kheddar A, Annick L, Sofiane R. Radius selection using kernel density estimation for the computation of nonlinear measures. Chaos: An Interdisciplinary Journal of Nonlinear Science 2021;31(8):083131.
  * [95] Garrappa R. Numerical Solution of Fractional Differential Equations: A Survey and a Software Tutorial. Mathematics 2018;6(16):1-23.


# Using Machine Learning and Natural Language Processing Techniques to Analyze
and Support Moderation of Student Book Discussions

Jernej Vivod
Faculty of Computer and Information Science, University of Ljubljana
Ljubljana, 1000, Slovenia
<EMAIL_ADDRESS>

###### Abstract

The increasing adoption of technology to augment or even replace traditional
face-to-face learning has led to the development of a myriad of tools and
platforms aimed at engaging the students and facilitating the teacher’s
ability to present new information. The IMapBook project aims at improving the
literacy and reading comprehension skills of elementary school-aged children
by presenting them with interactive e-books and letting them take part in
moderated book discussions. This study aims to develop and illustrate a
machine learning-based approach to message classification that could be used
to automatically notify the discussion moderator of a possible need for an
intervention and also to collect other useful information about the ongoing
discussion. We aim to predict whether a message posted in the discussion is
relevant to the discussed book, whether the message is a statement, a
question, or an answer, and in which broad category it can be classified. We
incrementally enrich our used feature subsets and compare them using standard
classification algorithms as well as the novel Feature stacking method. We use
standard classification performance metrics as well as the Bayesian correlated
t-test to show that the use of described methods in discussion moderation is
feasible. Moving forward, we seek to attain better performance by focusing on
extracting more of the significant information found in the strong temporal
interdependence of the messages.

##  1 Introduction

Recent decades have brought about an increase in the use of computer-based
tools in practically every field of human endeavor. The field of education is
no exception. Such tools can be used to augment or even completely replace
traditional face-to-face teaching methods. The emergence of online learning
platforms has necessitated the development of means to enable learning
activities, such as group discussions, to be performed through the use of
technology. One such example of a learning platform is the IMapBook software
suite aimed at increasing the literacy and reading comprehension skills of
elementary school-aged children through the use of web-based eBooks, embedded
games related to their contents, as well as moderated group discussions.
Keeping these discussions constructive and relevant can be difficult and
usually requires a discussion moderator to be present at all times. This can
limit the opportunities for such discussions to take place. Leveraging the
methods and insights from the fields of artificial intelligence and machine
learning, we can attempt to develop systems to automatically classify messages
into different categories and detect when the discussion has veered off course
and necessitates intervention. Our research tackles this problem using a
compilation of discussions obtained during pilot studies testing the
effectiveness of using the IMapBook software suite in 4th-grade classrooms.
The studies were performed in 8 different Slovene primary schools and, in
total, included 342 students. The discussions consist of 3541 messages along
with annotations specifying their relevance to the book discussion, type,
category, and broad category. The ID of the book being discussed and the time
of posting are also included, as are the poster’s school, cohort, user ID, and
username. Each message was also manually translated into English to aid non-
Slovene-speaking researchers. The use of the Slovene language presents unique
challenges in applying standard language processing methods, many of which are
not as readily available as for other, more widely spoken languages.

Given a sequence of one or more newly observed messages, we want to estimate
the relevance of each message to the actual topic of discussion. Namely, we
want to assign messages into two categories — relevant to the book being
discussed or not. Additionally, we want to predict whether the message is a
question, an answer, or a statement which we call the type of the message.
Finally, we want to assign a category label to each message where the possible
labels can be either ’chatting’, ’switching’, ’discussion’, ’moderating’, or
’identity’. Building a predictive model capable of performing such predictions
with acceptable performance would allow us to experiment with including this
new level of automation in the IMapBook software suite as well as in any
related products. The research insights are also applicable to areas such as
online user comments and content moderation.

##  2 Related Work

The objective of our research is closely related to tasks concerning online
content moderation which has been the subject of much research in recent
years. Perhaps one of the earliest studies done on this subject is the 2009
study by Yin et al. Yin et al. (2009) in which the authors used
sentiment/contextual features in tandem with the TF-IDF approach to detect
online harassment. An earlier study by Mclaren et al. specifically focuses on
the use of machine learning techniques to support the mediation of student
online discussions in a uniquely constrained network-like environment that
differs significantly from ours Mclaren et al. (2007). A 2016 study by Kadunc
focuses on using machine learning methods to analyze the sentiment of Slovene
online comments and provides an important contribution in the form of an
opinion lexicon Kadunc (2016). However, the specifics and unique challenges
presented by the problem of classifying short Slovene text produced by this
age group remains an area with little to no research currently done.

##  3 Methods

Achieving the goal of creating a working predictive model for the task of
message classification requires careful processing of the raw data in such a
way as to expose as much useful information as possible. This process of
feature extraction and feature engineering often results in very high
dimensional descriptions of our data that can be prone to problems arising as
part of the so-called curse of dimensionality Domingos (2012). This can be
mitigated by using classification models well-suited for such data as well as
performing feature ranking and feature selection.

###  3.1 Extracting and Engineering Features

Building a quality predictive model requires a good characterization of each
message in terms of discriminative and non-redundant features. Extracting such
features from raw text data is a non-trivial task that is subject to much
research in the field of natural language processing. Here, we describe the
feature extraction process used in our study. A detailed evaluation of the
features is presented in the Results section.

####  3.1.1 General Message Features

Looking at the messages in our dataset, we can immediately notice simple but
potentially important differences between the messages in terms of word count,
punctuation use, and other attributes that can be easily deduced by merely
inspecting the raw data without any need for context.

We extract the number of words in the message, the maximal, minimal, and
average word lengths, the number of digits in the message, the number of
punctuation marks in the text, the number of capital letters in the text, the
number of consecutive repeated characters and check whether the message starts
with a capital letter and if it ends with a period.

####  3.1.2 Important Words

We can gain valuable insight by observing the presence of members from
important word groups. We compiled lists of chat usernames used in the
discussions, common given names in Slovenia, common curse words used in
Slovenia as well as any proper names found in the discussed books. We also
created a list of nouns, verbs and adjectives which we observed to be highly
discriminative such as misliti (to think), knjiga (book), najljubše (favorite)
among others. We noted any presence or absence of words from these lists as
features describing the pertaining message.

####  3.1.3 Multiset-Based Features

The Bag-of-words model and its variations characterize documents by counting
the occurrences of each word from a pre-defined set. It is based on the simple
assumption that similar documents share a similar vocabulary. We can augment
the basic Bag-of-words model by weighing each word in the vocabulary
proportional to its rarity in the corpus, using the assumption that rarity
implies discriminativity.

We constructed a simple bag-of-words model using unigrams and bigrams. We
required both the unigrams and bigrams to appear in at least two messages. All
words were converted to lower case, any punctuation and emojis were
disregarded and the words were lemmatized before the construction of the
model. We also converted any consecutive repeated letters in the words to
single occurrences. The resulting model consisted of 2009 unique unigrams and
bigrams. We used a similar method to include part-of-speech tagging into our
set of considered features. We describe this process in detail in the
following subsection.

####  3.1.4 Part-of-Speech Tagging-Based Features

Part-of-speech tagging is the process of labeling the words in a text based on
their corresponding part of speech. Describing each message in terms of its
associated part-of-speech labels allows us to use another perspective from
which we can view and analyze the corpus. The use of non-standard Slovene and
misspellings make part-of-speech tagging a non-trivial task. We attempted to
solve the issue by constructing a simple dictionary mapping known non-standard
and colloquial versions of Slovene words into their standard equivalents. We
constructed the dictionary using the corpus available as part of the JANES
project Fišer (2020). Again, we removed any punctuation, repeated characters,
and converted the text to lower-case before applying the dictionary. We used a
part-of-speech tagger trained on the IJS JOS-1M corpus to perform the tagging
Virag (2014). We simplified the results by considering only the part of speech
and its type. We characterized each message by the number of occurrences of
each label which can be viewed as applying a bag-of-words model with ’words’
being the part-of-speech tags.

####  3.1.5 Time-Based Features and Models

The features described so far consider each message as an independent unit and
do not take into consideration the inherent time dependency. By observing the
temporal distribution of labels, we can observe that the labels are not
uniformly distributed. Messages relevant to the book seem to appear in
clusters and observing a message marked as a question naturally leads us to
expect an answer in the subsequent messages. We can use the sequence of labels
in the dataset to compute a label transition probability matrix defining a
Markov model. We can also compute the conditional probabilities of each label
based on the labels of the previous messages. We implemented both models with
conditional probabilities computed given the previous 4 labels. We also
assigned to each message the number of times the poster has posted in a row
and the number of messages authored by the posters in the last 20 messages.

###  3.2 Class Label Distribution and Resampling

It is important to inspect the distribution of class labels in any dataset and
note any severe imbalances that can cause problems in the model construction
phase as there may not be enough data to accurately represent the general
nature of the underrepresented group. Such imbalances also warrant care in
result interpretation.

Figure 1 shows the distribution of class labels for each of the prediction
objectives. We can see that the distribution of broad category labels is
notably imbalanced with 40.3% of messages assigned to the broad category of
’chatting’, but only 1%, 4.5% and 8% to ’switching’, ’moderation’ and ’other’
respectively.

Figure 1: Class label distribution for different prediction objectives.

Several methods have been developed to handle such imbalances, the most
notable being resampling and the use of cost-sensitive classification methods.
Resampling techniques aim to balance the dataset by removing data from the
overrepresented classes, generating synthetic data for the minority classes,
or both. We evaluate the use of SMOTE resampling in tandem with cleaning
unwanted overlap between classes using TOMEK links Batista et al. (2003).

###  3.3 Feature Ranking and Feature Selection

The relative values of features for building a quality predictive model often
vary significantly. Determining the importance of features is a non-trivial
task for which many metrics and methods have been proposed. The subset of
features used to build the model can have an important effect on its
performance and overall usefulness. A model induced on a well-chosen feature
subset will be more general and easier to interpret. A notable group of
algorithms well-suited for the task of ranking features is the family of
Relief-like algorithms which offer performance acceptable for use with large
datasets as well as the notable capability of detecting feature-feature
interactions. These algorithms work by sampling training data instances and
scoring the attributes based on how well they separate the sampled instances
from closest instances corresponding to a different class as well as on the
similarity to closest instances from the same class by this attribute
Kononenko et al. (1997). We use the SWRF* (Sigmoid Weighted ReliefF Star)
algorithm to perform the feature ranking Stokes and Visweswaran (2012). We
also estimate the importance of features by observing the coefficients of a
fitted logistic regression model and analyzing the Gradient boosting model
fitted to the training data.

###  3.4 Prediction Models

The rapidly advancing field of machine learning has produced a myriad of
methods that can be used to make predictions in a supervised learning setting.
We evaluate the use of well-known classification algorithms such as Random
forests, Support vector machines, Gradient boosting, and logistic regression.
We also implement and evaluate an ensemble method called Feature stacking that
was specially developed for such tasks. It is important to critically compare
any results obtained by such sophisticated methods to the outputs of baseline
models such as the so-called majority classifier, which always predicts the
most common label found in training data with maximal certainty as well as the
random-guess or uniform classifier which, predicts the labels uniformly at
random. To be useful, any implemented method should be statistically proven to
outperform these trivial baselines.

####  3.4.1 The Feature Stacking Method

A notable algorithm developed for the task of sentence classification that is
especially suited for high-dimensional data is the so-called Feature stacking
approach Lui (2012) which combines several classifiers each trained on a
feature subset and a final classifier that takes as its input the output of
these classifiers. The authors of the paper describing the Feature stacking
approach suggest using logistic regression as the feature classification
method since it makes the method resemble a neural network structure. We use
an SVM as the final classifier.

During training, the training data is converted to new features consisting of
logistic regression outputs for each feature subset. This is achieved using
k-fold cross-validation. Next, logistic regression is fitted to the entire
training data feature subsets and is used to encode the test data. A final
meta-classifier is fitted to the training data encoded using logistic
regression. Test data is first encoded using a trained logistic regression
model and finally classified with the meta-classifier. The features stacking
method is contrasted with the more typical feature concatenation method on the
diagram shown in figure 2.

Figure 2: The feature stacking method contrasted with the typical method of
concatenating features (taken from Lui (2012)).

####  3.4.2 Utilizing Predictions of Temporal Models

We combined the predictions of the classification model with the probabilities
computed using the Markov model and the conditional probabilities by weighing
the class probabilities obtained by each method as shown below.

$\displaystyle p(C)=$ $\displaystyle(1-\alpha-\beta)\cdot p_{c}(C)$
$\displaystyle+\alpha\cdot p_{m}(C)$ $\displaystyle+\beta\cdot
p(C|\textrm{previous n labels})$

Here, $p_{c}$ and $p_{m}$ represents the probability obtained by the
classification model and the Markov model respectively. We performed an
exhaustive cross-validated grid-search to tune the $\alpha$ and $\beta$
parameters of the combined ensemble model.

##  4 Experiments

###  4.1 Experimental Setup

We began by building classification models using only general features
obtained by observing the word counts, word lengths, and character properties
in individual messages. We compared the use of different models and performed
feature scoring to rank the perceived usefulness of each feature. All model
evaluations were performed using 10 repetitions of 10-fold cross-validation.
Explicit comparisons between different methods were made using the Bayesian
correlated t-test which can be used to compute probabilities of one method
being better than the other and avoids some common pitfalls associated with
the more typical frequentist approaches Benavoli et al. (2016). Subsequently,
the next batch of features was added to the feature extraction/engineering
process and the evaluation process repeated. for the initial feature subset
evaluations, we focused exclusively on the book relevance prediction
objective. Using the full feature set, we evaluate the best scoring models on
all prediction objectives.

###  4.2 Experimental Results

####  4.2.1 The Initial Feature Subset

Table 1 shows the results obtained by evaluating the support vector machine
model built using the starting set of features. Using the Bayesian correlated
t-test, we estimated the support vector machine model to be better than the
random forest model with a probability equal to $0.59$. The probability of the
models performing equally was estimated to be $0.39$. The support vector
machine was also estimated to be marginally better than the gradient boosting
model with a probability of 0.48 and a probability of 0.47 of them being
equal. We did not use the feature stacking method as it is not well defined
for data with few features. The model outperformed the baseline random and
majority models with an estimated probability of 1.0.

| not relevant | relevant
---|---|---
precision | 0.801 | 0.696
recall | 0.807 | 0.691
f1-score | 0.798 | 0.702
support | 215.5 | 138.6
Table 1: Results for the book relevance classification objective obtained by
evaluating the support vector machine model using the initial feature subset.

We used the gradient boosting algorithm, logistic regression coefficients, and
the SWRF* algorithm to estimate the discriminativity of features in the
initial subset. Averaging the estimations, the average word length and the
word count of the message were deemed most important, followed by the maximal
word length and the amount of punctuation in the message.

Figure 3 shows the separability of relevant and non-relevant messages by the
two top-rated features. We can see that a notable portion of non-relevant
messages has a considerably higher average word length. This can be explained
by observing that many non-relevant messages contain long gibberish words as
well as words with consecutively repeated letters.

Figure 3: Seperability of relevant and non-relevant messages by average word
length and word count.

####  4.2.2 Including Important Words and User’s Post History

We augmented the initial feature subset with counts of curse words, repeated
letters, counts of special verbs and nouns deemed important, such as ’misliti’
(to think), ’knjiga’ (book), counts of common Slovene given names, counts of
chat usernames, the number of times the poster posted in a row and the portion
of poster’s posts in the last 20 messages. We repeated the evaluation process
for the book relevance prediction objective. Using the Bayesian correlated
t-test, the SVM classifier was again evaluated to be better by a small margin.

Table 2 shows the results obtained by evaluating the support vector machine
model build using the augmented set of features.

| not relevant | relevant
---|---|---
precision | 0.825 | 0.731
recall | 0.828 | 0.727
f1-score | 0.826 | 0.728
support | 215.5 | 138.6
Table 2: Results for the book relevance classification objective obtained by
evaluating the support vector machine model using the augmented feature
subset.

Using the Bayesian correlated t-test, the support vector machine model built
using the augmented feature set was estimated to be better than the one built
using the initial feature set with a probability of 0.97. The probability of
them being equal was estimated to be 0.02.

####  4.2.3 Including Bag-of-Words Features

We proceeded by adding the unigram and bigram counts to the feature set. This
increased the dimensionality of the dataset substantially by adding 2009
additional features. Using the Bayesian correlated t-test, the feature
stacking method was determined as the most probable best classification model.

Table 3 shows the results obtained by evaluating the feature stacking method
model build using the enriched set of features.

| not relevant | relevant
---|---|---
precision | 0.872 | 0.851
recall | 0.918 | 0.706
f1-score | 0.894 | 0.772
support | 215.5 | 138.6
Table 3: Results for the book relevance classification objective obtained by
evaluating the feature stacking method model using additional bag-of-words
features.

The estimated probability that the new feature set produces a better model is
equal to 1.0.

####  4.2.4 Including Part-Of-Speech Tagging Features

Next, we included the Part-of-Speech tagging based features consisting of the
part of speech and its type pair counts. This added 30 additional features.
Again, the most probable best model was evaluated to be the feature stacking
method.

The comparison between feature stacking method models either using POS
tagging-based features or not indicates that the new features do not improve
the model for this prediction objective. The probability of the model built
using the previous feature set being better is estimated to be 0.25 while the
probability of the models being equal is estimated to be 0.67.

####  4.2.5 Including the Temporal Models

We used an exhaustive, cross-validated grid search to determine the optimal
weights of the Markov chain model and the conditional probability-based model,
which were estimated to be 0.06 and 0.07 respectively. We compared the feature
stacking model with and without the combined use of temporal models using the
Bayesian correlated t-test for the book relevance prediction objective. The
ensemble model combining the predictions of temporal models was evaluated to
be better with a probability of 0.16 while the probability of the models being
equal was estimated to be 0.74. The posterior distribution used to compute the
probabilities is shown in figure 4.

Figure 4: Posterior distribution used to compute the probabilities of one
classification model being better than the other. The region of practical
equivalence is denoted by the vertical bars.

####  4.2.6 Evaluating the Models using All Prediction Objectives

Figure 5 shows the confusion matrix for the book relevance prediction
objective using an 80/20 train-test split and the feature stacking method
model. It is useful to inspect the messages corresponding to false negatives
or false positives and compare them to correctly classified messages to try to
determine ways the feature set could be improved. Table 4 lists a subset of
false negatives and false positives. We can see that the actual label can be
extremely dependent on the context of the conversation which makes it very
difficult for a model with limited ability to process such context to
correctly classify messages shown in the table.

Message | Predicted | True
---|---|---
zakaj | no | yes
aja | no | yes
prestrašeno | no | yes
itak | no | yes
MMMMMMMMMMMMMMMMM NE VEM | no | yes
Jakob, kaj pomeni razdal? Nadaljuj z branjem, prosim. | yes | no
AMPAK POTEM SI UMAZAN | yes | no
Odlična ideja, Smartno B19, prosim nadaljuj z branjem. | yes | no
SmartnoB23, kaj pa ti meniš o odgovoru na vprašanje? | yes | no
jaz se tudi strinjam z smartno 11 | yes | no
Table 4: A subset of false negatives and false positives for the book
relevance prediction objective.

Figure 6 shows the ROC curve and the AUROC value obtained using the same
train-test split and model.

Figure 5: Confusion matrix obtained using a 80/20 train-test split. Figure 6:
ROC curve and the AUROC value obtained using a 80/20 train-test split.

We proceeded by evaluating the models using the type and broad category
prediction objectives using the full feature set. We report the results for
the feature stacking method which was estimated by the Bayesian correlated
t-test to have the highest probability of being the best model in the
evaluated set of models.

Table 5 shows the results obtained by evaluating the feature stacking method
for the message type prediction objective.

| A | Q | S
---|---|---|---
precision | 0.708 | 0.773 | 0.784
recall | 0.734 | 0.706 | 0.788
f1-score | 0.718 | 0.737 | 0.785
support | 115.7 | 67.3 | 171.1
Table 5: Results for the type classification objective obtained by evaluating
the feature stacking method model using the complete feature set.

Figure 7 shows the confusion matrix obtained using a train-test split for this
classification objective.

Figure 7: Confusion matrix obtained using a 80/20 train-test split. ’S’ refers
to the ’statement’ label, ’Q’ to the ’question’ label, and ’A’ to the ’answer’
label.

Table 6 shows the results obtained by evaluating the feature stacking method
for the broad category prediction objective.

| C | D | I | M | O | S
---|---|---|---|---|---|---
precision | 0.632 | 0.683 | 0.895 | 0.878 | 0.859 | 0.0
recall | 0.766 | 0.796 | 0.496 | 0.360 | 0.321 | 0.0
f1-score | 0.691 | 0.735 | 0.635 | 0.503 | 0.460 | 0.0
support | 143.0 | 119.7 | 41.7 | 17.7 | 28.2 | 3.8
Table 6: Results for the broad category classification objective obtained by
evaluating the feature stacking method model using the complete feature set.

Figure 8 shows the confusion matrix obtained using a train-test split for this
classification objective. As stated earlier, it should be noted that the label
distribution for this prediction objective is notably imbalanced. Consequently
the proportion of samples labeled ’moderation’, ’other’, and ’switching’ in
the test set was 6.3%, 4.7%, and 3.2% respectively.

Figure 8: Confusion matrix obtained using a 80/20 train-test split. ’C’ refers
to the ’chatting’ label, ’D’ refers to the ’discussion’ label, ’I’ refers to
the ’identity’ label, ’M’ refers to the ’moderation’ label, ’O’ refers to the
’other’ label, and ’S’ refers to the ’switching’ label.

All models outperformed the baselines with the estimated probability of 1.0.

##  5 Conclusion and Future Work

The best results were achieved by using the Feature stacking method model
built on the complete feature subset. The results indicate the performance to
be sufficient for the methods to be used in real-world tools and platforms. A
significant portion of the information needed for correct classifications is
hidden in the strong temporal interdependence of the messages which our
developed methods exploited only marginally.

## References

  * Batista et al. (2003) Gustavo Batista, Ana Bazzan, and Maria-Carolina Monard. 2003.  Balancing training data for automated annotation of keywords: a case study.  pages 10–18.
  * Benavoli et al. (2016) Alessio Benavoli, Giorgio Corani, Janez Demsar, and Marco Zaffalon. 2016.  Time for a change: a tutorial for comparing multiple classifiers through bayesian analysis.
  * Domingos (2012) Pedro Domingos. 2012.  A few useful things to know about machine learning.  _Commun. ACM_ , 55(10):78–87.
  * Fišer (2020) Darja Fišer. 2020.  The janes project.  http://nl.ijs.si/janes/.  Accessed: 2020-05-1.
  * Kadunc (2016) Klemen Kadunc. 2016.  Določanje sentimenta slovenskim spletnim komentarjem s pomočjo strojnega učenja.
  * Kononenko et al. (1997) Igor Kononenko, Edvard Šimec, and Marko Robnik-Šikonja. 1997.  Overcoming the myopia of inductive learning algorithms with relieff.  _Applied Intelligence_ , 7(1):39–55.
  * Lui (2012) Marco Lui. 2012.  Feature stacking for sentence classification in evidence-based medicine.  In _Proceedings of the Australasian Language Technology Association Workshop 2012_ , pages 134–138, Dunedin, New Zealand.
  * Mclaren et al. (2007) Bruce Mclaren, Oliver Scheuer, Maarten Laat, Rakheli Hever, Reuma Groot, and Carolyn Rosé. 2007.  Using machine learning techniques to analyze and support mediation of student e-discussions.  pages 331–338.
  * Stokes and Visweswaran (2012) Matthew E. Stokes and Shyam Visweswaran. 2012.  Application of a spatially-weighted Relief algorithm for ranking genetic predictors of disease.  _BioData mining_ , 5(1):20–20.
  * Virag (2014) Jernej Virag. 2014.  slo_pos.  https://github.com/izacus/SlovenianLemmatizer.  Accessed on 2019/05/21.
  * Yin et al. (2009) Dawei Yin, Zhenzhen Xue, Liangjie Hong, Brian D. Davison, and Lynne Edwards. 2009\.  Detection of harassment on web 2.0.  In _Proceedings of the Content Analysis in the WEB 2_.