Lampropoulos and Kavousanakis

*Michail E. Kavousanakis, School of Chemical Engineering, National Technical University of Athens, Iroon Polytechneiou 9, Zografou 15780, Athens, Greece.

# Assessment of intra-tumor heterogeneity in a two-dimensional vascular tumor
growth model

Ioannis Lampropoulos  Michail Kavousanakis  School of Chemical Engineering,
National Technical University of Athens, Iroon Polytechneiou 9, Zografou 157
80, Athens, Greece<EMAIL_ADDRESS>
(26 April 2016; 6 June 2016; 6 June 2016)

###### Abstract

[Summary] We present a two-dimensional continuum model of tumor growth, which
treats the tissue as a composition of six distinct fluid phases; their
dynamics are governed by the equations of mass and momentum conservation. Our
model divides the cancer cells phase into two sub-phases depending on their
maturity state. The same approach is also applied for the vasculature phase,
which is divided into young sprouts (products of angiogenesis), and fully
formed-mature vessels. The remaining two phases correspond to healthy cells
and extracellular material (ECM). Furthermore, the model foresees the
existence of nutrient chemical species, which are transferred within the
tissue through diffusion or supplied by the vasculature (blood vessels). The
model is numerically solved with the Finite Elements Method and computations
are performed with the commercial software Comsol Multiphysics ®. The
numerical simulations predict that mature cancer cells are well separated from
young cancer cells, which form a protective shield for the growing tumor. We
study the effect of different mitosis and death rates for mature and young
cancer cells on the tumor growth rate, and predict accelerated rates when the
mitosis rate of young cancer cells is higher compared to mature cancer cells.

###### keywords:

tumor heterogeneity, finite elements, cellular sub-populations, multiphase
model

††articletype: Article Type00footnotetext: Abbreviations: ITH, Intra-tumor
heterogeneity; ECM, extracellular material; FEM, finite elements method

##  1 Introduction

Malignant tumors are dynamic, self-organizing, and heterogeneous biological
structures. During their course, tumors become more heterogeneous and composed
of a diverse collection of cells with different phenotypic behaviors, whether
they belong in the same, or different tumors 1. Unfortunately, evidence
suggests that heterogeneity provides fuel for resistance to various therapies,
and it is often the cause for treatment failure and cancer relapse due to the
selective advantage it potentially provides 2, 3.

Tumor heterogeneity is mainly classified into two major categories: (a) Inter-
Tumor Heterogeneity, which refers to the divergent behavior between different
tumors of the same type and (b) Intra-Tumor Heterogeneity, which indicates
differences in cancer cell sub-populations within the same tumor, as a result
of an ongoing evolutionary process accelerated by the piled up genetic and
epigenetic mutations 1, 4, 5. A third heterogeneity classification, called
Inter-Site Heterogeneity, refers to the formation of heterogeneous tumors in
the body of the same patient 6. In this study, the focus is steered towards
Intra-Tumor Heterogeneity (ITH), i.e. the observation of different sub-
populations within the same malignant neoplasm.

ITH is a widely observed phenomenon in cancer research. All kinds of malignant
tumors consist of cells that are inherently unstable, prone to mutations and
highly diverse. ITH is caused by genetic and epigenetic mutations,
corresponding to differences in the phenotype of the cells constituting the
same tumor 7, 8, 9. Indeed, there is evidence suggesting that ITH dynamics are
not entirely attributed to gene mutations 7, 10. Regardless of the causes, ITH
has a radically negative impact on the success rate of cancer treatment, being
a major contributor to a tumor’s resistance to targeted therapy 8. Cells
belonging to different sub-populations antagonize each other and the healthy
cells, in a struggle to secure precious nutrients for their needs. Because of
that, tumors evolve in a fashion that could be considered Darwinian and
cellular heterogeneity is a decisive factor for the survival of such an
unstable ecosystem 11, 12. Heterogeneity increases the probability for
resistance to therapy and is often associated to various factors of the
tumor’s immediate microenvironment. Hinohara and Polyak 13 have studied both
genetic and transcriptomic heterogeneity and discovered that drug tolerance in
cancer cells can be characterized by specific and reversible epigenetic
states. In this context, they have shown that targeting heterogeneity can
increase treatment efficacy, while fighting the rise of new and resilient cell
populations.

Cancer heterogeneity has increasingly attracted the scientific community’s
attention during the last few years, as shown by the continuously increasing
number of relative theoretical and experimental studies. Baslan et al. 14
detected cancer cells of the same tumor carrying different numbers of copies
of the same genetic material. Losic et al. 7 studied ITH in liver cancer,
concluding that ITH has critical impact in cancer treatment, since sampling a
particular sub-population from a patient’s tumor can lead to false assumptions
about the state of the illness and derail treatment 7, 10. McDonald et al. 15
concluded that ITH correlates with a lesser immune response and an overall
smaller survival rate in breast cancer. Gay et al. 8 suggest that cancer
heterogeneity is heavily influencing cancer treatment, as it is inevitable
that tumors will show heterogeneous sub-populations of cells with different
responses to treatments. Unless this is taken into account, the resistant sub-
populations may survive target treatment and repopulate, breeding a new tumor
more resistant compared to its predecessor. It goes without saying that,
quantifying ITH is a big opportunity to improve the success rate of cancer
treatments.

ITH is also linked with the tendency of cells to metastasize. Miller et al. 16
injected cells from the same mouse mammary tumor into different mice, in an
effort to distinguish variance in the injected tumors’ behavior. Indeed, their
research showed heterogeneous stability levels in the metastatic properties of
the inserted cancerous sub-populations. For further insights into the concept
of tumor heterogeneity, the reader is referred to the reviews of Fisher et al
1, Gay et al. 8, Hinohara and Polyak 13, Zeng and Dai 17, Koltai 18, and
Meacham and Morrison 19.

In order to enhance our quantitative understanding of the clinical and
experimental settings, as well as to illuminate the underlying bio-medical
phenomena, several mathematical models have been developed for the simulation
of heterogeneous tumor growth. Zhao et al. 20 developed an optimization
algorithm in an effort to estimate the effect of various drug combinations on
heterogeneous tumors. Their predictions were based on the already measured
performance of the drugs being administered one at a time on various cell sub-
populations. The efficacy of a drugs combination was calculated as a weighted
average of its single-drug counterparts on a homogeneous tumor. Alvarez et al.
21 studied the tumor-immune dynamics of heterogeneous tumors through a non-
linear mathematical model of cancer immuno-surveillance. The model takes into
account ITH by distinguishing two cancer cell sub-populations based on their
immunogenicity, i.e. their ability to trigger an immune response. On the topic
of immunosurveillance one can also refer to the reviews of Dunn et al. 22, 23.
Harris et al. 24 contemplated upon the merits of different approaches of
heterogeneous tumor modeling, namely agent-based models, population dynamics
and multiscale models. Stamatelos et al. 25 presented a 3D hybrid image-based
modeling framework and demonstrated the correlation between the underlying
irregularities in vascular growth and the emergence of ITH.

Heterogeneity is also intertwined with stemness. In a tumor where ITH is
expressed, there have been confirmed cases of cellular sub-populations
expressing a phenotype resembling that of stem cells 26, 27, 28, 29, 30, 31.
Tripathi et al. 32 developed a model targeting epithelial-mesenchymal
heterogeneity, a phenomenon observed amongst the cells produced through
mitotic division of a cancer cell. Their computations succeeded in capturing
the changes observed in the volume fractions of different cancerous phenotypes
coexisting in tumors of the prostate. Furthermore they were able to predict
the tumor’s composition of different phenotypes in cells sub-populations.
Interestingly enough, they concluded that taking simultaneous action on both
epithelial and mesenchymal cells will possibly increase the likelihood for
tumor growth restriction.

In this study, the heterogeneity element we focus on is the distinction of
cancer cells based on their maturity state. As it is going to be further
explained in later paragraphs, the cells that are presented in the model are
divided into two sub-populations: the young and mature cells. The state of
maturity in cancerous populations is an emerging topic in the medical science.
Differentiation therapy is an aspiring method of treating cases of aggressive
tumors by inducing cancer cells maturation and differentiation through the
insertion of the appropriate chemical agents into the tumor’s microenvironment
33, 34. Of course, maturity in cancer cells has not been studied only in
relation with the above technique. Sontag et al. 35 studied the possible
correlation between a cancer cell’s thermal resistance and its maturity state,
posing the initial hypothesis that cancer stem cells found in human
hepatocellular carcinoma are exhibiting higher resilience to high
temperatures, which are produced in hyperthermic techniques often used to
treat hepatocellular carcinoma 35, 36, than that of the non-stem and mature
cells.

It is beyond doubt that vasculature also holds a very important role in the
phenomenon of cancer growth. Its presence is essential in malignant tumors as
vessels generated through angiogenesis are the prelude to metastasis 3738.
Heterogeneity is also encountered in the vasculature generated during
angiogenesis 21, 39, 40. Heterogeneity in the surrounding vasculature of the
tumor can be considered a leading cause for ITH and the irregularities it can
cause may hinder drug delivery in the tumor area 25, 41. Vasculature produced
internally by the cancer cells differ significantly from usual vasculature 37,
42, 43, 44. Pries et al. 37 performed a series of simulations in an effort to
correlate vessels structural properties to oxygen distribution in the tumor
area and response to different types of stimuli. Vessels produced through
angiogenesis are immature and permeable, among other characteristics 37, 45,
46.

In the present study, we attempt to incorporate the different stages of
maturity by dividing the vascular phase into two sub-phases: the young and
mature vessels. As it will be presented in the next paragraphs, these two
phases exhibit different properties in an attempt to simulate the experimental
observations about microvasculature and hypoxic regions within the tumor. The
model presented in this paper is an extension of the model of Hubbard and
Byrne 47, which simulates the growth of a vascular tumor and incorporates
interactions between cancer cells, healthy cells, vasculature, extracellular
material and the supplied nutrient. Here, we incorporate ITH by modeling the
existence of distinct sub-populations of cancer cells, as well as of
heterogeneous vasculature, and quantify their effect on the overall dynamics
of a growing tumor. We also study the spatial tumor heterogeneity, by
inspecting the relative position of different sub-populations within the
tissue, which hosts healthy, cancer cells, vessels and extracellular matrix
material (ECM). The model is numerically solved with the method of Finite
Elements, and computations are performed with Comsol Multiphysics ®.

The present paper is structured as follows: in Section 2, we present the
details of the developed computational model, which incorporates cancer cell,
and vessel heterogeneity during the growth of a vascular tumor. In Section 3,
we present the results of simulations for different scenarios, where cancer
cell heterogeneity is studied under the prism of different mitosis, and death
rates, as well as for different material properties of the modelled phases.
The final Section 4 is devoted to the discussion of the produced results and
the suggestions for future research.

##  2 Model development

In this section, we describe a multiphase model of vascular tumor growth,
which incorporates cancer cell, and vasculature heterogeneity. It is a two-
dimensional, continuous model treating the tissue as a composition of six
distinct and interacting phases. For each of the phases, we model the
evolution of their volume fraction, and in particular the volume fraction of:
(i) healthy cells, which is denoted with $\theta_{h}$, (ii) young cancer
cells, denoted with $\theta_{yc}$, (iii) mature cancer cells, $\theta_{mc}$,
(iv) young vessels, $\theta_{yv}$, (v) mature vessels, $\theta_{mv}$ and (vi)
extracellular material (ECM), $\theta_{ecm}$. The aforementioned phases are
modelled as viscous fluids exhibiting macroscopic flow, and we denote their
velocities and pressures with $\vec{u}_{i}=\left(u_{i},v_{i}\right)$ and
$p_{i}$, respectively for $i=h,yc,mc,yv,mv,ecm$. The model also stipulates the
presence of a diffusible nutrient, $c$, which is supplied by vessels (young
and mature), is transported through diffusion within the tissue and gets
consumed by normal and cancer cells for survival and proliferation purposes.

Volume fractions of all six phases, as well as their velocities and pressures
are derived by applying to each phase the principle of mass and momentum
balance with appropriate constitutive equations; the concentration of the
diffusible nutrient is computed by solving a reaction-diffusion equation.

###  2.1 Mass balance equations

Under the assumption that density in a living tissue is uniform and constant,
the mass balances for the six phases are simplified as follows:

$\frac{\partial\theta_{i}}{\partial
t}+\nabla\cdot\left(\vec{u}_{i}\theta_{i}\right)=q_{i},\mbox{ for
}i=h,yc,mc,yv,mv,ecm.$ (1)

The term, $\nabla\cdot\left(\vec{u}_{i}\theta_{i}\right)$ formulates mass
transfer through convection, $\theta_{i}$ is the volume fraction of phase $i$
and $\vec{u}_{i}$ its velocity. The source term, $q_{i}$, is associated with
processes including cell mitosis (proliferation), death, angiogenesis, vessel
occlusion (in general processes triggering transfer of mass between phases).

####  2.1.1 Healthy cells

The mass balance for the healthy cells is formulated as follows:

$\frac{\partial\theta_{h}}{\partial
t}+\nabla\cdot\left(\vec{u}_{h}\theta_{h}\right)=k_{m,h}\theta_{h}\theta_{ecm}\left(\frac{c}{c_{p}+c}\right)-k_{d,h}\theta_{h}\left(\frac{c_{c_{1}}+c}{c_{c_{2}}+c}\right),$
(2)

where $k_{m,h}$ is the healthy cells mitosis rate constant, $k_{d,h}$ is the
healthy cells death rate constant. The first term on the right hand side of
Equation (2) models mitosis, with $c_{p}$ denoting the nutrient concentration
at which the mitosis rate becomes half maximal. The second term on the right
hand side of Equation (2) represents the death rate of healthy cells;
$c_{c_{1}},c_{c_{2}}$ denote threshold nutrient concentration values
regulating the cellular death rate. Finally, in order to ensure that the death
rate increases as the nutrient concentration decreases, we set
$c_{c_{1}}>c_{c_{2}}$ .

####  2.1.2 Cancer cells

The cancer cells are grouped into two sub-populations depending on their state
of maturity (age). The first age group includes young cancer cells and the
second group includes mature cells, which originate from young cells through a
maturation process. In our model, mitosis and death rates of cancer cells can
be altered depending on their age.

The mass balance equations for young and mature cancer cells are:

$\frac{\partial\theta_{yc}}{\partial
t}+\nabla\cdot\left(\vec{u}_{yc}\theta_{yc}\right)=\left(k_{m,yc}\theta_{yc}+k_{m,mc}\theta_{mc}\right)\theta_{ecm}\left(\frac{c}{c_{p}+c}\right)-k_{d,yc}\theta_{yc}\left(\frac{c_{c_{1}}+c}{c_{c_{2}}+c}\right)-k_{mat}^{cell}\theta_{yc},$
(3)

$\frac{\partial\theta_{mc}}{\partial
t}+\nabla\cdot\left(\vec{u}_{mc}\theta_{mc}\right)=-k_{d,mc}\theta_{mc}\left(\frac{c_{c_{1}}+c}{c_{c_{2}}+c}\right)+k_{mat}^{cell}\theta_{yc},$
(4)

where $k_{m,yc}$ and $k_{d,yc}$ denote the young cancer cells mitosis and
death rate constant, respectively. For mature cancer cells, the mitosis and
death rate constants are denoted with, $k_{m,mc}$ and $k_{d,mc}$,
respectively. The first term on the right hand side of Equation (3) represents
the mitosis rate, which depends on the nutrient concentration, $c$ and the
volume fraction of ECM, $\theta_{ecm}$. The mitosis rate is proportional to
the volume fraction of cancer cells (both young and mature); through the cell
birth process young cancer cells will be produced.

For simplicity purposes, we set the threshold nutrient concentration values,
$c_{p},c_{c_{1}},c_{c_{2}}$ for mitosis and death rates to be identical for
normal and cancer cells. In order to take into account that cancer cells
proliferate more rapidly, and are less prone to death compared to healthy
cells, we set $k_{m,yc},k_{m,mc}>k_{m,h}$ and $k_{d,yc},k_{d,mc}<k_{d,h}$.
Finally, we formulate the maturation process through which young cells become
mature ones, with the last term of the right hand side of Equations (3) and
(4); maturation is formulated as a first order reaction proportional to the
volume fraction of young cells, $\theta_{yc}$, and $k_{mat}^{cell}$ is the
cell maturation constant rate. We define the value of the cell maturation rate
constant, $k_{mat}^{cell}$ through:

$k_{mat}^{cell}=\beta\cdot k_{m,yc},$ (5)

where $\beta$ is a constant, and $0<\beta<<1$, i.e., the maturation rate
constant is significantly smaller compared to the mitosis rate constant of
young cancer cells.

####  2.1.3 Vasculature

An additional source of intrinsic heterogeneity is the differentiation of
blood vessels. We adopt a similar approach with the division of the cancerous
phase into young and mature cells, and separate the blood vessels phase into
two distinct phases: young vessels (sprouts formed due to angiogenesis) and
mature vessels. The mass balance equations for young and mature vessels
respectively are:

$\begin{split}\frac{\partial\theta_{yv}}{\partial
t}+\nabla\cdot\left(\vec{u}_{yv}\theta_{yv}\right)&=k_{ang}\left(\theta_{h}+\theta_{yc}+\theta_{mc}\right)\left(\theta_{yv}+\theta_{mv}\right)\left(\frac{\theta_{ecm}}{\epsilon+\theta_{ecm}}\right)\frac{c}{\left(c_{a}+c\right)^{2}}\\\
&-k_{occ}^{yv}\theta_{yv}\mathcal{H}\left(\theta_{h}p_{h}+\theta_{yc}p_{yc}+\theta_{mc}p_{mc}-p_{crit}^{yv},h\right)-k_{mat}^{ves}\theta_{yv},\end{split}$
(6)

$\frac{\partial\theta_{mv}}{\partial
t}+\nabla\cdot\left(\vec{u}_{mv}\theta_{mv}\right)=-k_{occ}^{mv}\theta_{mv}\mathcal{H}\left(\theta_{h}p_{h}+\theta_{yc}p_{yc}+\theta_{mc}p_{mc}-p_{crit}^{mv},h\right)+k_{mat}^{ves}\theta_{yv},$
(7)

The first term of the right hand side in Equation (6) corresponds to the
angiogenesis process, during which young vessels are produced. We denote with,
$k_{ang}$, the angiogenesis rate constant and the angiogenesis rate is
modelled to be proportional to the total volume fraction of healthy, young and
mature cancer cells, as well as proportional to the total volume fraction of
young and mature vessels. The vessel growth is also assumed to get promoted by
the presence of ECM and the diffusible nutrient, and reaches a maximum rate
for high values of ECM ($\theta_{ecm}$) and nutrient concentration values;
$\epsilon$ corresponds to the ECM volume fraction at which the angiogenesis
rate becomes half maximal, and $c_{a}$ corresponds to the nutrient
concentration at which the angiogenesis rate becomes maximal. We note that the
angiogenesis term appears only as a source term for young vessels.

Vessel occlusion is present in both young (second term of the right hand side
in Equation (6)) and mature vessels (first term of the right hand side in
Equation (7). For both phases the occlusion rate depends proportionally to the
volume fraction of vessels (young and mature, respectively), and gets
triggered when the pressure exerted by the surrounding phases
($\theta_{h}p_{h}+\theta_{yc}p_{yc}+\theta_{mc}p_{mc}$) exceeds a threshold
value, denoted with $p_{crit}^{yv}$ for young and $p_{crit}^{mv}$ for mature
vessels. For this purpose, we model $\mathcal{H}$ to be a smooth transition
function:

$\mathcal{H}\left(\theta_{h}p_{h}+\theta_{yc}p_{yc}+\theta_{mc}p_{mc}-p_{crit},h\right)=\frac{1}{2}\left[1+\tanh\left(\frac{\theta_{h}p_{h}+\theta_{yc}p_{yc}+\theta_{mc}p_{mc}-p_{crit}}{h_{sm}}\right)\right],$
(8)

where $h_{sm}$ is a smoothness parameter, and $p_{crit}=p_{crit}^{yv}$ or
$p_{crit}^{mv}$, for young and mature vessels, respectively.

The major differences in the functionality of vessels belonging to different
age groups refer to their ability to withstand vascular occlusion, with mature
vessels being more resilient compared to young sprouts. Thus, we set different
pressure threshold values, namely: $p_{crit}^{mv}$ for mature vessels, and
$p_{crit}^{yv}$ for young vessels with $p_{crit}^{yv}<p_{crit}^{mv}$. In
addition, we consider that the rate constant of occlusion for mature vessels,
$k_{occ}^{mv}$ is smaller compared to the occlusion rate constant of young
vessels, $k_{occ}^{yv}$, i.e.: $k_{occ}^{mv}<k_{occ}^{yv}$ Finally, we model
the maturation process of blood vessels as a first order reaction, which is
proportional to the volume fraction of young vessels, and $k_{mat}^{ves}$
being the maturation rate constant.

####  2.1.4 Extracellular Material

ECM is a passive medium hosting material that is used for proliferation
purposes, or produced from cellular death and vascular occlusion (necrotic
material). By assuming that the studied system is closed (no external
replenishment of resources) the sum of source terms, $q_{i}$ for all phases
must sum up to $0$:

$\sum_{i}q_{i}=0,\mbox{ for }i=h,yc,mc,yv,mv,ecm.$ (9)

In accordance to the assumption above, the source term for ECM is defined as:

$\frac{\partial\theta_{ecm}}{\partial
t}+\nabla\cdot\left(\vec{u}_{ecm}\theta_{ecm}\right)=q_{ecm}=-\sum_{i}q_{i},\textrm{
for }i=h,yc,mc,yv,mv.$ (10)

In practice, we adopt the assumption that the tissue has no voids:

$\sum_{i}\theta_{i}=1,i=h,yc,mc,yv,mv,ecm,$ (11)

and calculate the volume fraction of ECM, through Equation (11), instead of
computing $\theta_{ecm}$ from Equation (10).

####  2.1.5 Diffusible nutrient

We assume that the timescales of transport and reaction processes for nutrient
are considerably shorter (minutes) compared to the timescales of processes
associated with fluid phase changes (days or weeks); under this assumption,
the nutrient is in a quasi-steady state and the general form of its mass
balance equation can be written as follows:

$-D_{c}\nabla^{2}c=q_{c}.$ (12)

$D_{c}$ is the diffusion coefficient for the nutrient and $q_{c}$ its source
term equal to:

$q_{c}=k_{rep}(\theta_{yv}+\theta_{mv})(c_{v}-c)-(k_{c,h}\theta_{h}c+k_{c,yc}\theta_{yc}c+k_{c,mc}\theta_{mc}c)-\left[k_{cm,h}\theta_{h}+k_{cm,yc}\theta_{yc}+k_{cm,mc}\theta_{mc}\right]\theta_{ecm}\left(\frac{c}{c_{p}+c}\right);$
(13)

$k_{rep}$ is the nutrient replenishment rate constant for the replenishment
performed by the local vasculature, $c_{v}$ is the nutrient concentration
inside said vasculature $k_{c,i},i=h,yc,mc$, denote the consumption rate
constant for cell sustenance of healthy and cancer cells, respectively and
$k_{cm,i},i=h,yc,mc$, the rate constants for the consumption of nutrient for
the purposes of mitosis for healthy and cancer cells, respectively. The value
of $k_{cm,i}$ is given by the formula:

$k_{cm,i}=\frac{k_{cm,h}\cdot k_{m,i}}{k_{m,h}},i=h,yc,mc.$ (14)

Again, we stress that the nutrient contribution to the overall volume of the
tissue is neglected and that its mass transport is induced by diffusion.

###  2.2 Momentum balance equations

Assuming creeping flow, we write the general form of momentum balance equation
for each phase,$i$:

$\nabla\cdot\left(\theta_{i}\boldsymbol{\sigma}_{i}\right)+\vec{F}_{i}=\vec{0},\mbox{
for }i=h,yc,mc,yv,mv,ecm,$ (15)

where $\boldsymbol{\sigma}_{i}$ denotes the stress tensor of phase, $i$, and
$\vec{F}_{i}$ denotes the momentum source terms, which include the effects of
pressure and inter-phase drag:

$\vec{F}_{i}=p_{i}\boldsymbol{I}\nabla\theta_{i}+\sum_{j,j\neq
i}d_{i,j}\theta_{i}\theta_{j}\left(\vec{u}_{j}-\vec{u}_{i}\right),\mbox{ for
}i=h,yc,mc,yv,mv,ecm.$ (16)

Considering each fluid phase as a viscous and compressible fluid, the stress
tensor is:

$\boldsymbol{\sigma}_{i}=-p_{i}\boldsymbol{I}+\mu_{i}\left(\nabla\vec{u}_{i}+\left(\nabla\vec{u}_{i}\right)^{T}\right)-\frac{2}{3}\mu_{i}\left(\nabla\cdot\vec{u}_{i}\right)\boldsymbol{I}.$
(17)

In Equations (16)-(17), $p_{i}$ is the phase pressure, $d_{i,j}$ is the drag
coefficient due to the relative movement of the phases, $i$ and $j$, surfaces;
$\mu_{i}$ is the dynamic viscosity of phase $i$. For the calculation of both
pressures and velocity fields for each phase, a continuity equation for the
phase mixture is required, and can be produced by the summation of the mass
balance equations of all phases (Equation (1) for $i=h,yc,mc,yv,mv,ecm$):

$\sum_{i}\nabla\cdot\left(\vec{u}_{i}\theta_{i}\right)=0,i=h,yc,mc,yv,mv,ecm.$
(18)

Furthermore, appropriately defined constitutive relations correlating the
different phase pressures are also required for the closure of the system of
equations. In particular, we set the vascular phase pressures,
$p_{yv},p_{mv}$, to be equal to a reference pressure, $p_{ref}=0$, (externally
imposed pressure in the vasculature, which is set to 0), and the rest of the
pressures are related through:

$p_{h}=p_{yc}=p_{mc}=p_{ecm}+\Sigma\left(\theta_{h}+\theta_{yc}+\theta_{mc}\right).$
(19)

$\Sigma(\theta)$ is a function calculating the increase in pressure exerted by
cells whenever their local density exceeds their natural value $\theta^{*}$
(the cellular volume fraction value in a healthy tissue):

$\Sigma(\theta)=\begin{cases}\frac{\Lambda(\theta-\theta^{*})}{\left(1-\theta\right)^{2}},&\mbox{if
}\theta\geq\theta^{*}\\\ 0,&\mbox{otherwise}.\end{cases}$ (20)

where $\Lambda$ is a tension constant measuring the tendency of cells to
restore their natural density.

###  2.3 Non-dimensionalised model

Before proceeding to the presentation of the results of the model, we convert
the variables to a dimensionless form. Volume fractions $\theta_{i}$ are
already dimensionless so the rest of the variables are non-dimensionalised as
follows:

$t^{\prime}=k_{m,h}t,\vec{x^{\prime}}=\frac{\vec{x}}{R_{o}},\vec{u^{\prime}}_{i}=\frac{\vec{u}_{i}}{k_{m,h}R_{o}},p^{\prime}_{i}=\frac{p_{i}}{\Lambda},c^{\prime}=\frac{c}{c_{v}},$
(21)

where $R_{o}$ is a typical length scale; here $R_{o}$ is the radius of the
initial cancerous seed. The resulting dimensionless mass balances read:

$\frac{\partial\theta_{h}}{\partial
t}+\nabla\cdot\left(\theta_{h}\vec{u}_{h}^{\prime}\right)=\theta_{h}\theta_{ecm}\left(\frac{c^{\prime}}{c_{p}^{*}+c^{\prime}}\right)-k_{d,h}^{*}\theta_{h}\left(\frac{c_{c_{1}}^{*}+c^{\prime}}{c_{c_{2}}^{*}+c^{\prime}}\right),$
(22)

$\frac{\partial\theta_{yc}}{\partial
t}+\nabla\cdot\left(\theta_{yc}\vec{u}_{yc}^{\prime}\right)=(k_{m,yc}^{*}\theta_{yc}+k_{m,mc}^{*}\theta_{mc})\theta_{ecm}\left(\frac{c^{\prime}}{c_{p}^{*}+c^{\prime}}\right)-k_{d,yc}^{*}\theta_{yc}\left(\frac{c_{c_{1}}^{*}+c^{\prime}}{c_{c_{2}}^{*}+c^{\prime}}\right)-k_{mat}^{cell^{*}}\theta_{yc},$
(23)

$\frac{\partial\theta_{mc}}{\partial
t}+\nabla\cdot\left(\theta_{mc}\vec{u}_{mc}^{\prime}\right)=-k_{d,mc}^{*}\theta_{mc}\left(\frac{c_{c_{1}}^{*}+c^{\prime}}{c_{c_{2}}^{*}+c^{\prime}}\right)+k_{mat}^{cell^{*}}\theta_{yc},$
(24)

$\begin{split}\frac{\partial\theta_{yv}}{\partial
t}+\nabla\cdot\left(\theta_{yv}\vec{u}_{yv}^{\prime}\right)&=k_{ang}^{*}(\theta_{h}+\theta_{yc}+\theta_{mc})(\theta_{yv}+\theta_{mv})\left(\frac{\theta_{ecm}}{\epsilon+\theta_{ecm}}\right)\frac{c^{\prime}}{(c_{a}^{*}+c^{\prime})^{2}}\\\
&-k_{occ}^{yv^{*}}\theta_{yv}\mathcal{H}(\theta_{h}p_{h}^{\prime}+\theta_{yc}p_{yc}^{\prime}+\theta_{mc}p_{mc}^{\prime}-p_{crit}^{yv^{*}},h_{sm}^{*})-k_{mat}^{ves^{*}}\theta_{yv},\end{split}$
(25)

$\frac{\partial\theta_{mv}}{\partial
t}+\nabla\cdot\left(\theta_{mv}\vec{u}_{mv}^{\prime}\right)=-k_{occ}^{mv^{*}}\theta_{mv}\mathcal{H}(\theta_{h}p_{h}^{\prime}+\theta_{yc}p_{yc}^{\prime}+\theta_{mc}p_{mc}^{\prime}-p_{crit}^{mv^{*}},h_{sm}^{*})+k_{mat}^{ves^{*}}\theta_{yv},$
(26)

$\displaystyle\theta_{ecm}=1-\theta_{h}-\theta_{yc}-\theta_{mc}-\theta_{yv}-\theta_{mv},$
(27)

where:

$k_{m,yc}^{*}=\frac{k_{m,yc}}{k_{m,h}},k_{m,mc}^{*}=\frac{k_{m,mc}}{k_{m,h}},$
$k_{d,h}^{*}=\frac{k_{d,h}}{k_{m,h}},k_{d,yc}^{*}=\frac{k_{d,yc}}{k_{m,h}},k_{d,mc}^{*}=\frac{k_{d,mc}}{k_{m,h}},$
$k_{occ}^{mv^{*}}=\frac{k_{occ}^{mv}}{k_{m,h}},k_{occ}^{yv^{*}}=\frac{k_{occ}^{yv}}{k_{m,h}},k_{ang}^{*}=\frac{k_{ang}}{k_{m,h}c_{v}},$
$k_{mat}^{cell^{*}}=\frac{k_{mat}^{cell}}{k_{m,h}},k_{mat}^{ves^{*}}=\frac{k_{mat}^{ves}}{k_{m,h}},$
$c_{p}^{*}=\frac{c_{p}}{c_{v}},c_{c_{1}}^{*}=\frac{c_{c_{1}}}{c_{v}},c_{c_{2}}^{*}=\frac{c_{c_{2}}}{c_{v}},c_{a}^{*}=\frac{c_{a}}{c_{v}},$
$h_{sm}^{*}=\frac{h_{sm}}{\Lambda},p_{crit}^{yv^{*}}=\frac{p_{crit}^{yv}}{\Lambda},p_{crit}^{mv^{*}}=\frac{p_{crit}^{mv}}{\Lambda}.$

Since the mass balance equations are hyperbolic, boundary conditions are only
defined on segments of the domain’s, $\Omega$ boundaries (denoted with
$\partial\Omega$) where: $\vec{u}_{i}\cdot\vec{n}<0$ (inflow segments), with
$\vec{n}$ denoting the outward-pointing unit normal vector of the boundary.
For boundaries with $\vec{u}_{i}\cdot\vec{n}<0$, the imposed boundary
conditions are:

$\theta_{i}=\theta_{i}^{\infty},i=h,yc,mc,yv,mv.$ (28)

The dimensionless mass balance equation for the nutrient is:

$\begin{split}-D_{c}^{*}\nabla^{2}c^{\prime}&=(\theta_{mv}+\theta_{yv})(1-c^{\prime})-(k_{c,h}^{*}\theta_{h}c^{\prime}+k_{c,yc}^{*}\theta_{yc}c^{\prime}+k_{c,mc}^{*}\theta_{mc}c^{\prime})\\\
&-\left[k_{cm,h}^{*}\theta_{h}+k_{cm,yc}^{*}\theta_{yc}+k_{cm,mc}^{*}\theta_{mc}\right]\theta_{ecm}\left(\frac{c^{\prime}}{c_{p}^{*}+c^{\prime}}\right),\end{split}$
(29)

where:
$D_{c}^{*}=\frac{D_{c}}{k_{rep}R_{o}^{2}},k_{c,h}^{*}=\frac{k_{c,h}}{k_{rep}},k_{c,yc}^{*}=\frac{k_{c,yc}}{k_{rep}},k_{c,mc}^{*}=\frac{k_{c,mc}}{k_{rep}},k_{cm,h}^{*}=\frac{k_{cm,h}}{k_{rep}c_{v}},k_{cm,yc}^{*}=\frac{k_{cm,yc}}{k_{rep}c_{v}},k_{cm,mc}^{*}=\frac{k_{cm,mc}}{k_{rep}c_{v}}$.
The mass balance of nutrient is an elliptic partial differential equation,
thus there is no need to identify inflow segments of $\partial\Omega$. In this
work, we impose Neumann type boundary condition for the nutrient
concentration:

$\nabla c^{\prime}\cdot\vec{n}=0.$ (30)

The dimensionless form of momentum balances read:

$\sum_{j,j\neq
i}d_{i,j}^{*}\theta_{i}\theta_{j}\left(\vec{u^{\prime}}_{j}-\vec{u^{\prime}}_{i}\right)-\theta_{i}\nabla\cdot\left(\Lambda^{*}p^{\prime}_{i}\boldsymbol{I}\right)+\nabla\cdot\left[\theta_{i}\left[\mu_{i}^{*}\left(\nabla\vec{u^{\prime}}_{i}+\left(\nabla\vec{u^{\prime}}_{i}\right)^{T}\right)-\frac{2}{3}\mu_{i}^{*}\left(\nabla\cdot\vec{u^{\prime}}_{i}\right)\boldsymbol{I}\right]\right]=\vec{0},$
(31)

for $i,j=h,yc,mc,yv,mv,ecm$. Here
$d_{i,j}^{*}=\frac{d_{i,j}}{d_{h,yc}},\Lambda^{*}=\frac{\Lambda}{d_{h,yc}k_{m,h}R_{o}^{2}},\mu_{i}^{*}=\frac{\mu_{i}}{d_{h,yc}R_{o}^{2}}$.
The boundary conditions for Equations (31) are:

$\boldsymbol{\sigma}_{i}\cdot\vec{n}=\vec{0},\mbox{ for }i=h,yc,mc,yv,mv$ (32)

Since the specification of the normal stress along the whole boundary,
$\partial\Omega$, is not possible for all phases,we specify the velocity for
one phase (here the ECM) in order to obtain a unique solution:

$\vec{u}_{ecm}=\vec{0}.$ (33)

Finally, the dimensionless form of the continuity equation (Equation (18))
reads:

$\sum_{i}\nabla\cdot\left(\theta_{i}\vec{u^{\prime}}_{i}\right)=0,\mbox{ for
}i=h,yc,mc,yv,mv,ecm,$ (34)

and the dimensionless equations of constitutive equations for pressures
(Equation (19) are:

$p^{\prime}_{h}=p^{\prime}_{yc}=p^{\prime}_{mc}=p^{\prime}_{ecm}+\Sigma^{\prime}(\theta_{h}+\theta_{yc}+\theta_{mc}),$
(35)

with $\Sigma^{\prime}(\theta)$:

$\Sigma^{\prime}(\theta)=\begin{cases}\frac{(\theta-\theta^{*})}{\left(1-\theta\right)^{2}},&\mbox{if
}\theta\geq\theta^{*}\\\ 0,&\mbox{otherwise}.\end{cases}$ (36)

##  3 Results

All computations in the present study are performed in the environment of
Comsol Multiphysics®, which is based on the Finite Elements Method (FEM). The
tissue is modelled as circular domain with (dimensionless) radius,
$R_{tissue}=16$. The computational mesh is an unstructured mesh generated with
the method of Delaunay Triangulation, and contains a total of approximately
$20,000$ elements, which results in approximately $570,000$ degrees of
freedom. The required computational time on an Intel®CoreTM i5-3360M CPU @
2.80 GHz is approximately 11 hrs on average for each run presented in this
section.

For simplification purposes, we drop the "$\prime$" and "*" symbols for all
dimensionless variables. All simulations are initialized by seeding a number
of cancer cells on a healthy tissue, which can be considered being at an
equilibrium state 47 and satisfies the following assumptions:

  * •

the fluid is a perfect mixture; as a result, the volume fraction of each phase
is uniformly distributed and shows no macroscopic velocity
($\vec{u}_{i}=\vec{0},i=h,yv,mv,ecm$);

  * •

the tissue shows no sign of contamination from cancer cells (young or mature)
($\theta_{yc}=\theta_{mc}=0$);

  * •

healthy cells maintain a constant, natural cell density: in that state, no
stress is exerted between cells, as this volume fraction value serves as the
limit beyond which cell to cell interactions begin to occur
($\theta_{h}=\theta^{*}=0.6$);

  * •

the pressure of each fluid phase is equal to zero (as a direct consequence of
the above assumption)

Following the assumptions above, the system of Equations (22)-(27) is reduced
to the following set of nonlinear equations:

$\displaystyle\theta^{*}\theta_{ecm}\frac{c}{c_{p}+c}-k_{d,h}\theta^{*}\frac{c_{c_{1}}+c}{c_{c_{2}}+c}$
$\displaystyle=$ $\displaystyle 0,$ (37) $\displaystyle-
k_{occ}^{mv}\theta_{mc}\mathcal{H}\left(-p_{crit}^{mv},h_{sm}\right)+k_{mat}^{ves}\theta_{yv}$
$\displaystyle=$ $\displaystyle 0,$ (38)
$\displaystyle\theta^{*}+\theta_{mv}+\theta_{yv}+\theta_{ecm}$
$\displaystyle=$ $\displaystyle 1,$ (39) $\displaystyle
k_{ang}\theta^{*}(\theta_{mv}+\theta_{yv})\frac{\theta_{ecm}}{\epsilon+\theta_{ecm}}\frac{c}{(c_{a}+c)^{2}}-k_{occ}^{yv}\theta_{yv}\mathcal{H}\left(-p_{crit}^{yv},h_{sm}\right)-k_{mat}^{ves}\theta_{yv}$
$\displaystyle=$ $\displaystyle 0,$ (40)
$\displaystyle(\theta_{mv}+\theta_{yv})(1-c)-k_{c,h}\theta^{*}c-k_{cm,h}\theta^{*}\theta_{ecm}\frac{c}{c_{p}+c}$
$\displaystyle=$ $\displaystyle 0.$ (41)

Solving the model for the volume fractions of
$\theta_{mv},\theta_{yv},\theta_{ecm}$ and the non-dimensionalized
concentration of $c$ yields the initial conditions of the studied system. In
order to obtain a unique solution (four unknown variables for a system of five
equations), we also consider $k_{mat}^{ves}$ as an extra unknown variable.
Naturally, the results produced have to satisfy
$\theta_{mv},\theta_{yv},\theta_{ecm},c\in(0,1)$ and $k_{mat}^{ves}>0$ in
order to be considered as valid.

Using the parameter values presented in Table 1, the obtained solution is the
following:

$\theta_{mv}(x,y,0)=0.015335,\theta_{ecm}(x,y,0)=0.38250,\theta_{yv}(x,y,0)=0.0021625,c(x,y,0)=0.25320,k_{mat}^{ves}=0.029041.$
(42)

###  3.1 Growing tumor with uniformly behaving sub-populations

In this paragraph, we simulate the case of a growing tumor in which both young
and mature cancer cells behave in an homogeneous manner, i.e. mature and young
cancer cells share the same proliferation and death rate constants. The
initial population of cancer cells are seeded in a circular domain of radius,
$R_{o}=1$ around the origin of the computational domain. In particular, the
initial cancerous seed is implanted according to the formula:

$\theta_{yc}(x,y,0)=\begin{cases}0.05\cos^{2}\left(\frac{\pi
r}{2}\right),&\mbox{if }\sqrt{x^{2}+y^{2}}\leq R_{o}=1\\\
0,&\mbox{otherwise},\end{cases}$ (43)

It is assumed that the initial cancerous seed consists exclusively of young
cancer cells. Consequently, the initial distribution of mature cancer cells
is:

$\theta_{mc}(x,y,0)=0.$ (44)

The young cancer cells replace healthy cells, thus the initial condition for
$\theta_{h}$ reads:

$\theta_{h}(x,y,0)=0.6-\theta_{yc}(x,y,0).$ (45)

The initial condition of $\theta_{yv},\theta_{mv},\theta_{emc}$ and $c$ is
provided from (42). In this numerical experiment, the constants describing the
two sub-populations of cancer cells are kept equal, so as to simulate a tumor
with uniform growth and death rates. The values of parameters are presented in
Table 1. For the rest of simulations presented in this paper, we use the same
set of parameter values unless explicitly stated otherwise.

Table 1: Parameter values used in the simulations. Parameter name | Value or expression | Description
---|---|---
$k_{m,yc}$ | $2.0$ | Young cancer cell mitosis rate constant
$k_{m,mc}$ | $2.0$ | Mature cancer cell mitosis rate constant
$k_{d,h}$ | $0.15$ | Healthy cell death rate constant
$k_{d,yc}$ | $0.075$ | Young cancer cell death rate constant
$k_{d,mc}$ | $0.075$ | Mature cancer cell death rate constant
$k_{mat}^{cell}$ | $0.05\cdot k_{m,yc}$ | Cancer cell maturation rate constant
$c_{p}$ | $0.25$ | Threshold nutrient concentration for mitosis
$c_{c_{1}}$ | $0.2$ | Nutrient concentration regulating cellular death rate
$c_{c_{2}}$ | $0.1$ | Threshold nutrient concentration for cellular death rate
$k_{occ}^{mv}$ | $0.095$ | Mature vessel occlusion rate constant
$k_{occ}^{yv}$ | $0.135$ | Young vessel occlusion rate constant
$k_{ang}$ | $0.0029449$ | Angiogenesis rate constant
$c_{a}$ | $0.05$ | Threshold nutrient concentration for angiogenesis rate
$p_{crit}^{mv}$ | $0.31$ | Critical pressure for mature vessel occlusion
$p_{crit}^{yv}$ | $0.26$ | Critical pressure for young vessel occlusion
$\epsilon$ | $0.01$ | Threshold ECM volume fraction for angiogenesis rate
$h_{sm}$ | $0.2$ | Smoothness parameter for vessel occlusion function
$D_{c}$ | $1.0$ | Nutrient diffusion coefficient
$k_{c,h}$ | $0.01$ | Consumption rate constant for sustenance of healthy cells
$k_{c,yc}$ | $0.01$ | Consumption rate constant for sustenance of young cancer cells
$k_{c,mc}$ | $0.01$ | Consumption rate constant for sustenance of mature cancer cells
$k_{cm,h}$ | $0.1$ | Consumption rate constant of nutrient used for mitosis of healthy cells
$k_{cm,yc}$ | $k_{cm,h}\cdot k_{m,yc}$ | Consumption rate constant of nutrient used for mitosis of young cancer cells
$k_{cm,mc}$ | $k_{cm,h}\cdot k_{m,mc}$ | Consumption rate constant of nutrient used for mitosis of mature cancer cells
$\Lambda$ | $0.1$ | Tension constant
$\mu_{i}$ | $10.0$ | Dynamic viscosity
$d_{i,j}$ | $1.0$ | Drag coefficient
$R_{tissue}$ | $16$ | Radius of modelled tissue (circular domain, $\Omega$)

We impose Dirichlet boundary conditions for Equations (22)-26:

$\theta_{i}=\theta_{i}^{\infty},$ (46)

where $\theta_{i}^{\infty}$ is the initial volume fraction of phase $i$, i.e.:
$\theta_{i}^{\infty}=\theta_{i}(t=0)$.

The boundary condition for nutrient, $c$, is:

$\nabla c\cdot\vec{n}=0.$ (47)

Finally, the boundary conditions imposed for the momentum balance Equations
(31) are:

$\boldsymbol{\sigma}_{i}\cdot\vec{n}=\vec{0},i=h,yc,mc,yv,mv$ (48)

and Dirichlet (no-slip) condition for ECM:

$\vec{u}_{ecm}=\vec{0}.$ (49)

In this simulation, both cancer cell sub-populations feature double mitosis
rate constant, and half death rate constant compared to healthy cells. As
shown in Figure 1, the population of young cancer cells rapidly starts to
expand outwards, while in the same time the maturation process converts part
of the young sub-population into mature cancer cells. The increased
proliferation rate of cancer cells (compared to healthy cells) leads to local
increases of the cancerous phase density, overcoming natural levels and
enabling interactions between cells; these interactions are the driving force
for tumor growth.

Figure 1: Spatial volume fraction distribution of young tumor cells (left
panel), and mature tumor cells (right panel) at $t=50,100$ and $150$.  Figure
2: (a) Spatial volume fraction of ECM and (b) spatial concentration
distribution of nutrient, $c$ at $t=50,100$ and $150$.

Furthermore, the relentless proliferation of cancer cells starves the inner
region of the tumor from both nutrient and ECM as illustrated in Figure 2.
However, one can observe the existence of large volume fraction values of ECM
at $t=150$. This is attributed to the formation of a necrotic zone in the
interior of the tumor. In an environment lacking the nutrients essential to
maintain live cells, the cells in the interior of the tumor succumb to
necrosis forming the tumor’s necrotic zone.

The increased cancer cell volume fractions lead to increases in the exerted
pressure around the tumor. As a result, the vasculature in the vicinity of the
tumor breaks down, thus further aggravating the problem of nutrient transport
to the tumor’s interior. As reported previously, mature vessels feature the
ability to withstand increased pressure values compared to young sprouts. This
increase in the resilience of mature vessels is illustrated in Figure 3. In
particular, we depict snapshots of the average ratio of mature vessels over
the total volume fraction of vessels along the radial direction (left vertical
axis), complemented by the average radial distribution of cancer cells volume
fractions (right vertical axis). One can observe that at the vicinity of the
tumor’s front the vessels are predominately mature, as a result of their
increased resilience. Furthermore, throughout the exterior region of the
tumor, the balance shifts in favor of the younger vessels, due to the
angiogenesis process.

Figure 3: Average percentage of mature vessels (grey lines) at $t=0$, $t=50$
(circled grey line), $t=100$ (crossed grey line) and $t=150$ (squared grey
line) along the radial direction of the tissue domain. Black curves depict the
average radial distribution of total cancer cells volume fraction at
$t=50,100,150$, (circled, crossed and squared black line, respectively).

Figure 4 depicts the average radial distribution of young and mature cancer
cells at different time instances. Interestingly enough, younger cells form an
outer layer acting as a protective shield for mature cells. Figure 4 further
establishes that the different age group cells are unevenly distributed along
the radial direction of the developed tumor. Mature cells mostly inhabit the
quiescent and necrotic zones of the tumor, while the proliferating zone’s
cancer cells are predominantly young.

Figure 4: Snapshots of the mean radial distribution of young cancer cells
(black lines) at $t=50$ (circled black line), $t=100$ (crossed black line) and
$t=150$ (squared black line) along the radial direction of the circular
domain, $\Omega$. Grey lines depict the mean radial distribution of mature
cancer cells volume fraction at $t=50,100,150$, (circled, crossed and squared
grey line, respectively).

Due to the uncontrolled proliferation rate of cancer cells, one can observe
the formation of a macroscopic velocity field. In general, they initially
present an outward-pointing motion in an effort to relieve the built up
pressure. After the formation of the quiescent and the necrotic core, cancer
cells start to perform a second motion. Due to the volume fraction gradient
and the consequent pressure gradient developed between the proliferation zone
and the rest of the tumor zones, cancer cells in the interior of the tumor
rush inwards to smoothen that gradient and inevitably succumb to necrosis. The
flux $\theta_{i}\cdot\vec{u}_{i}$ of young and mature cancer cells ($i=yc,mc$,
respectively) is illustrated in Figure 5 at different time instances. During
the initial stages, both phases exhibit an outward-pointing motion; young
cancer cells preserve this behavior even at later stages, whereas mature
cancer cells gradually present the tendency of moving towards the tumor’s
inner core as time progresses. It is thus evident that the absence of
uniformity in spatial distribution of young and mature cancer cells originates
from the different motility of the two sub-populations.

Figure 5: Evolution of cancer cell phases
fluxes,$\theta_{i}\cdot\vec{u}_{i},i=yc,mc$ for $t=50,100$ and $150$. The
vector length is proportional to its magnitude.

###  3.2 Effect of mitosis rate heterogeneity

A measure to assess the growth rate of a tumor is by computing the average
distance of the tumor’s front from the origin of the domain, $R_{mean}$. In
particular, we determine $R_{mean}$ as the average distance of contour,
$\theta_{yc}+\theta_{mc}=\theta^{*}=0.6$ from the origin of the domain. In
Figure 6, we depict the evolution of $R_{mean}$ for different mitosis rate
constants of young and mature cancer cells. Tumors with young cancer cells
featuring higher rates of mitotic activity, appear to grow faster.

This behavior comes in agreement with Figure 5 showing that young cells are
primarily responsible for the tumor’s expansion.

Figure 6: Evolution of $R_{mean}$ for various combinations of
$k_{m,yc}-k_{m,mc}$ values.

Figure 7: Evolution of average cancer cells volume fraction in the area of the
initial cancerous seed ($\sqrt{x^{2}+y^{2}}\leq R_{o}=1$) for different
combinations of $k_{m,yc}$ and $k_{m,mc}$ values.

It is then to be expected from what is shown in Figure 5 that tumors with more
active mature cells will show a delay in the formation of the necrotic core
compared to the rest. Indeed, as shown in Figure 7, these tumors exhibit a
significant delay in the necrotic core formation. Furthermore, Figure 7
illustrates the evolution of cancer cells volume fraction in the area of the
initial seed ($r=\sqrt{x^{2}+y^{2}}\leq 1$). Cells in the tumor’s interior
evade necrosis for longer periods of time, specifically in cases where the
overall mitosis rate has relatively lower value and cancer cells do not starve
their environment from nutrients at the same rate. As expected from the above
observations, the fastest growing tumor is the one with young cells featuring
more intense mitotic activity. A more detailed view of the tumor’s morphology
is presented in Figure 8 for three representative scenarios: (a) equal mitosis
rate constant values for young and mature cells, (b) higher mitosis rate
constant for younger cells, and (c) lower mitosis rate constant for younger
cells. The tumor with $k_{m,yc}=2.6$ and $k_{m,mc}=1.3$ grows faster compared
to the case: $k_{m,yc}=1.5-k_{m,mc}=3.0$, and develops earlier a necrotic
core.

Figure 8: Spatial volume fraction distribution for tumor cells
($\theta_{yc}+\theta_{mc}$) with (a) $k_{m,yc}=k_{m,mc}=2.0$ (top), (b)
$k_{m,yc}=2.6-k_{m,mc}=1.3$ (middle) and (c) $k_{m,yc}=1.5-k_{m,mc}=3.0$
(bottom) for $t=50,100$ and $150$.

Finally, we quantify the aggressiveness of the growing tumors for different
combination values of mitosis rate constants, by computing the total surface
covered by cancer cells at different time instances. In particular, we compute
the surface integral of the cancer cell volume fractions on the computational
domain over the surface of the said domain, $S=\pi R_{tissue}^{2}$:

$a\%=\frac{100\%}{\pi
R_{tissue}^{2}}\iint_{S}\left(\theta_{yc}+\theta_{mc}\right)\mathrm{d}x\mathrm{d}y.$
(50)

Figure 9 shows the evolution of $a\%$ for three representative cases of
mitosis constant rate values for young and mature cancer cells, and reveals
the higher aggressiveness of tumors with young cells exhibiting the highest
mitosis rate.

Figure 9: Temporal evolution of the domain’s area percentage covered by cancer
cells for tumor with equivalently proliferating cancer cells sub-populations
($k_{m,yc}=k_{m,mc}=2.0$), tumor with more proliferative young cells
($k_{m,yc}>k_{m,mc}$), and tumor with more proliferative mature cancer cells
(($k_{m,yc}<k_{m,mc}$).

###  3.3 Effect of cellular death rate heterogeneity

In the present paragraph, we present the effect of the death rate constant of
young and mature cancer sub-populations. Here, the ratio of death rate
constants $k_{d,yc}$ and $k_{d,mc}$, for young and mature cancer cells,
respectively, is kept constant. $\frac{k_{d,yc}}{k_{d,mc}}=2$ for tumors where
young cells are more susceptible to necrosis and
$\frac{k_{d,yc}}{k_{d,mc}}=\frac{1}{2}$ for tumors with young cancer cells
being more resilient compared to mature cells. $R_{mean}$ is calculated for a
series of simulations covering a range of $k_{d,yc}$ and $k_{d,mc}$
combinations satisfying both ratios. The results of these calculations is
presented in Figure 10. Tumors with more resilient young cells show a tendency
to expand at a higher rate, however one can observe that for the studied cases
the differences are not as apparent as in the results presented in the
previous paragraph, where the effect of mitosis rate heterogeneity is
examined.

Figure 10: Evolution of $R_{mean}$ for various combinations of
$k_{d,yc}-k_{d,mc}$ values.

Despite the fact that the growth rates do not show significant differences for
the different combinations of death rate constants, the examination of the
resulting tumor morphology shows interesting dissimilarities. Figure 11
illustrates snapshots of the cancer cell volume fraction for different death
rates of young and mature cancer cells. All three tumors grow at a similar
rate; however, significant differences in terms of the formation of a necrotic
core can be observed. In particular, tumors with more resilient mature cancer
cells tend to form a necrotic core of smaller size compared to the rest of the
studied cases at $t=150$. This comes in agreement with the observation that
mature cells are primarily focused in replenishing the dead cells in the inner
region of the tumor.

Figure 11: Spatial volume fraction distribution for tumor cells with (a)
$k_{d,yc}=k_{d,mc}=0.075$ (top), (b) $k_{d,yc}=0.1-k_{d,mc}=0.05$ (middle) and
$k_{d,yc}=0.06-k_{d,mc}=0.12$ (bottom) for $t=50,100$ and $150$.

Finally, we quantify the aggressiveness of the developing tumor, by comparing
the fraction of surface that is covered by cancer cells, $a\%$, using Equation
(50). The results are presented in Figure 12 for three representative
scenarios: (a) equal death rate constants, (b) more resilient mature cells,
and (c) more resilient young cells confirming our previous ascertainment that
the growth rates of the heterogeneous tumors are similar when varying the
death rate constants of young and mature cancer cells.

Figure 12: Temporal evolution of the domain’s area percentage, $a\%$, covered
by cancer cells for heterogeneous tumors with $k_{d,yc}=k_{d,mc}=0.075$,
$k_{d,yc}=0.1-k_{d,mc}=0.05$ and $k_{d,yc}=0.06-k_{d,mc}=0.12$.

###  3.4 Effect of material properties

In this paragraph, we assess the tumor growth dependency on the material
parameters of cancerous fluid phase; in particular, we examine the effect of
the viscosity and drag coefficients values for the mature cancer cell sub-
population. This assessment is performed by computing $R_{mean}$, which is
defined in paragraph 3.2. Figure 13(a) depicts the growth rate of tumors for
different values of the drag coefficients, $d_{mc,j}$, in the interphase
between mature cancer cells and the rest of the fluid phases
($j=h,yc,yv,mv,ecm$). By increasing the interphase drag results in a delay in
the tumor growth; such an effect is to be expected since increases in
$d_{mc,j}$ bolster the force, which opposes to the cells expansion movement.
Figure 13(b) shows the evolution of $R_{mean}$ for various values of the
dynamic viscosity $\mu_{mc}$ of mature cancer cells. Here, the trend of tumor
growth rate with changes in $\mu_{mc}$ is not as clear as in the case of
interphase drag. When the mature cancer fluid phase is behaving as a highly
viscous one, an increase in viscosity does not yield a clear trend in the
obtained tumor dynamics. The highest growth rate is attained for
$\mu_{mc}=50$, whereas the lowest growth rate for $\mu_{mc}=10$. The
simulation of the highest studied viscosity value $\mu_{mc}=200$ results the
second lowest tumor growth rate (exceeds only the $\mu_{mc}=10$ case). In
addition, by observing the equal growth rates at the initial stages of the
simulations for different interphase drag coefficients and viscosity values,
one can infer that the time needed by the tumor to surpass natural cell
density and start expanding is not affected by the mature cells material
properties.

Figure 13: Evolution of $R_{mean}$ for various values of mature cancer cells
(a) drag coefficients $d_{mc,j}$ (left) and (b) dynamic viscosity $\mu_{mc}$
(right).

We also study the effect of the cellular material properties of mature
cancerous sub-population on the tumor’s morphology, and in particular the size
of the developed necrotic core. Figure 14 presents the average volume fraction
of cancer cells within a radius $r=\sqrt{x^{2}+y^{2}}\leq 1$ from the origin
of the tissue domain. For values of $d_{mc,j}$, ranging from $1$ to $10$,
there is no clear effect on the time of formation of a necrotic region. By
further increasing its value to $d_{mc,j}=100$ one can observe a delay in the
formation of the necrotic core. Again, an increase of the viscosity of mature
cancer cells fluid phase does not seem to alter significantly the time of
formation, as well as the size of the necrotic core.

Figure 14: Average cancer cells volume fraction in the area of the initial
cancerous seed ($r=\sqrt{x^{2}+y^{2}}\leq 1$) for different values of
$d_{mc,j}$ (left) and $\mu_{mc}$ (right).

##  4 Discussion and Conclusions

In this work, we present a continuum-level computational model, which
simulates the complex and unstable environment of a growing heterogeneous
malignant tumor by treating the host tissue as a multi-phase fluid mixture.
The model features flexibility, since it can be easily modified to incorporate
additional elements of heterogeneity, with relatively low computational
requirements. In particular, the present study adopts the modeling concept
presented in 47, and incorporates the effect of ITH by introducing sub-
populations of cancer cells, which represent different stages of cellular
maturity. In addition, the microvasculature is categorized into mature blood
vessels and young sprouts equipped with different levels of ability to
withstand vascular occlusion.

The resulting system of partial differential equations formulating mass and
momentum balances for the different phases of the mixture is solved with the
Finite Elements Method utilizing the commercial software Comsol Multiphysics
®. By observing the dynamics of the different sub-populations involved in the
model, the young cancer cells are mostly responsible for the expansion of the
tumor; during the course of simulations the young cells form a protective
layer shielding the inner layers of the tumor, where mature cells are more
abundant. Meanwhile, mature cells move towards the interior of the tumor due
to the pressure exerted by the cells in the proliferating rim.

The concept of ITH was further explored by testing different combination
values of kinetic parameters so as to determine whether a tumor with
heterogeneous properties proliferates at a different rate compared to a tumor
whose kinetic parameters are uniform throughout its cancer cell sub-
populations. Indeed, our simulations reveal that introducing heterogeneity in
the cellular mitosis rate holds a crucial role in the resulting dynamics of
the tumor in question. When young cancer cells exhibit higher mitotic activity
compared to mature cells, the overall dynamics of the growing tumor are
considerably faster as opposed to cases where the mature cells are considered
to proliferate at higher rates. This observation can be considered in
congruence with the approach of differentiation therapy for malignant tumors,
a therapy based on forcefully differentiating cancer cells into a specific
mature state 33, 34. We adopt a similar approach in order to associate ITH
with cellular death rates. In this case, we observe again changes in the
macroscopic tumor dynamics; however, the computed growth rates for different
combinations of cellular death rates are not as significant as in the case of
modifying their proliferation rate.

We believe that our study is amenable on further modifications in order to
assess various topics related to ITH, and its implications to therapeutic
strategies. A possible enhancement to our model can be made by substituting
the existing kinetics in order to model the cellular metabolism to incorporate
a larger array of phenomena, such as the interaction of the cells with
different chemicals and nutrients or separating apoptosis from necrosis. In
addition, our model can be extended to incorporate possible treatments and
therapeutic strategies, including chemotherapy, radiotherapy, virotherapy or
hyperthermic techniques. Apart from evaluating the efficacy of these
techniques individually, their synergy with certain chemical agents or
elements of the tumor’s micro-environment can also be simulated (e.g.,
interaction with macrophages). Our modeling approach to sub-populate cancer
cells depending on their state of maturity, can be trivially modified in order
to incorporate the effect of differentiation therapies 33, 34, with the
inclusion of appropriate chemical agents that induce the maturation process.
Furthermore, the inclusion of cancer stem cells (CSC) as being one of the
cancerous phases during the study of a therapeutic path is an exciting new
prospect since, CSC are often responsible for provoking cancer relapse due to
their –at least partial- immunity to certain therapies 35, 48. Another
possible extension of the present model is to incorporate the effect of
macrophages, i.e., the immune cells most common in the immediate tumor
environment 49. Cancer cells have the ability to secrete chemical compounds
that mediate macrophage chemotaxis, i.e., direct the movement of macrophages
based on these chemicals’ concentration gradient.

### Conflict of interest

The authors declare no potential conflict of interests.

## References

  * 1 Fisher R, Pusztai L, Swanton C. Cancer heterogeneity: implications for targeted therapeutics. Br J Cancer 2013; 108(3): 479–485.  doi: 10.1038/bjc.2012.581
  * 2 Marusyk A, Polyak K. Tumor heterogeneity: causes and consequences. Biochim Biophys Acta Rev Cancer 2010; 1805(1): 105–117.  doi: 10.1016/j.bbcan.2009.11.002
  * 3 Marusyk A, Almendro V, Polyak K. Intra-tumour heterogeneity: a looking glass for cancer?. Nat Rev Cancer 2012; 12(5): 323–334.  doi: 10.1038/nrc3261
  * 4 Loeb LA. Human cancers express mutator phenotypes: origin, consequences and targeting. Nat Rev Cancer 2011; 11(6): 450–457.  doi: 10.1038/nrc3063
  * 5 Buck A, Aichler M, Huber K, Walch A. Chapter Five - In Situ Metabolomics in Cancer by Mass Spectrometry Imaging. In: Drake RR, McDonnell LA. , eds. Applications of Mass Spectrometry Imaging to Cancer. 134 of Advances in Cancer Research. Academic Press. 2017 (pp. 117 - 132)
  * 6 Piraino SW, Thomas V, O’Donovan P, Furney SJ. Mutations: Driver Versus Passenger. In: Boffetta P, Hainaut P. , eds. Encyclopedia of Cancer (Third Edition)Oxford: Academic Press. third edition ed. 2019 (pp. 551-562)
  * 7 Losic B, Craig AJ, Villacorta-Martin C, et al. Intratumoral heterogeneity and clonal evolution in liver cancer. Nat Commun 2020; 11(1): 1–15.  doi: 10.1038/s41467-019-14050-z
  * 8 Gay L, Baker AM, Graham TA. Tumour cell heterogeneity. F1000Res 2016; 5.  doi: 10.12688/f1000research.7210.1
  * 9 McGranahan N, Swanton C. Cancer evolution constrained by the immune microenvironment. Cell 2017; 170(5): 825–827.  doi: 10.1016/j.cell.2017.08.012
  * 10 Jamal-Hanjani M, Wilson GA, McGranahan N, et al. Tracking the evolution of non–small-cell lung cancer. N Engl J Med 2017; 376(22): 2109–2121.  doi: 10.1056/nejmoa1616288
  * 11 Nowell P. The clonal evolution of tumor cell populations. Science 1976; 194(4260): 23–28.  doi: 10.1126/science.959840
  * 12 Greaves M, Maley CC. Clonal evolution in cancer. Nature 2012; 481(7381): 306–313.  doi: 10.1038/nature10762
  * 13 Hinohara K, Polyak K. Intratumoral heterogeneity: more than just mutations. Trends Cell Biol 2019; 29(7): 569–579.  doi: 10.1016/j.tcb.2019.03.003
  * 14 Baslan T, Kendall J, Volyanskyy K, et al. Novel insights into breast cancer copy number genetic heterogeneity revealed by single-cell genome sequencing. eLife 2020; 9: e51480.  doi: 10.7554/elife.51480
  * 15 McDonald KA, Kawaguchi T, Qi Q, et al. Tumor heterogeneity correlates with less immune response and worse survival in breast cancer patients. Ann Surg Oncol 2019; 26(7): 2191–2199.  doi: 10.1245/s10434-019-07338-3
  * 16 Miller F, Miller B, Heppner G. Characterization of metastatic heterogeneity among subpopulations of a single mouse mammary tumor: heterogeneity in phenotypic stability.. Invasion Metastasis 1983; 3(1): 22–31.  PMID: 6677618.
  * 17 Zeng T, Dai H. Single-cell RNA sequencing based computational analysis to describe disease heterogeneity. Front Genet 2019; 10: 629\.  doi: 10.3389/fgene.2019.00629
  * 18 Koltai T. Tumor heterogeneity: An overlooked issue in cancer therapeutics. Integr Cancer Sci Ther 2017; 4(2): 1–4.  doi: 10.15761/icst.1000236
  * 19 Meacham CE, Morrison SJ. Tumour heterogeneity and cancer cell plasticity. Nature 2013; 501(7467): 328–337.  doi: 10.1038/nature12624
  * 20 Zhao B, Pritchard JR, Lauffenburger DA, Hemann MT. Addressing genetic tumor heterogeneity through computationally predictive combination therapy. Cancer Discov 2014; 4(2): 166–174.  doi: 10.1158/2159-8290.cd-13-0465
  * 21 Alvarez RF, Barbuto JA, Venegeroles R. A nonlinear mathematical model of cell-mediated immune response for tumor phenotypic heterogeneity. J Theor Biol 2019; 471: 42–50.  doi: 10.1016/j.jtbi.2019.03.025
  * 22 Dunn G, Bruce A, Ikeda H, Old L. J.; Schreiber, RD Cancer immunoediting: From immunosurveillance to tumor escape. Nat Immunol 2002; 3(11): 991–998.  doi: 10.1038/ni1102-991
  * 23 Dunn GP, Old LJ, Schreiber RD. The three Es of cancer immunoediting. Annu Rev Immunol 2004; 22: 329–360.  doi: 10.1146/annurev.immunol.22.012703.104803
  * 24 Harris LA, Beik S, Ozawa PM, Jimenez L, Weaver AM. Modeling heterogeneous tumor growth dynamics and cell–cell interactions at single-cell and cell-population resolution. Curr Opin Syst Biol 2019; 17: 24–34.  doi: 10.1016/j.coisb.2019.09.005
  * 25 Stamatelos SK, Bhargava A, Kim E, Popel AS, Pathak AP. Tumor ensemble-based modeling and visualization of emergent angiogenic heterogeneity in breast cancer. Sci Rep 2019; 9(1): 1–14.  doi: 10.1038/s41598-019-40888-w
  * 26 Lapidot T, Sirard C, Vormoor J, et al. A cell initiating human acute myeloid leukaemia after transplantation into SCID mice. Nature 1994; 367(6464): 645–648.  doi: 10.1038/367645a0
  * 27 Al-Hajj M, Wicha MS, Benito-Hernandez A, Morrison SJ, Clarke MF. Prospective identification of tumorigenic breast cancer cells. Proc Natl Acad Sci USA 2003; 100(7): 3983–3988.  doi: 10.1073/pnas.0530291100
  * 28 O’Brien CA, Pollett A, Gallinger S, Dick JE. A human colon cancer cell capable of initiating tumour growth in immunodeficient mice. Nature 2007; 445(7123): 106–110.  doi: 10.1038/nature05372
  * 29 Ricci-Vitiani L, Lombardi DG, Pilozzi E, et al. Identification and expansion of human colon-cancer-initiating cells. Nature 2007; 445(7123): 111–115.  doi: 10.1038/nature05384
  * 30 Singh SK, Hawkins C, Clarke ID, et al. Identification of human brain tumour initiating cells. Nature 2004; 432(7015): 396–401.  doi: 10.1038/nature03128
  * 31 Collins AT, Berry PA, Hyde C, Stower MJ, Maitland NJ. Prospective identification of tumorigenic prostate cancer stem cells. Cancer Res 2005; 65(23): 10946–10951.  doi: 10.1158/0008-5472.can-05-2018
  * 32 Tripathi S, Chakraborty P, Levine H, Jolly MK. A mechanism for epithelial-mesenchymal heterogeneity in a population of cancer cells. PLoS Comput Biol 2020; 16(2): e1007619.  doi: 10.1371/journal.pcbi.1007619
  * 33 Wang Y, Zhang Y, Qian C, et al. GSK3$\beta$/$\beta$-catenin signaling is correlated with the differentiation of glioma cells induced by wogonin. Toxicol Lett 2013; 222(2): 212–223.  doi: 10.1016/j.toxlet.2013.07.013
  * 34 Leszczyniecka M, Roberts T, Dent P, Grant S, Fisher PB. Differentiation therapy of human cancer: basic science and clinical applications. Pharmacol Ther 2001; 90(2-3): 105–156.  doi: 10.1016/s0163-7258(01)00132-2
  * 35 Sontag D, Miles DA, Uhanova J, Van Iderstine MG, Yang J, Minuk GY. The effects of hyperthermia on human hepatocellular carcinoma stem and mature cancer cells. Ann Hepatol 2020; 19(3): 265–268.  doi: 10.1016/j.aohep.2019.12.006
  * 36 Naeem E, Abid S. Ablative techniques in hepatocellular carcinoma treatment. Hepatoma Res 2018; 4: 18.  doi: 10.20517/2394-5079.2018.22
  * 37 Pries AR, Cornelissen AJ, Sloot AA, et al. Structural adaptation and heterogeneity of normal and tumor microvascular networks. PLoS Comput Biol 2009; 5(5): e1000394.  doi: 10.1371/journal.pcbi.1000394
  * 38 Folkman J. Tumor angiogenesis: therapeutic implications. N Engl J Med 1971; 285(21): 1182–1186.  doi: 10.1056/nejm197111182852108
  * 39 Bergers G, Benjamin LE. Tumorigenesis and the angiogenic switch. Nat Rev Cancer 2003; 3(6): 401–410.  doi: 10.1038/nrc1093
  * 40 Folkman J, Watson K, Ingber D, Hanahan D. Induction of angiogenesis during the transition from hyperplasia to neoplasia. Nature 1989; 339(6219): 58–61.  doi: 10.1038/339058a0
  * 41 Junttila MR, Sauvage dFJ. Influence of tumour micro-environment heterogeneity on therapeutic response. Nature 2013; 501(7467): 346–354.  doi: 10.1038/nature12626
  * 42 Jain RK. Normalizing tumor vasculature with anti-angiogenic therapy: a new paradigm for combination therapy. Nat Med 2001; 7(9): 987–989.  doi: 10.1038/nm0901-987
  * 43 Gazit Y, Berk DA, Leunig M, Baxter LT, Jain RK. Scale-invariant behavior and vascular network formation in normal and tumor tissue. Phys Rev Lett 1995; 75(12): 2428.  doi: 10.1103/physrevlett.75.2428
  * 44 Baish JW, Gazit Y, Berk DA, Nozue M, Baxter LT, Jain RK. Role of tumor vascular architecture in nutrient and drug delivery: an invasion percolation-based network model. Microvasc Res 1996; 51(3): 327–346.  doi: 10.1006/mvre.1996.0031
  * 45 Dewhirst MW, Navia IC, Brizel DM, Willett C, Secomb TW. Multiple etiologies of tumor hypoxia require multifaceted solutions. Clin Cancer Res 2007; 13(2): 375–377.  doi: 10.1158/1078-0432.ccr-06-2629
  * 46 Sorg BS, Moeller BJ, Donovan O, Cao Y, Dewhirst MW. Hyperspectral imaging of hemoglobin saturation in tumor microvasculature and tumor hypoxia development. J Biomed Opt 2005; 10(4): 044004.  doi: 10.1117/1.2003369
  * 47 Hubbard M, Byrne H. Multiphase modelling of vascular tumour growth in two spatial dimensions. J Theor Biol 2013; 316: 70–89.  doi: 10.1016/j.jtbi.2012.09.031
  * 48 Clevers H. The cancer stem cell: premises, promises and challenges. Nat Med 2011; 17(3): 313–319.  doi: 10.1038/nm.2304
  * 49 Unver N. Macrophage chemoattractants secreted by cancer cells: Sculptors of the tumor microenvironment and another crucial piece of the cancer secretome as a therapeutic target. Cytokine Growth Factor Rev 2019; 50: 13–18.  doi: 10.1016/j.cytogfr.2019.05.010


Proof —  This follows from Lemma 3.31, Theorem 4.32 and Theorem 4.41. $\hbox
to0.0pt{$\sqcap$\hss}\sqcup$

We have already scatteredly dealt with (local) central nilpotency and
hypercentrality in the previous results. This was unavoidable due to the many
connections they have with all kinds of properties of a brace. But there is
more to be said. For example, we saw in Theorem 4.19 that for a locally
supersoluble brace $B$, the concepts of locally centrally-nilpotent ideal and
locally $B$-nilpotent ideal coincide. Of course, this also holds for
hypercyclic braces, but in this latter case we can add something more.

###### Theorem 4.43

Let $B$ be a hypercyclic brace, and let $I$ be a (non-zero) locally centrally-
nilpotent ideal of $B$. Then $I$ is $B$-hypercentral.

Proof —  Since $B$ is hypercyclic, we can find a non-zero ideal $J$ of $B$
such that

  * •

$J\leq I$, and

  * •

either $J\leq\operatorname{Soc}(B)$ and $(J,+)$ is an infinite cyclic group,
or $|J|$ is prime.

We claim that $J\\!\leq\\!\zeta(I)$. Since $I$ is locally centrally-nilpotent,
we have that $L/M\leq\zeta(I/M)$ whenever $L/M$ is a chief factor of $I$ (see
[5], Theorem 4.6). Thus, if $|J|$ is prime, then $J\leq\zeta(I)$. On the other
hand, if $J$ is infinite, then $J/J^{p,+}\leq\zeta(B/J^{p,+})$ for every prime
$p$, so clearly $J\leq\zeta(B)$ also in this case. This shows that
$\zeta_{1}(I)_{B}\neq\\{0\\}$. Finally, an easy transfinite induction shows
that $I=\zeta_{\alpha}(I)_{B}$ for some ordinal $\alpha$, and completes the
proof of the statement. $\hbox to0.0pt{$\sqcap$\hss}\sqcup$

###### Corollary 4.44

Let $B$ be a hypercyclic brace. If $I$ is a hypercentral ideal of $B$, then
$I$ is $B$-hypercentral.

###### Corollary 4.45

Let $B$ be a hypercyclic brace. If $B$ is locally a finite brace of prime
power order, then $B$ is hypercentral.

Proof —  This follows from Theorem 4.43 and Theorem 3.7. $\hbox
to0.0pt{$\sqcap$\hss}\sqcup$

It has been proved in [5] that the join of finitely many $B$-hypercentral
ideals of a brace $B$ is still $B$-hypercentral. As a consequence of the above
results, we can say that in a hypercyclic brace, the join of arbitrarily many
hypercentral ideals is hypercentral.

###### Corollary 4.46

Let $B$ be a hypercyclic brace. Then the join of all hypercentral ideals of
$B$ is hypercentral, and coincides with the Hirsch–Plotkin ideal of $B$.

Proof —  This follows from Corollary 4.20, Theorem 4.43, and the fact that
every hypercentral brace is locally centrally-nilpotent. $\hbox
to0.0pt{$\sqcap$\hss}\sqcup$

In case of a locally supersoluble brace $B$ whose additive group is torsion-
free, we can say something more about the quotient of $B$ by the
Hirsch–Plotkin ideal.

###### Theorem 4.47

Let $B$ be a locally supersoluble brace, and let $H$ be the Hirsch–Plotkin
ideal of $B$. If $U_{2}^{+}(B)=\\{0\\}$, then $U_{2}^{+}(B/H)=\\{0\\}$.

Proof —  Suppose $U$ is not locally centrally-nilpotent, where
$U/H:=U_{2}^{+}(B/H)$. Then there is a finitely generated subbrace $F$ of $U$
that is not locally centrally-nilpotent. Consequently $|F/F\cap H|$ is odd, so
Theorem 3.45 yields that $F$ is centrally nilpotent, a contradiction. $\hbox
to0.0pt{$\sqcap$\hss}\sqcup$

If a locally supersoluble brace $B$ contains a large $B$-hypercentral ideal,
then the brace is bound to be hypercyclic.

###### Theorem 4.48

Let $B$ be a locally supersoluble brace, and let $I$ be a $B$-hypercentral
ideal of $B$ such that $B/I$ is finitely generated. Then $B$ is hypercyclic.

Proof —  Let $Z=\zeta(I)_{B}$, and let $x_{1},\ldots,x_{n}$ generating $B$
modulo $I$. If $a\in Z$, then $C=\langle a,x_{1},\ldots,x_{n}\rangle$ is
supersoluble, which means that $\langle a\rangle^{C}\leq I$ contains a non-
zero ideal $J$ of $C$ such that either $J\leq\operatorname{Soc}(C)$ and
$(J,+)$ infinite cyclic, or $|J|$ is a prime. Since $B=IC$, we have that
$J\trianglelefteq B$. Then the result follows by Theorem 4.4. $\hbox
to0.0pt{$\sqcap$\hss}\sqcup$

We end this paper by noting the following analog of Lemma 3.43 (that could be
of an independent interest), the proof being the same.

###### Theorem 4.49

Let $B$ be a hypercyclic brace, and let $I$ be a hypercentral ideal such that
$(I,+)$ is torsion-free. Then there is an ascending chain

$\\{0\\}=I_{0}\leq I_{1}\leq\ldots I_{\alpha}\leq I_{\alpha+1}\leq\ldots
I_{\lambda}=I$

of ideals of $B$ such that, for all $\beta<\lambda$,
$I_{\beta+1}/I_{\beta}\leq\operatorname{Soc}(B/I_{\beta})\cap\zeta(I/I_{\beta})$
and $(I_{\beta+1}/I_{\beta},+)$ is infinite cyclic.

##
5
Acknowledgements

We wish to thank the referee for their useful comments and remarks that really
improved the paper.

## References

  * [1]
  * [2] D. Bachiller: “Solutions of the Yang–Baxter equation associated to skew left braces, with applications to racks”, J. Knot Theory Ramifications 27 (2018), no. 8, 1850055, 36pp.
  * [3] D. Bachiller – F. Cedo – E. Jespers: “Solutions of the Yang–Baxter equation associated with a left brace”, J. Algebra 463 (2016), 80–102.
  * [4] A. Ballester-Bolinches – R. Esteban-Romero: “Triply factorised groups and the structure of skew left braces”, Commun. Math. Stat., 10 (2022), 353–370.
  * [5] A. Ballester-Bolinches – R. Esteban-Romero – M. Ferrara – V. Pérez-Calabuig – M. Trombetti: “Central nilpotency of left skew braces and solutions of the Yang–Baxter equation”, submitted; ArXiv:2310.07474v2.
  * [6] A. Ballester-Bolinches – R. Esteban-Romero – M. Ferrara – V. Pérez-Calabuig – M. Trombetti: “The word problem and the conjugacy problem for the structure skew brace of a solution of the Yang–Baxter Equation are solvable”, submitted.
  * [7] A. Ballester-Bolinches – R. Esteban-Romero – V. Pérez-Calabuig: “A Jordan–Hölder theorem for skew left braces and their applications to multipermutation solutions”, Proc. Roy. Soc. Edinburgh Sect. A (2023), to appear; doi:10.1017/prm.2023.37.
  * [8] A. Ballester-Bolinches – R. Esteban-Romero – P. Jiménez-Seral – V. Pérez-Calabuig: “Soluble skew left braces and soluble solutions of the Yang-Baxter solution”, submitted; ArXiv:2304.13475v2.
  * [9] A. Ballester-Bolinches – M. Ferrara – V. Perez-Calabuig – M. Trombetti: “A note on right-nil and strong-nil skew braces”, submitted; ArXiv:2310.11123.
  * [10] R. Baer: “Classes of finite groups and their properties”, Illinois J. Math. 1 (1957) 115–187.
  * [11] R. Baer: “Überauflösbare Gruppen”, Abh. Math. Sem. univ. Hamburg 23 (1957), 11–28.
  * [12] R.J. Baxter: “Partition function of the eight-vertex lattice model”, Ann. Physics 70, No. 1 (1972), 193–228.
  * [13] M. Bonatto – P. Jedlička: “Central nilpotency of skew braces”, J. Algebra Appl. (2022), to appear; doi:10.1142/S0219498823502559.
  * [14] M. Castelli – F. Catino – G. Pinto: “Indecomposable involutive set-theoretic solutions of the Yang–Baxter equation”, J. Pure Appl. Algebra, 223, No. 10 (2019), 4477–4493.
  * [15] F. Cedó – J. Okniński: “Constructing finite simple solutions of the Yang–Baxter equation”, Adv. Math., 391 (2021), Paper No. 107968, 39 pp.
  * [16] F. Cedó – J. Okniński: “Indecomposable solutions of the Yang–Baxter equation of square-free cardinality”, Adv. Math. 430 (2023), Paper No. 109221, 26pp.
  * [17] F. Cedó – E. Jespers – L. Kubat – A. Van Antwerpen – C. Verwimp: “On various types of nilpotency of the structure monoid and group of a set-theoretic solution of the Yang–Baxter equation”, J. Pure Appl. Algebra 227 (2023), 107194.
  * [18] F. Cedó – A. Smoktunowicz – L. Vendramin: “Skew left braces of nilpotent type”, Proc. Lond. Math. Soc. (3) 118 (2019), no. 6, 1367–1392.
  * [19] I. Colazzo – M. Ferrara – M. Trombetti: “On derived-indecomposable solutions of the Yang–Baxter equation”, Publ. Mat. (2023); ArXiv:2210.08598v2.
  * [20] V. Drinfel’d: “Hopf Algebras and the Yang–Baxter Equation”, in Yang–Baxter Equation in Integrable Systems, World Scientific (1990), 264–268.
  * [21] P. Etingof – T. Schedler – A. Soloviev: “Set theoretical solutions to the quantum Yang–Baxter equation, Duke Math. J. 100 (1999), 169–209.
  * [22] L. Faddeev – N. Reshetikhin – L. Takhtajan: “Quantum groups”, in Braid Group, Knot Theory and Statistical Mechanics, World Scientific, Teaneck, NJ (1989), 97–110.
  * [23] N. Fuster-Corral: “Left Braces and the Yang–Baxter Equation”, Ph.D. Thesis, Universitat de Valencia, Spain (2021).
  * [24] The GAP Group. GAP – Groups, Algorithms, and Programming, Version 4.12.2, 2022\. http://www.gap-system.org.
  * [25] T. Gateva-Ivanova: “Set-theoretic solutions of the Yang-Baxter equation, braces and symmetric groups”, Adv. Math. 338 (2018), 649–701.
  * [26] F. de Giovanni – M. Trombetti: “Infinite minimal non-hypercyclic groups”, J. Algebra Appl. 14 (2015), 1550143, 15pp.
  * [27] T. Gateva-Ivanova – P. Cameron. “Multipermutation solutions of the Yang–Baxter equation”, Comm. Math. Phys. 309, No. 3 (2012), 583–621.
  * [28] L. Guarnieri – L. Vendramin: “Skew-braces and the Yang-Baxter equation”, Math. Comp. 86 (2017), 2519–2534.
  * [29] P. Hall: “Some sufficient conditions for a group to be nilpotent”, Illinois J. Math. 2 (1958), 787–801.
  * [30] E. Jespers – L. Kubat – A. Van Antwerpen – L. Vendramin: “Radical and weight of skew braces and their applications to structure groups of solutions of the Yang–Baxter equation”, Adv. Math. 385 (2021), Paper No. 107767, 20pp.
  * [31] E. Jespers – A. Van Antwerpen – L. Vendramin: “Nilpotency of skew braces and multipermutation solutions of the Yang-Baxter equation”, Comm. Contemp. Math. 25, No. 09 (2023), 2250064.
  * [32] M. Newell – M. Trombetti: “A local approach to stability groups”, to appear.
  * [33] J. Orza: “A construction of the free skew-brace”; ArXiv:2002.12131.
  * [34] D.J.S. Robinson: “Finiteness Conditions and Generalized Soluble Groups”, Springer, Berlin (1972).
  * [35] D.J.S. Robinson: “A Course in the Theory of Groups”, Springer, New York (1996).
  * [36] W. Rump: “Braces, radical rings, and the quantum Yang–Baxter equation”, J. Algebra 307 (2007), 153–170.
  * [37] L. Stefanello – S. Trappeniers: “On bi-skew braces and brace blocks”, J. Pure Appl. Algebra 227 (2023), 107295, 22pp.
  * [38] A. Smoktunowicz: “A note on set-theoretic solutions of the Yang–Baxter equation”, J. Algebra 500 (2018), 3–18.
  * [39] M. Trombetti: “The structure skew brace associated with a finite non-degenerate solution of the Yang–Baxter equation is finitely presented”, Proc. AMS (2023); doi:10.1090/proc/16566.
  * [40] L. Vendramin: “Extensions of set-theoretic solutions of the Yang–Baxter equation and a conjecture of Gateva-Ivanova”, J. Pure Appl. Algebra 220 (2016), 2064–2076.
  * [41] L. Vendramin and O. Konovalov: “YangBaxter: combinatorial solutions for the Yang-Baxter equation”, December 2022. version 0.10.2, https://gap-packages.github.io/YangBaxter/.
  * [42] C.N. Yang: ”Some exact results for the many-body problem in one dimension with repulsive delta-function interaction”, Phys. Rev. Lett. 19, No. 23 (1967), 1312–1315.
  * [43] B.A.F. Wehrfritz: “Supersoluble and locally supersoluble linear groups”, J. Algebra 17 (1971), 41–58.

Adolfo Ballester-Bolinches, Ramon Esteban-Romero, Vicent Pérez-Calabuig

Departament de Matemàtiques

Universitat de València, Dr. Moliner, 50, 46100 Burjassot, València (Spain)

<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>

Maria Ferrara

Dipartimento di Matematica e Fisica

Università degli Studi della Campania “Luigi Vanvitelli”

viale Lincoln 5, Caserta (Italy)

e-mail<EMAIL_ADDRESS>

Marco Trombetti

Dipartimento di Matematica e Applicazioni “Renato Caccioppoli”

Università degli Studi di Napoli Federico II

Complesso Universitario Monte S. Angelo

Via Cintia, Napoli (Italy)

e-mail<EMAIL_ADDRESS>