# Super-Linear Growth Reveals the Allee Effect in Tumors

Youness Azimzade<EMAIL_ADDRESS>Department of Physics, University of
Tehran, Tehran 14395-547, Iran  Abbas Ali Saberi<EMAIL_ADDRESS>Department of Physics, University of Tehran, Tehran 14395-547, Iran  Institut
für Theoretische Physik, Universitat zu Köln, Zülpicher Strasse 77, 50937
Köln, Germany  Robert A. Gatenby<EMAIL_ADDRESS>Cancer Biology
and Evolution Program, Moffitt Cancer Center, Tampa, FL, USA.  Integrated
Mathematical Oncology Department, Moffitt Cancer Center, Tampa, FL, USA.
Diagnostic Imaging Department, Moffitt Cancer Center, Tampa, FL, USA

###### Abstract

Abstract: Integrating experimental data into ecological models plays a central
role in understanding biological mechanisms that drive tumor progression where
such knowledge can be used to develop new therapeutic strategies. While the
current studies emphasize the role of competition among tumor cells, they fail
to explain recently observed super-linear growth dynamics across human tumors.
Here we study tumor growth dynamics by developing a model that incorporates
evolutionary dynamics inside tumors with tumor-microenvironment interactions.
Our results reveal that tumor cells’ ability to manipulate the environment and
induce angiogenesis drives super-linear growth —a process compatible with the
Allee effect. In light of this understanding, our model suggests that for
high-risk tumors that have a higher growth rate, suppressing angiogenesis can
be the appropriate therapeutic intervention.

##  I Introduction

Tumors comprise a huge number of cells that evolve and expand through a
diverse range of interactions between themselves and their environment. Such a
complex structure resembles an ecosystem in which different evolutionary
processes and ecological interactions exist simultaneously Korolev et al.
(2014). A relevant, still open problem is to find what processes are more
relevant and play the dominant role compared to others Turajlic et al. (2019).
Finding the answer to this question is of theoretical interest and is
necessary for designing new promising therapies Bozic et al. (2013); Gatenby
and Brown (2020). Despite numerous research carried on in this field, possible
cooperation and its role in the growth and clinical implications have remained
less understood, partially due to emphasis on competition between sub-
populations in tumor progression Wang et al. (2014); Parker et al. (2020);
Zahir et al. (2020); Lakatos et al. (2020); Williams et al. (2020).

Tumors grow in initially healthy tissues with a limited oxygen supply rate.
Such a restriction imposes a limit (up to a few mm3) on growing tumors
Carmeliet and Jain (2000). To expand further, tumors need to increase oxygen
supply rate by inducing angiogenesis —a process so common and essential that
has been introduced as a ”Hallmark of Cancer” Hanahan and Weinberg (2000,
2011). During angiogenesis, each hypoxic cell secretes growth factors that
induce angiogenesis, leading to mutualism among tumor cells Archetti and
Pienta (2019). Cooperation to manipulate the environment in favor of the whole
population is analogous to the ”Allee Effect” in which individuals cooperate
to increase population growth rate Allee and Bowen (1932). While the Allee
effect has been suggested to play a role in tumor growth Korolev et al.
(2014); Böttger et al. (2015); Brown et al. (2016); Gatenby et al. (2019),
direct evidence for such a role has been scarce Johnson et al. (2019).

Tumor growth dynamics can reveal the biological processes in cellular scale
Rodriguez-Brenes et al. (2013); West and Newton (2019). In this regards, if
all cells can duplicate, for tumor volume (or cell number) one can write:
$\frac{dC_{T}}{dt}\propto C_{T}^{\beta}$ with $\beta=1$ which represents an
exponential growth Collins (1956). For tumors growing in a 3D environment
(assuming a well-fed tumor in which cells obey contact inhibition), only cells
living on the tumor border can duplicate. Respectively, considering the tumor
as a sphere with a radius of $r$, the number of cells on its surface is
proportional to $r^{2}$ and we have $\beta=2/3$ Mayneord (1932). Additional
limitations, such as oxygen shortage, can decrease this fraction even further
and finally set an upper limit on population size West et al. (2001). However,
recent analysis of different cancers types revealed that tumor growth rate
versus volume increases super-linearly as $\frac{dC_{T}}{dt}\propto
C_{T}^{\beta}$ with $\langle\beta\rangle=1.25$ Pérez-García et al. (2020).

The mechanism behind super-linear growth was attributed to the evolutionary
dynamic of competition between cell types. Based on this analogy, once a new
cell type with a higher duplication rate emerges, it starts to dominate the
tumor and we should observe an increase in tumor growth rate. The continuous
emergence of new cell types with higher fitness advantages and their
domination increases the tumor growth rate. This process can cause tumors to
grow faster than exponentially growing populations and respectively have
$\beta>1$ Durrett et al. (2010); Pérez-García et al. (2020). Since the
considerably small probability of driver mutation and small fitness advantages
they bring are unlikely to have such a drastic effect on $\beta$, the search
for the driving mechanism, which is common among different types of cancer,
continues. Additionally, tumors with higher growth rates were associated with
lower survival times Pérez-García et al. (2020), and finding the driver
mechanism may have therapeutic implications. It should be noted that the
emergence of new types with higher fitness advantages can increase the growth
rate and accelerate invasion Phillips et al. (2010); Benichou et al. (2012);
Hallatschek and Fisher (2014). However, it is not clear that if such
acceleration is sufficient to significantly increase $\beta$, particularly
because even extreme driver mutations can remain unnoticed due to their small
effect on overall growth Bozic et al. (2019).

Here, we developed a multiscale model that includes driver mutations and
fitness advantages they bring and cellular oxygen consumption and the ability
to induce angiogenesis. Our model suggests that the evolutionary dynamics of
competition between sub-populations do not lead to super-linear growth.
Instead, super-linear growth emerges during angiogenesis and as a result of
the change in the fraction of actively duplicating cells. Thus, super-linear
growth dynamics with $\beta>1$ is the direct evidence for mutualism among
tumor cells and the Allee effect across human tumors. Additionally, since
higher proliferation has come due to more prosperous angiogenesis, suppressing
angiogenesis seems to be a proper therapeutic target in high-risk tumors.

##  II Model

Mathematical modeling has been used to understand and describe different
aspects of cancer incidence and evolution over the last century Byrne (2010);
Altrock et al. (2015). Interestingly, mathematical description of tumor growth
predates our understanding of cancer as a genetic disease Mayneord (1932);
Schrek (1935). More recently, mechanistic models have been developed to
explain the underlying biological processes that drive tumor growth by
including the most fundamental properties of tumor cells such as proliferation
and diffusion Anderson and Quaranta (2008). One can consider individual cells
and their behavior and study dynamics of a collection of cells where
corresponding models are known as individual-based models (IBMs) Metzcar et
al. (2019); Azimzade and Saberi (2019). Alternatively, one can start from a
mesoscopic scale and define averaged variables and parameters for a group of
cells by taking a mean-field approach (such models are known as continuum). In
this approach, a collection of partial differential equations (PDEs) describe
densities of cellular populations and relevant chemicals Swanson et al.
(2003); Mandonnet et al. (2003); Azimzade et al. (2019). When the number of
cells is small, IBMs are preferred due to a more detailed description of
individual cells. However, tumors normally contain billions of cells where it
is impossible to use IBMs. On such a large scale, continuum models that use
PDEs are preferred. These two approaches can be integrated to generate hybrid
models that simulate larger cell numbers with more biological details Rejniak
and Anderson (2011); Jiménez-Sánchez et al. (2021). With similar assumptions
for population dynamics, IBMs and PDEs may not lead to the same results
Azimzade et al. (2020) and the connection between the two approaches is not
necessarily straightforward. For example, obtaining the corresponding PDE for
an IBM can be demanding and lead to various PDEs depending on parameter
choices for a particular IBM Johnston et al. (2017). A common approach to
ensure that results are generic and model-independent is to use different
continuum and IBM models alongside each other and compare their output, as we
will do in this work.

As our main model, we develop a hybrid model that integrates cells’ internal
features into external limitations across different scales. We first explain
the basic model that rules tumor cells’ dynamics and different aspects will be
added to the model as we proceed. To include the central role of metabolism in
cellular activities Buchakjian and Kornbluth (2010), we define an internal
energy $u$, which is the intensive quantity and is assigned to ”each cell”.
This parameter is regulated by a simple conservation rule as Scalerandi and
Sansone (2002); Azimzade et al. (2018):

$\frac{\partial u}{\partial t}=-a_{u}u+b_{u}n(r,t),$ (1)

where $n(r,t)$ is the oxygen density at the position of cell ($(r)$) and time
$t$, with $a_{u}$ and $b_{u}$ being positive constants related to energy
consumption and accumulation rate.

Figure 1:  (a) Linear (Darwinian) Evolution Model. We start with a single cell
type and during each duplication there is chance ($P_{m}$ which is estimated
previously Bozic et al. (2010)) for a driver mutation and appearance of more
aggressive cell type. Each driver mutation brings a fitness advantage, $\Delta
R$ Bozic et al. (2010). (b) Punctuated Evolution Model. In this model,
different cell types with different fitness advantages live together but there
is no transition between cell types. (c) Relation between tumor and oxygen
supply system. Tumor cells can increase oxygen supply rate (induce
angiogenesis) by $b_{n1}$.

We consider a 3D lattice with a unit size of 100 $\mu$m which, considering a
typical 10 $\mu$m diameter for cancer cells, has the capacity for 1000 cancer
cells. Cells live in these units and increase their biomass until they reach
the threshold to duplicate. The density of cells of type $i$ at $(r,t)$ is
defined by, $C_{i}(r,t)=\frac{\text{number of cells at}\;(r,t)}{\text{capacity
of each site}}$. Tumor cells consume the existing oxygen and the vascular
system supplies it for the tissue. As such, the evolution of oxygen is
regulated by:

$\displaystyle\frac{\partial n(r,t)}{\partial t}$ $\displaystyle=$
$\displaystyle D_{n}\nabla^{2}n(r,t)$ $\displaystyle-
a_{n}n(r,t)+b_{n0}-c_{n}C_{t}(r,t)n(r,t),$

where $D_{n}$ is oxygen diffusion constant (equal to $10^{-5}$ cm2/s), $c_{n}$
is the oxygen consumption rate which is reported to be $6.65\times 10^{-17}$
mol cell-1s-1 Anderson (2005), $a_{n}$ is oxygen decay rate and $b_{n}$ is
oxygen supply rate. $C_{t}$ is total density of different cell populations and
can be obtained as $C_{t}(r,t)=\sum^{g}_{i=1}C_{\rm i}(r,t)$ in which $g$ is
the number of cell types.

The evolution of the cells’ internal energy, $u$, depends on the local density
of oxygen through Eq. (7). If cells have access to enough oxygen, $u$
increases to $u_{p}$ and the cells start to proliferate. Cells that surpass
the internal energy threshold follow the standard Fisher-Kolmogorov-Petrovsky-
Pishkonov Murray (2003) as:

$\frac{\partial C(r,t)}{\partial t}=D_{c}\nabla^{2}C(x,y,t)+RC(r,t)(1-C(r,t))$
(3)

where $D_{c}$ is the diffusion coefficient for cells ($\approx 10^{-10}$ cm2/s
Anderson (2005)) and $R$ is duplication rate for cells ($\approx 10^{-5}$ s-1
Anderson (2005)). For an environment with enough oxygen, cells always have the
required internal energy to duplicate with the rate of $R$ and our model
automatically becomes the standard FKPP equation. Such models have been used
widely in literature to model tumor growth Swanson et al. (2003); Mandonnet et
al. (2003); Harpold et al. (2007). They consider two main aspects of cellular
behavior: duplication and diffusion. While these parameters are of crucial
importance in tumor progression, they are highly sensitive to environmental
properties such as oxygen density Hanahan and Weinberg (2000, 2011); Whiteside
(2008). Our model integrates essential features of cells into environmental
limitations by assigning a central role to cells’ metabolic states. We start
all simulations with an initial density of cells in the center of a defined
environment by setting $C(r,0)=C_{0}$ for $r<r_{0}$. In the following parts,
we will add more features to the model to incorporate different possibilities.

###  II.1 Evolutionary Dynamics

Here, we consider two main hypotheses that have been used to describe the
frequency of sub-populations in tumors, Linear (Darwinian) Evolution and
Punctuated Evolution Wang et al. (2014); Davis et al. (2017); Turajlic et al.
(2019); Williams et al. (2020). The linear evolution suggests that as the
tumor grows, new cell types emerge due to mutations and those with higher
fitness advantages become dominant, similar to what Darwinian evolution
suggests. However, some tumors show high diversity in initial steps, which
suggests that they have experienced an abrupt increase in mutations at early
stages. A similar process happens in ”Punctuated Equilibria,” where adaptation
occurs in small time intervals and newly adapted individuals rapidly expand
out of the niche and through the wider population Gould and Eldredge (1972).
Due to such similarity, punctuated tumor growth evolution is used where
different cell types appear at the early stages of tumor growth. Those with
higher fitness advantages become dominant as the tumor grows Turajlic et al.
(2019). It should be noted that other evolutionary scenarios are possible
Davis et al. (2017). Due to a large number of cells and their nonlinear
behavior (cell plasticity, for example), the actual dynamics can be more
complex Turajlic et al. (2019). Still, we used these two approaches to check
if evolutionary dynamics of competition between sup-populations can lead to
super-linear growth.

####  II.1.1 Linear Evolution

As tumors progress, new mutations appear Hanahan and Weinberg (2011). However,
the detailed mechanism and rate by which mutations happen vary over time and
across different cell types. A typical aspect of mutation is that it can be
highly correlated. Once a mutation appears in a specif part of DNA, it can
lead to instabilities and create a hot spot prone to further mutations Nesta
et al. (2020). On the other hand, mutations can bring different fitness
advantages Vermeulen et al. (2013); McFarland et al. (2014). However, instead
of single mutations with huge fitness advantages, it is more accepted that
they happen more frequently and gradually lead to variations in cellular
features Castro-Giner et al. (2015). Driver mutations, in general, can lead to
higher fitness advantages in different properties of cells. For example, they
can lead to higher duplication rates, higher migration rates, or more
effective nutrient consumption. However, this process can not go on forever
and the differences between fitness advantages are not huge. On the other
hand, higher fitness advantages in one aspect are normally associated with
loss of other properties Persi et al. (2020); Hausser and Alon (2020). A well-
known example of this case is the go-or-grow paradigm Hatzikirou et al. (2012)
that suggests cells decide to either migrate or duplicate faster and they can
not accomplish both tasks simultaneously. Despite these complexities, simple
mathematical models can effectively describe mutations and fitness advantages
they bring Waclaw et al. (2015); Paterson et al. (2020). As with any random
process, mutations are considered to follow a Poison distribution. Thanks to
this simplification, the average rate by which driver mutations appear and the
corresponding fitness advantages they bring can be calculated. Analysis of
experimental results revealed that driver mutations appear with the rate of
$3.4\times 10^{-5}$ per duplication and each diver mutation leads to an
average 0.4 percent increase in fitness advantage Bozic et al. (2010). Here we
study the effect of driver mutations on duplication and migration rates.

In linear evolution, we start with one cell type, $C_{1}$, by defining the
initial condition as $C_{1}(r,0)=C_{0}$ for $r<r_{0}$ ($C_{i}(r,0)=0$ for
$i>1$ and $C_{1}(r,0)=0$ for $r>r_{0}$). $C_{1}$ cells increase their internal
energy and duplicate. During each duplication, one daughter cell will have the
same type as the parent. The other daughter cell will have the same type with
probability of $1-P_{m}$ and belongs to the next cell type with the
probability of $P_{m}$ with $P_{m}=3.4\times 10^{-5}$ ( Under this assumption,
mutations incidence follows Poison distribution. To include a wider range of
possibilities, we include larger values for $P_{m}$ up to $10^{-1}$ which is
huge compared to any experimental data. ). To take to account such transition
between cell types, we need to replace Eq. (3) with a new version as below:

$\displaystyle\frac{\partial C_{\rm i}(r,t)}{\partial t}=D_{c}\nabla^{2}C_{\rm
i}(r,t)$ $\displaystyle+C_{\rm i}(r,t)\delta_{1i}R_{i}(1-P_{m})[1-C_{t}(r,t)]$
$\displaystyle+P_{m}(1-\delta_{1i})R_{i-1}C_{\rm i-1}(r,t)[1-C_{t}(r,t)]$ (4)

where $\delta_{ij}$ denotes the Kronecker delta, i.e., $\delta_{ij}=1$ for
$i=j$ and $0$ otherwise, with $1\leq i,j\leq g$. The second term on the right
side of Eq. (4) is the creation of each type of the cells due to duplication
subtracted by those that gain new driver mutations and change their type
during duplication. The third term represents the creation of the $i$th type
(for $i\neq 1$) of the cells from mutations of the prior type.

Each driver mutation leads to $0.4$ percent increase in fitness advantage
Bozic et al. (2010). One can take into account increase in fitness advantages
due to such mutations by setting $\Delta R=1-\frac{R_{i+1}}{R_{i}}$ where
$R_{i}$ and $R_{i+1}$ are the duplication rates of types $i$ and $i+1$,
respectively. Fitness advantages are considered to be small Heide et al.
(2018). However, to include a wider range in fitness advantages Castro-Giner
et al. (2015); Vermeulen et al. (2013); McFarland et al. (2014), we also
consider higher values for $\Delta R$ up to $10^{-1}$ which also is huge
compared to available experimental data. It should be noted that driver
mutations may change other properties in cells as well. For example, similar
to other populations, tumor cells may gain higher motility due to driver
mutations Benichou et al. (2012). In this regard, we also considered driver
mutations to increase diffusion constant by $\Delta D$.

####  II.1.2 Punctuated Evolution

For Punctuated Evolution, we consider different cell types to live together
and compete over limited sources. We set their initial values to be equal but
their duplication rate to be different as: $R_{i+1}=R_{i}\times(1+\Delta R)$
for $1\leq i\leq g-1$. Again, to include a wide range of variations in fitness
advantages, we will consider different values for $\Delta R$. Since in this
model mutations are considered to be negligible, the evolution of populations
follows the equation 4 with $P_{m}$=0. We consider 100 cell types (this number
does not change the outcome as long as the highest fitness advantage
difference remains in a reasonable range — we have checked up to 800 percent).
The frequency of subpopulations with different numbers of mutations is subject
to mutation rate and possible fitness advantages they bring Turajlic et al.
(2019); Hausser and Alon (2020). We should consider different possibilities
for the initial frequency of mutations. For the first case, as the initial
condition we set $C_{i}(r)=C_{0}$ $C_{i}(r,0)=C_{0}$ for $r<r_{0}$
($C_{i}(r,0)=0$ for $r>r_{0}$) and let the populations evolve based on
aforementioned rules. It is widely known for neutrally evolving populations
that population number decreases versus mutation frequency, $i$, as $1/i$
Williams et al. (2016). We can implement such a scenario by considering
$C_{i}(r,0)\propto C_{0}/i$ as the initial condition for $r<r_{0}$. These two
evolutionary models will be used to check if evolutionary dynamics of
competition between different cell types can lead to the emergence of super-
linear growth.

###  II.2 Angiogenesis

Angiogenesis occurs when starving cells start to secrete chemicals such as
Vessel Endothelial Growth Factor (VEGF) to induce endothelial cell duplication
and formation of vessels. Endothelial cells in response provide starving cells
with nutrients. A precise model can take into account the density of
endothelial cells, the density of starving cells, the rate by which chemicals
are secreted/up taken, as well as their temporal and spatial evolution
Scalerandi and Sansone (2002); Phillips et al. (2020). However, such an
approach leads to an excessive number of parameters Phillips et al. (2020);
Zheng and Sweidan (2018). For the sake of simplicity, we take a different
approach. Rationally, one can assume the angiogenesis (and subsequent increase
in oxygen supply rate) to be a monotonic function of tumor volume d’Onofrio
(2007). Additionally, oxygen supply rate can be a monotonic function of
vascular density. Putting together, a single parameter can effectively
describe how angiogenesis affects oxygen supply rate. As such, oxygen
evolution (previously described by Eq. (II)) will have an additional term as:

$\displaystyle\frac{\partial n(r,t)}{\partial t}$ $\displaystyle=$
$\displaystyle D_{n}\nabla^{2}n(r,t)$ $\displaystyle-
a_{n}n(r,t)+b_{n0}+b_{n1}C_{T}-c_{n}C_{t}(r,t)n(r,t),$

where $b_{n1}$ (in mol cell-1 s-1mm-3) depicts the increase in oxygen supply
as a result of angiogenesis. This approach is conceptually in line with the
previous studies on the dynamics of vasculature d’Onofrio (2007); Rieger and
Welter (2015); Tracqui (2009); Xu et al. (2016). However, contrary to what one
may naively assume, angiogenesis is not a linear process. A vascular tumor may
fail to grow larger due to failure in inducing more angiogenesis, a process
known as post-vascular dormancy Hahnfeldt et al. (1999); Naumov et al. (2008).

###  II.3 Solving the model

We have two kinds of inter-related variables that evolve over time, oxygen
density and cellular densities. Oxygen density is constantly changing due to
supply and consumption by cells. On the other hand, cellular density evolves
if cells have accumulated enough internal energy. Respectively, equations that
describe these variables are not coupled as one may see in conventional
approaches and allows us to solve equations separately. At each time step, we
solve oxygen density using Forward Time Centered Space (FTCS). For simplicity,
we explain the approach for the 1D version of Eq. (2) and extension to higher
dimensions is straightforward. Discretization of the reaction terms
(consumption, decay and supply) is straightforward. For Laplacian-like term in
a uniformly-space grid with a set of grid points $x_{i}=i\delta$ ($i=0\l,\pm
1\;,\pm 2,\cdots$), let $n(x_{i},t)=n_{i}(t)$, and use central difference-
discretization, $\mbox{\boldmath$\nabla$}n\approx(n_{i+1}-n_{i-1})/(2\delta)$,
where $j$ belongs to nearest neighbors of $i$ and $\delta$ is the distance
between $i$ and $i+1$, we obtain:

$\frac{\partial n_{i}}{\partial
t}=D_{n}(n_{i}-n_{i-1})+D_{n}(n_{i}-n_{i+1})-a_{n}n_{i}(t)+b_{n0}-c_{n}C_{t}(x_{i},t)n_{i}(t)\;.$
(6)

As mentioned before, we set unit size (and respectively distance between two
neighboring grid points) at $100\;\mu$m. For each grid site we solve for $n$
by writing:

$n_{i}(t+\Delta t)=n_{i}(t)+\Delta
t\sum_{j\in\\{i\\}}D_{n}[n_{j}(t)-n_{i}(t)]+\Delta
t\big{(}-a_{n}n_{i}(t)+b_{n0}-c_{n}C_{t}(x_{i},t)n_{i}(t)\big{)}\;,$

where $\Delta t$ is the time step. FTCS method typically converges to the
solution and is stable. Stability requires that $D\Delta t/(4\delta^{2})<0.5$,
suggesting $\Delta t\sim 10s$. As the initial condition for oxygen, we set
$n(r,t=0)=b_{n0}/a_{n}$. For cellular densities that surpass the energy
threshold, the corresponding FKPP equation was solved using FTCS method
Azimzade et al. (2019) and codes are available at
https://gitlab.com/YounessAzimzade/superlinearallee.

###  II.4 Parameters and Variables

To give a comprehensive view to readers, here we present all variables that we
have defined (Table I) and parameters we have used (Table II) in our model.

Table 1:  Variables defined in the model and their definitions and values.  Variable |  Definition
---|---
$u$ |  Internal energy of each cell
$n$ |  Oxygen density
$C$ |  Cellular density defined as $C_{i}(r,t)=\frac{\text{number of cells at}\;(r,t)}{\text{capacity of each site (1000)}}$
$C_{i}$ |  Density of type $i$
$C_{t}$ |  Total cellular density defined by $C_{t}(r,t)=\sum^{g}_{i=1}C_{\rm i}(r,t)$
$C_{T}$ |  Tumor volume $C_{T}(t)=\int C_{t}(r,t)dV$
$\frac{dC_{T}(t)}{dt}$ |  Tumor growth rate
Table 2:  Parameters used in the model and their definitions and values.  Parameter |  Definition |  Value
---|---|---
$a_{u}$ |  Energy consumption rate |  $0.1$ (estimated)
$b_{u}$ |  Energy accumulation rate |  $1$ (estimated)
$u_{p}$ |  Proliferation energy threshold |  $8$ (estimated)
$R$ |  Cellular duplication rate |  $0.05$ Haass et al. (2014)
$D_{C}$ |  Cell diffusion constant |  $10^{-10}$cm2/s Anderson (2005); Bray (2001)
$D_{n}$ |  Oxygen diffusion constant in tissue |  $10^{-5}$ cm2/s Anderson (2005)
$c_{n}$ |  Oxygen consumption rate |  $6.65\times 10^{-17}$ mol cell-1s-1 Anderson (2005)
$b_{n0}$ |  Oxygen supply rate |  variable - estimated
$b_{n1}$ |  Increase in oxygen supply rate |  variable - estimated
$a_{n}$ |  Oxygen decay rate |  10-4 \- estimated
$P_{m}$ |  Driver mutation rate |  $3.4\times 10^{-5}$ per division Bozic et al. (2010)
$\Delta R$ |  Fitness Advantage of each driver mutation |  $4\times 10^{-3}$ Bozic et al. (2010)

##  III Results

We first study the basic model and then will analyze how adding new aspects
affect the growth dynamics. To quantify growth dynamics, we measure tumor
volume, $C_{T}(t)=\int C_{t}(r,t)dV$ and the accumulated duplication of all
cells, $\frac{dC_{T}(t)}{dt}=\int\frac{\partial C_{t}(r,t)}{\partial t}dV$,
where the latter gives the growth rate for whole tumor. We first check two
possibilities: (a) Where the oxygen supply rate is high ($b_{n0}=10^{-4}$ mol
s-1mm-3) and all cells will obtain the requisite oxygen to duplicate. (b)
Where the oxygen supply rate is low ($b_{n0}=5\times 10^{-6}$ mol s-1mm-3) and
its shortage halt the growth when tumor volume approaches $\sim$2 mm3 — about
the maximum size tumors can grow without angiogenesis.

Figure 2: (a) Volume versus time for different oxygen supply rate ($b_{n0}$ in
mol s-1mm-3) in the half-logarithmic diagram. Both high and low oxygen supply
initially show exponential growth. (b) Growth rate versus volume for different
values of oxygen supply. Oxygen shortage halts growth and induces smaller
values for $\beta$. Limited oxygen supply forces an upper limit where the
population size can not surpass it.

For high oxygen supply rate, our model converges to a traditional 3D FKPP
equation. As such, initially all cells can duplicate and we observe an
exponential growth. Thus, one has $\frac{dC_{T}(t)}{dt}\sim C_{T}(t)$ and
tumor volume grows as $C_{T}(t)=C_{T}(0)e^{Rt}$ (see early times dynamics in
FIG. 2 (a) and FIG. 2 (b) where we have $\beta=1$). For a considerably large
tumor, only cells that are living on the surface will have a chance to
duplicate and fraction of actively duplication cells for a spherical tumor of
radius $r$ becomes proportional to $r^{2}$. With minor simplification, we
have: $\frac{dC_{T}(t)}{dt}\sim r^{2}$. Thus, in the log-log diagram we have
the $\frac{dC_{T}(t)}{dt}\sim C_{T}(t)^{\beta}$ with $\beta=2/3$, as FIG. 2
(b) shows.

In the second case, which is more realistic, considering a wide range of
evidence for hypoxia, the oxygen supply rate is lower and cells duplicate
slower. Yet, as FIG. 2 (a) shows, tumor grows exponentially at earlier stages
and we have $\beta=1$ (see FIG. 2 (b)). As the tumor grows larger, more and
more cells face a nutrient shortage, and duplication becomes slower. Limited
oxygen supply imposes an upper limit on population size. When tumor size
reaches the limit, growth stops and $\frac{dC_{T}(t)}{dt}$ approaches zero.

###  III.1 Evolutionary Dynamics

Here we check if evolutionary processes can lead to the emergence of super-
linear dynamics. Theoretically, emergence of new cell lines with higher
fitness advantages can accelerate the growth Durrett et al. (2010); Pérez-
García et al. (2020) (see FIG. 3 (a)).

####  III.1.1 Linear Evolution

We start the simulation with a single cell type and allow new cell types to
arise due to mutations. We include different values for $\Delta R$ and $P_{m}$
and assume that driver mutations increase duplication rate. As FIG. 3 (b)
shows tumor volume versus time, higher mutation rate (fitness advantage) lead
to faster growth. However, this change is very small for experimentally
relevant parameters. On the other hand, we need to check if this acceleration
in growth is big enough to cause super-linear growth. As an analysis of growth
dynamics in FIG. 3 (c) shows, even considerably large values of $\Delta R$ and
$P_{m}$ do not increase $\beta$.

Similar to previous part, we consider driver mutations to occur with the rate
of $P_{m}$ and each driver mutation diffusion constant by the factor of
$\Delta D$ as $D_{i+1}=D_{i}(1+\Delta D)$. As FIG. 3 (d) shows, even
considerably huge rates for driver mutations ($P_{m}=3.4\times 10^{-2}$) and
high fitness advantages ($\Delta D=0.04$) increases growth rate slightly. More
importantly, such an increase in growth rate does not affect growth rate and
as FIG. 3 (e) shows the growth rate remains almost the same.

Emerging new cell types with higher duplication (migration) rate increase
invasion velocity and accelerates the growth, but it does not necessarily lead
to $\beta>1$. To have large values of $\beta$, the newly emerged types should
take over the population in short time intervals. New cell types have small
fitness advantages to overcome existing populations in short times.
Additionally, for compactly growing populations, the clonal sweep (domination
of one type) is a slow process Hallatschek and Fisher (2014). That is why we
do not observe the fast domination of new types that hypothetically could lead
to super-linear growth. Similar results on the negligible effect of driver
mutations on growth rate have been reported previously with an entirely
different model Waclaw et al. (2015). Due to such small effects, it is
difficult to detect driver mutations Bozic et al. (2019).

Figure 3: (a) Emergence of cell types with higher proliferation rate can lead
to accelerated growth which is faster than exponential growth and exhibit
$\beta>1$ Durrett et al. (2010); Pérez-García et al. (2020). (b) Tumor volume
($C_{T}$) versus time for different values of $P_{m}$ and $\Delta R=0.04$. (c)
Growth rate versus tumor volume for the same tumors. (d) Tumor volume versus
time for $\Delta D=0.04$ and different values of $P_{m}$. (e) Growth rate
versus tumor size for the same tumors. Analysis of growth dynamics for linear
evolution suggests that emergence of new sub-populations with higher fitness
advantage is not able to impose a drastic change on $\beta$.

####  III.1.2 Punctuated Evolution

Depending on the value of $\Delta R$, fitness advantage varies across cell
types with a different number of mutations. With $\Delta R=0.004$, $\Delta
R=0.01$ and $\Delta R=0.02$, we have $\frac{R_{100}}{R_{1}}=1.48$,
$\frac{R_{100}}{R_{1}}=2.68$ and $\frac{R_{100}}{R_{1}}=7.24$, respectively.
Such differences are huge compared to experimental results (10-20 percent
difference Heide et al. (2018)), but we include them to check for extreme
cases. We start simulations with uniform distribution across cell types ( by
$C_{i}(0)=C_{0}$) and allow cell populations to evolve based on the
aforementioned rules and Eq. 4 with $P_{m}=0$. As FIG. 4 (a) shows, increasing
$\Delta R$ leads to considerably faster growth. Doubling times for small tumor
sizes ($C_{T}<0.1$ mm3 where the growth is exponential) with $\Delta R=0$,
$\Delta R=0.004$, $\Delta R=0.01$ and $\Delta R=0.02$ are of 21, 17.5, 11 and
6 days, respectively. For small tumor sizes where the growth is exponential,
growth can exceed exponential growth as more fit types start to dominate the
tumor, as suggested by Durrett et al. (2010) and FIG. 3 (a) shows. This
process can lead to $\beta>1$ in small tumor sizes, but it does not lead to
super-linear growth in larger, clinically relevant, tumor sizes as FIG. 4 (b)
shows. Two aspects of this process needs more clarifications: First, for
experimentally relevant parameters ($\Delta R\leq 0.01$) even small tumors do
not show super-linear growth ( $\Delta R=0.004$ at FIG. 4 (b)). Super linear
growth emerges for $\Delta R\sim 0.02$ where the $C_{1}00$ can duplicate 7.25
times faster than $C_{1}$. Clearly such a huge difference in tumor cells
proliferation rate has not been reported yet and seems impossible considering
limitations that cells face in their way to significantly increase their
fitness advantage Hausser and Alon (2020). More importantly, even such a huge
difference can only lead to super-linear growth in small tumor sizes. As FIG.
4 (b) shows, for larger and clinically relevant sizes ($C_{T}\sim 10^{3}$ mm3)
where super-linear growth has been observed Pérez-García et al. (2020), even
these huge differences do not lead to super linear growth.

Figure 4:  (a) Volume versus time for Punctuated evolution with uniform
distribution for cellular densities ($C_{i}(0)=C_{0}$). (b) Growth rate versus
tumor size for the same tumors. Experimentally relevant values of differences
in cellular fitness does not affect growth dynamics. Large values of $\Delta
R$ lead to huge differences in duplication rates and leads to $\beta>1$ for
very small tumor sizes where the growth is exponential for $\Delta R=0$. Even
such a huge differences in duplication rates do not change $\beta$ in larger
and clinically relevant sizes.

For non-uniform density of cell types at initial point ($C_{i}(0)=C_{0}/i$),
we again check for growth dynamics. Similar to the previous part, larger
values of $\Delta R$ lead to faster growth (see FIG. 5 (a)). Interestingly,
non-uniform distribution leads to larger values of $\beta$ compared to uniform
distribution (see FIG. 5 (b)). However, even with such high values for $\beta$
at initial times, $\beta$ approaches 2/3 as the tumors grow larger. As FIG. 5
(c) shows, for $\Delta R>0$, the composition of tumor changes and the
frequency of sub-populations with higher duplication rate increases. Depending
on the value of $\Delta R$, two main scenarios can happen: for larger values
of $\Delta R$, population(s) with higher duplication rate take over the whole
population at initial times and once the tumor becomes slightly larger, the
growth approaches to growth with $\beta=2/3$. On the other hand, if $\Delta$
is smaller, populations with higher duplication rates can not take over the
whole tumor in short time intervals and if finally, they do, the process is
slow and does not affect $\beta$. To cross-check our results, we developed an
individual-based model for exponential growth with the same mutation rates and
fitness advantage variation and the results confirm our findings here author
(supplemtal year).

Figure 5:  (a) Volume versus time for different values of $\Delta R$ and with
$C_{i}(0)=C_{0}/i$ as the initial condition. (b) Growth rate versus tumor size
for the same tumors. While in smaller sizes we can observe $\beta>1$, larger
tumors exhibit $\beta=2/3$. (c) Fraction of each cell type for tumors that
have been growing for 90 days. In tumors with high $\Delta R$, cells with
higher duplication rate are able to take over whole population which leads to
$\beta>1$. later on, such tumor exhibit $\beta=2/3$.

###  III.2 Angiogenesis

As mentioned, tumors normally face an oxygen shortage in the early stages, and
considering a low oxygen supply rate, for example, $b_{n0}=5\times 10^{-6}$
mol s-1mm-3, is rational. With no angiogenesis, such a tumor can grow up to
about 2 $mm^{3}$. Imagine a tumor that is well-vascularized in early times.
Such a tumor is compatible with high values for $b_{n1}$ in our model and
similar to the standard FKPP with $\beta=1$ at early times followed up by
$\beta=2/3$ in longer periods. It should be noted that experiencing hypoxia
triggers angiogenesis and tumors barely develop angiogenesis and this result
has been included for comparison.

Tumors with smaller values for $b_{n1}$, grow until they become close to the
limit imposed by oxygen supply, where most cells become inactive due to oxygen
shortage. Such a halt in growth in spit of initial angiogenesis is known as
post-vascular dormancy. Through a positive feedback loop, the ability to
induce angiogenesis increases oxygen supply and grows larger gradually.
Through this process, even previously hypoxic cells start to duplicate,
resulting in a high acceleration in tumor growth rate and $\beta>1$ (see FIG.
4). In other words, a hypoxic tumor experiences a high acceleration in growth
rate during angiogenesis and exhibits $\beta>1$. The growth rate for such a
tumor can be well below the growth rate of a well-fed tumor of the same size.
Once a tumor became well-oxygenated, in long times, it approaches a standard
FKPP growth with $\beta=2/3$.

Figure 6: Growth rate versus volume for different values of $b_{n1}$ (in mol
cell-1 s-1mm-3). A well-fed tumor initially grows exponentially and then
exhibits a surface growth with $\beta=2/3$. However, hypoxic tumors exhibit
accelerated growth and $\beta>1$ during angiogenesis. The red area shows where
previous models suggested having $\beta>1$ due to growing faster than
exponential growth Durrett et al. (2010); Pérez-García et al. (2020). The blue
area is where our model reveals to exhibit super-linear growth with $\beta>1$.
In this area, growth is not faster than exponential growth, but the
acceleration in growth due to increased oxygen supply leads to $\beta>1$.
Inset: Volume versus time for different values of $b_{n1}$.

##  IV Discussion

Healthy tissues sustain the environment’s properties, such as cellular
duplication and oxygen density, by incorporating a diverse range of regulatory
mechanisms. Once a group of cells avoid growth regulations and duplicate
spontaneously, the rest of the activities (including oxygen regulating
mechanisms) remain normal. As a result, the growing population consumes
existing oxygen and soon faces oxygen shortage and that is why most tumors
experience hypoxia at early stages. To grow larger, tumors need to manipulate
the environment to provide more oxygen, among other alterations required for
survival and progression Hanahan and Weinberg (2000, 2011). Tumor cells
secrete chemicals and growth factors to manipulate the environment, which
leads to a mutualism among tumor cells Archetti and Pienta (2019). While this
process has been known for years, it has not been considered as the dominant
and determinant factor in tumor growthTurajlic et al. (2019), partially due to
ample evidence on tumor cells diversity and the central role that has been
assigned to evolutionary dynamics of competition among tumor cells Parker et
al. (2020); Zahir et al. (2020); Lakatos et al. (2020); Watson et al. (2020).

For a tumor that experiences hypoxia, major parts (if not all) of tumor cells,
even those living on the surface, remain inactive due to oxygen shortage
Rodriguez-Brenes et al. (2013). Respectively, the growth rate for tumors is
only a fraction of its capacity, staying in dormancy Hahnfeldt et al. (1999);
Naumov et al. (2008). If the oxygen supply rate increase due to further
angiogenesis, more and more cells become active, unleashing tumor capability
to grow. Previously, it was suggested that the growth at such a phase might be
exponential Naumov et al. (2008). However, such a transition in the active
fraction of tumor cells can drastically accelerate growth and increase
$\beta>1$ (see FIG. 6). Thus, instead of an increase in individual duplication
rates, which, as our model suggests, is too small to have a relevant effect,
the transition in the number of actively duplicating cells accelerates tumor
growth and leads to $\beta>1$. Respectively, the previously known Hallmark of
cancer —that is common among different cancer types— drives super-linear
dynamics among a wide range of human tumors. The idea of accelerated growth
has been under investigation in ecological literature for a long time Phillips
et al. (2010); Benichou et al. (2012); Hallatschek and Fisher (2014). Once a
sub-population gains driver mutation, it may become dominant. However, an
accelerated growth does not lead to super-linear growth necessarily, as we
showed here.

Inducing angiogenesis through mutually secreting growth factors resembles the
Allee effect Gatenby et al. (2019, 2020), because it leads to an increase in
individual duplication rate by providing more oxygen for tumor cells Korolev
et al. (2014). Our results suggest that observing the super-linear growth,
which is a direct result of mutualism, across human tumors is evidence for the
determinant role of the Allee effect in tumor progression. Additionally, for
populations that exhibit the Allee effect, one can write
$\frac{dC}{dt}=D\nabla^{2}C+RC(1+C)(1-C)$ Fadai et al. (2020). An additional
term as $1+C$ depicts the positive dependency of duplication rate (individual
fitness) on population local density. Depending on the nature of cooperation
among individuals, it is possible to substitute the local density by whole
population size through substituting $C$ with $C_{T}$ ($C_{T}=\int C_{t}dv$).
In the light of this understanding, we can interpret previously presented
model, $\frac{\partial C}{\partial t}=D\nabla^{2}C+C(R_{0}+R_{1}C_{T})(1-C)$
Pérez-García et al. (2020), as an Allee effect. Thus, while the proposed model
was suggested to represent an increase in fitness advantage due to driver
mutations, it is compatible with the well-known Allee effect.

FIG. 6 revealed that super-linear growth emerges in a different parameter zone
and is not necessarily faster than exponential growth. Instead, it is a result
of acceleration. Such an understanding may challenge previous results in
different aspects. One of the major results presented in previous work is that
super-linear growth is an inherently 3D phenomenon Pérez-García et al. (2020).
We checked if 2D environments support super-linear growth. We used a 2D
version of our model to see if the acceleration in growth due to angiogenesis
leads to super-linear growth. Interestingly, even in the 2D model,
angiogenesis leads to super-linear growth author (supplemtal year). As
complementary models, we used a 2D FKPP model Azimzade et al. (2019) and a 2D
individual-based model Azimzade and Saberi (2019) to study growth dynamics in
the presence of an increasing duplication rate. All models show that super-
linear growth emerges in 2D environments as well author (supplemtal year). As
such, it seems that mutualism (and the Allee effect) in 2D environments can
lead to super-linear dynamics, suggesting that it is not an inherently 3D
phenomenon.

###  IV.1 Therapeutic Applications

Understanding the mechanism behind super-linear growth is not only of
theoretical interest. Analysis of growth rate data suggested that tumors with
higher growth rate worth outcomes than those with lower rates Pérez-García et
al. (2020). Previously, the difference between these two types (high-risk and
low risk) was suggested to be initiated from differences in their evolutionary
stage. Since tumors’ evolutionary status is hard (if not impossible) to
interfere with, there were no therapeutic implications. Our model suggests
that angiogenesis is responsible for changes in growth rate and super-linear
dynamics. Respectively, those tumors with higher activity are, in fact, more
successful in inducing angiogenesis. Anti-angiogenetic drugs are well
developed and have been used for different cases and our findings suggest that
such drugs can be effective for high-risk tumors. It should be noted that
anti-angiogenetic drugs may fail due to genotypic changes on endothelial
cells, the vascular normalization phenomenon and the vasculogenic mimicry.
More importantly, anti-angiogenic molecules’ usage leads to a hypoxic tumor
microenvironment, enhancing tumor invasiveness Ribatti et al. (2019). Due to
such complexities, a practical intervention needs to incorporate additional
interactions to predict possible failure.

##  V Conclusion

Recent data revealed super-linear growth across human tumors and acceleration
in growth due to the evolutionary dynamics of competition was introduced as
the driving mechanism behind it. Our model suggests that while the
acceleration of growth leads to super-linear growth, the emergence of cell
types with higher fitness advantages does not lead to such an acceleration.
Instead, hypoxic tumor cells accelerate their growth by manipulating the
environment and inducing angiogenesis. This mutualism in manipulation is
compatible with the Allee effect suggested being involved in tumor growth
previously.

##  VI Acknowledgment

We would like to thank Víctor M. Pérez García for reading our draft and his
fruitful comments.

## References

  * Korolev et al. (2014) K. S. Korolev, J. B. Xavier, and J. Gore, Nature Reviews Cancer 14, 371 (2014).
  * Turajlic et al. (2019) S. Turajlic, A. Sottoriva, T. Graham, and C. Swanton, Nature Reviews Genetics 20, 404 (2019).
  * Bozic et al. (2013) I. Bozic, J. G. Reiter, B. Allen, T. Antal, K. Chatterjee, P. Shah, Y. S. Moon, A. Yaqubie, N. Kelly, D. T. Le, et al., Elife 2, e00747 (2013).
  * Gatenby and Brown (2020) R. A. Gatenby and J. S. Brown, Nature Reviews Clinical Oncology pp. 1–12 (2020).
  * Wang et al. (2014) J. Wang, H. Khiabanian, D. Rossi, G. Fabbri, V. Gattei, F. Forconi, L. Laurenti, R. Marasca, G. Del Poeta, R. Foà, et al., Elife 3, e02869 (2014).
  * Parker et al. (2020) T. Parker, E. Madan, K. Gupta, E. Moreno, and R. Gogna, Trends in Cancer (2020).
  * Zahir et al. (2020) N. Zahir, R. Sun, D. Gallahan, R. A. Gatenby, and C. Curtis, Nature Genetics 52, 759 (2020).
  * Lakatos et al. (2020) E. Lakatos, M. J. Williams, R. O. Schenck, W. C. Cross, J. Househam, L. Zapata, B. Werner, C. Gatenbee, M. Robertson-Tessi, C. P. Barnes, et al., Nature Genetics 52, 1057 (2020).
  * Williams et al. (2020) M. J. Williams, L. Zapata, B. Werner, C. P. Barnes, A. Sottoriva, and T. A. Graham, Elife 9, e48714 (2020).
  * Carmeliet and Jain (2000) P. Carmeliet and R. K. Jain, Nature 407, 249 (2000).
  * Hanahan and Weinberg (2000) D. Hanahan and R. A. Weinberg, cell 100, 57 (2000).
  * Hanahan and Weinberg (2011) D. Hanahan and R. A. Weinberg, Cell 144, 646 (2011).
  * Archetti and Pienta (2019) M. Archetti and K. J. Pienta, Nature Reviews Cancer 19, 110 (2019).
  * Allee and Bowen (1932) W. Allee and E. S. Bowen, Journal of Experimental Zoology 61, 185 (1932).
  * Böttger et al. (2015) K. Böttger, H. Hatzikirou, A. Voss-Böhme, E. A. Cavalcanti-Adam, M. A. Herrero, and A. Deutsch, PLoS Comput Biol 11, e1004366 (2015).
  * Brown et al. (2016) J. S. Brown, J. J. Cunningham, and R. A. Gatenby, IEEE Transactions on Biomedical Engineering 64, 512 (2016).
  * Gatenby et al. (2019) R. A. Gatenby, J. Zhang, and J. S. Brown, Cancer Research 79, 3174 (2019).
  * Johnson et al. (2019) K. E. Johnson, G. Howard, W. Mo, M. K. Strasser, E. A. Lima, S. Huang, and A. Brock, PLoS Biology 17, e3000399 (2019).
  * Rodriguez-Brenes et al. (2013) I. A. Rodriguez-Brenes, N. L. Komarova, and D. Wodarz, Trends in Ecology & Evolution 28, 597 (2013).
  * West and Newton (2019) J. West and P. K. Newton, Proceedings of the National Academy of Sciences 116, 1918 (2019).
  * Collins (1956) V. P. Collins, Am J Roentgenol 76, 988 (1956).
  * Mayneord (1932) W. V. Mayneord, The American Journal of Cancer 16, 841 (1932).
  * West et al. (2001) G. B. West, J. H. Brown, and B. J. Enquist, Nature 413, 628 (2001).
  * Pérez-García et al. (2020) V. M. Pérez-García, G. F. Calvo, J. J. Bosque, O. León-Triana, J. Jiménez, J. Pérez-Beteta, J. Belmonte-Beitia, M. Valiente, L. Zhu, P. García-Gómez, et al., Nature Physics pp. 1–6 (2020).
  * Durrett et al. (2010) R. Durrett, J. Foo, K. Leder, J. Mayberry, and F. Michor, Theoretical Population Biology 78, 54 (2010).
  * Phillips et al. (2010) B. L. Phillips, G. P. Brown, and R. Shine, Journal of evolutionary biology 23, 2595 (2010).
  * Benichou et al. (2012) O. Benichou, V. Calvez, N. Meunier, and R. Voituriez, Physical Review E 86, 041908 (2012).
  * Hallatschek and Fisher (2014) O. Hallatschek and D. S. Fisher, Proceedings of the National Academy of Sciences 111, E4911 (2014).
  * Bozic et al. (2019) I. Bozic, C. Paterson, and B. Waclaw, PLoS computational biology 15, e1007368 (2019).
  * Byrne (2010) H. M. Byrne, Nature Reviews Cancer 10, 221 (2010).
  * Altrock et al. (2015) P. M. Altrock, L. L. Liu, and F. Michor, Nature Reviews Cancer 15, 730 (2015).
  * Schrek (1935) R. Schrek, The American Journal of Cancer 24, 807 (1935).
  * Anderson and Quaranta (2008) A. R. Anderson and V. Quaranta, Nature Reviews Cancer 8, 227 (2008).
  * Metzcar et al. (2019) J. Metzcar, Y. Wang, R. Heiland, and P. Macklin, JCO clinical cancer informatics 2, 1 (2019).
  * Azimzade and Saberi (2019) Y. Azimzade and A. A. Saberi, Journal of Statistical Mechanics: Theory and Experiment 2019, 103502 (2019).
  * Swanson et al. (2003) K. R. Swanson, C. Bridge, J. Murray, and E. C. Alvord Jr, Journal of the neurological sciences 216, 1 (2003).
  * Mandonnet et al. (2003) E. Mandonnet, J.-Y. Delattre, M.-L. Tanguy, K. R. Swanson, A. F. Carpentier, H. Duffau, P. Cornu, R. Van Effenterre, E. C. Alvord Jr, and L. Capelle, Annals of Neurology 53, 524 (2003).
  * Azimzade et al. (2019) Y. Azimzade, A. A. Saberi, and M. Sahimi, Physical Review E 100, 062409 (2019).
  * Rejniak and Anderson (2011) K. A. Rejniak and A. R. Anderson, Wiley Interdisciplinary Reviews: Systems Biology and Medicine 3, 115 (2011).
  * Jiménez-Sánchez et al. (2021) J. Jiménez-Sánchez, Á. Martínez-Rubio, A. Popov, J. Pérez-Beteta, Y. Azimzade, D. Molina-García, J. Belmonte-Beitia, G. F. Calvo, and V. M. Pérez-García, PLOS Computational Biology 17, e1008266 (2021).
  * Azimzade et al. (2020) Y. Azimzade, M. Sasar, and I. Maleki, Scientific Reports 10, 1 (2020).
  * Johnston et al. (2017) S. T. Johnston, R. E. Baker, D. S. McElwain, and M. J. Simpson, Scientific reports 7, 1 (2017).
  * Buchakjian and Kornbluth (2010) M. R. Buchakjian and S. Kornbluth, Nature Reviews Molecular Cell Biology 11, 715 (2010).
  * Scalerandi and Sansone (2002) M. Scalerandi and B. C. Sansone, Physical Review Letters 89, 218101 (2002).
  * Azimzade et al. (2018) Y. Azimzade, A. A. Saberi, and M. Sahimi, Scientific Reports 8, 1 (2018).
  * Bozic et al. (2010) I. Bozic, T. Antal, H. Ohtsuki, H. Carter, D. Kim, S. Chen, R. Karchin, K. W. Kinzler, B. Vogelstein, and M. A. Nowak, Proceedings of the National Academy of Sciences 107, 18545 (2010).
  * Anderson (2005) A. R. Anderson, Mathematical Medicine and Biology 22, 163 (2005).
  * Murray (2003) J. D. Murray, _Mathematical Biology I: An Introduction_ (Springer-Verlag, 2003).
  * Harpold et al. (2007) H. L. Harpold, E. C. Alvord Jr, and K. R. Swanson, Journal of Neuropathology & Experimental Neurology 66, 1 (2007).
  * Whiteside (2008) T. Whiteside, Oncogene 27, 5904 (2008).
  * Davis et al. (2017) A. Davis, R. Gao, and N. Navin, Biochimica et Biophysica Acta (BBA)-Reviews on Cancer 1867, 151 (2017).
  * Gould and Eldredge (1972) S. J. Gould and N. Eldredge, Models in paleobiology pp. 82–115 (1972).
  * Nesta et al. (2020) A. V. Nesta, D. Tafur, and C. R. Beck, Trends in Genetics (2020).
  * Vermeulen et al. (2013) L. Vermeulen, E. Morrissey, M. Van Der Heijden, A. M. Nicholson, A. Sottoriva, S. Buczacki, R. Kemp, S. Tavaré, and D. J. Winton, Science 342, 995 (2013).
  * McFarland et al. (2014) C. D. McFarland, L. A. Mirny, and K. S. Korolev, Proceedings of the National Academy of Sciences 111, 15138 (2014).
  * Castro-Giner et al. (2015) F. Castro-Giner, P. Ratcliffe, and I. Tomlinson, Nature Reviews Cancer 15, 680 (2015).
  * Persi et al. (2020) E. Persi, Y. I. Wolf, D. Horn, E. Ruppin, F. Demichelis, R. A. Gatenby, R. J. Gillies, and E. V. Koonin, Nature Reviews Genetics pp. 1–12 (2020).
  * Hausser and Alon (2020) J. Hausser and U. Alon, Nature Reviews Cancer 20, 247 (2020).
  * Hatzikirou et al. (2012) H. Hatzikirou, D. Basanta, M. Simon, K. Schaller, and A. Deutsch, Mathematical Medicine and Biology 29, 49 (2012).
  * Waclaw et al. (2015) B. Waclaw, I. Bozic, M. E. Pittman, R. H. Hruban, B. Vogelstein, and M. A. Nowak, Nature 525, 261 (2015).
  * Paterson et al. (2020) C. Paterson, H. Clevers, and I. Bozic, Proceedings of the National Academy of Sciences 117, 20681 (2020).
  * Heide et al. (2018) T. Heide, L. Zapata, M. J. Williams, B. Werner, G. Caravagna, C. P. Barnes, T. A. Graham, and A. Sottoriva, Nature Genetics 50, 1633 (2018).
  * Williams et al. (2016) M. J. Williams, B. Werner, C. P. Barnes, T. A. Graham, and A. Sottoriva, Nature genetics 48, 238 (2016).
  * Phillips et al. (2020) C. M. Phillips, E. A. Lima, R. T. Woodall, A. Brock, and T. E. Yankeelov, Plos One 15, e0231137 (2020).
  * Zheng and Sweidan (2018) X. Zheng and M. Sweidan, Journal of Mathematical Biology 77, 1589 (2018).
  * d’Onofrio (2007) A. d’Onofrio, Physical Review E 76, 031920 (2007).
  * Rieger and Welter (2015) H. Rieger and M. Welter, Wiley Interdisciplinary Reviews: Systems Biology and Medicine 7, 113 (2015).
  * Tracqui (2009) P. Tracqui, Reports on Progress in Physics 72, 056701 (2009).
  * Xu et al. (2016) J. Xu, G. Vilanova, and H. Gomez, PloS One 11, e0149422 (2016).
  * Hahnfeldt et al. (1999) P. Hahnfeldt, D. Panigrahy, J. Folkman, and L. Hlatky, Cancer Research 59, 4770 (1999).
  * Naumov et al. (2008) G. N. Naumov, J. Folkman, O. Straume, and L. A. Akslen, Apmis 116, 569 (2008).
  * Haass et al. (2014) N. K. Haass, K. A. Beaumont, D. S. Hill, A. Anfosso, P. Mrass, M. A. Munoz, I. Kinjyo, and W. Weninger, Pigment Cell & Melanoma Research 27, 764 (2014).
  * Bray (2001) D. Bray, _Cell Movements: from molecules to motility_ (Garland Science, 2001).
  * author (supplemtal year) See Supplemental Material at [url] for complementrary indivdual based model that confirm our results and different 2D models which exhibit super-linear growth dynamics.
  * Watson et al. (2020) C. J. Watson, A. Papula, G. Y. Poon, W. H. Wong, A. L. Young, T. E. Druley, D. S. Fisher, and J. R. Blundell, Science 367, 1449 (2020).
  * Gatenby et al. (2020) R. A. Gatenby, Y. Artzy-Randrup, T. Epstein, D. R. Reed, and J. S. Brown, Cancer Research 80, 613 (2020).
  * Fadai et al. (2020) N. T. Fadai, S. T. Johnston, and M. J. Simpson, bioRxiv p. 774000 (2020).
  * Ribatti et al. (2019) D. Ribatti, T. Annese, S. Ruggieri, R. Tamma, and E. Crivellato, Translational Oncology 12, 981 (2019).

##  VII  Supplementary Material

###  VII.1  Agent-based exponential growth

We present an agent-based model for exponential growth. In this model, we
consider two scenarios as well: linear (Darwinian evolution) and punctuated
evolution.

We first have a basic model in which all cells belong to the same type with no
mutations. We randomly select cells for duplication with a probability of $R$.
At each time step, the total number of selections is equal to the cell number.
Since there is no restriction for duplication, the whole number of cells, $N$,
exponentially increases over time.

For Darwinian evolution, we start with 100 cells that belong to the same type.
We randomly select cells for duplication. The selected cell will duplicate
with a probability of $R_{1}$. During each duplication, one of the daughter
cells belongs to the same type and the next cell belongs to the next
generation with a probability of $P_{m}$. We run simulations until the total
number of cells reach $10^{10}$, a typical cell number in a tumor. Duplication
rate for each generation increases as: $R_{i+1}=R_{i}(1+\Delta R)$ with
$\Delta R=0.004$ Bozic et al. (2010).

For punctuated evolution we consider 100 cells each belonging to a different
generations and while the mutation rate is zero, duplication rate increase for
different generations as: $R_{i+1}=R_{i}(1+\Delta R)$ with $\Delta R=0.004$.
We run simulations until the total number of cells increases up to $10^{10}$,
a typical cell number in a tumor with clinically relevant size. We check for
the number of cells and the number of new-born cells in each time step. As FIG
7 (a) shows, cell number exponentially grows for all cases. Number of new-born
cells linearly increases versus cell number (see FIG 7 (b)). As such, even in
exponentially growing populations, evolutionary dynamics do not lead to super-
linear growth.

Figure 7: Exponential Growth.

###  VII.2 2D models

To check if 2D environments also exhibit super-linear growth, we use three
different models. First, we develop a 2D version of our model. To this end, we
consider the simulation box to be 2D. As the second model, we use a FKPP
equation with increasing reaction term as below:

$\displaystyle\frac{\partial C}{\partial
t}=D_{c}\nabla^{2}C+(R_{0}+R_{1}C_{T})C(1-C)$ (7)

with $C_{T}=\int C(x,y)dxdy$.

Finally we developed an agent-based model Azimzade and Saberi (2019) in which
instead of a fixed duplication rate of $R_{0}$, as tumors becomes larger,
duplication rate increases as $R=R_{0}+R_{1}N$ with $N=\Sigma C$. As the FIGs
8-10 show, all these models can lead to emergence of super-linear growth.

Figure 8: 2D version of our model. Tumor volume versus time (a) and growth
activity versus volume (b).

Figure 9: 2D FKPP.

Figure 10: 2D Discrete.

##  VIII Data Availability

Codes for different models are provided at:
https://gitlab.com/YounessAzimzade/superlinearallee


# A Note on Effects of Generalized and Extended Uncertainty Principles on
Jüttner Gas

Hooman<EMAIL_ADDRESS>Sara
<EMAIL_ADDRESS>Amir Hadi
<EMAIL_ADDRESS>1 Research Institute for Astronomy and
Astrophysics of Maragha (RIAAM), University of Maragheh, Maragheh P.O. Box
55136-553, Iran
2 Department of Physics, Faculty of Sciences, Yasouj University, Yasouj
75918-74934 Iran

###### Abstract

In recent years, the implications of the generalized (GUP) and extended (EUP)
uncertainty prin-ciples on Maxwell-Boltzmann distribution have been widely
investigated. However, at high energy regimes, the validity of Maxwell-
Boltzmann statistics is under debate and instead, the Jüttner distribution is
proposed as the distribution function in relativistic limit. Motivated by
these considerations, in the present work, our aim is to study the effects of
GUP and EUP on a system that obeys the Jüttner distribution. To achieve this
goal, we address a method to get the distribution function by starting from
the partition function and its relation with thermal energy which finally
helps us in finding the corresponding energy density states.

##  I Introduction

A general prediction of any quantum gravity theory is the possibility of the
existence of a minimal length in nature, known as the Planck length, below
which no other length can be observed. It is commonly believed that in the
vicinity of the Planck length, the smooth structure of spacetime is replaced
by a foamy structure due to quantum gravity effects minimallength . Therefore,
the Planck scale can be regarded as a separation line between classical and
quantum gravity regimes. There is a general consensus that in the scale of
this minimal size, the characteristics of different physical systems would be
altered. For instance, the introduction of a minimal length scale results in a
generalization of the Heisenberg un-certainty principle (HUP) in such a way
that it incorporates gravitationally induced un-certainty, postulated as the
generalized uncertainty principle (GUP) GUP1993 . In fact, the HUP breaks down
for energies near the Planck scale, i.e., when the Schwarz-schild radius is
comparable to the Compton wavelength and both are close to the Planck length.
This deficiency is removed by revising the characteristic scale through the
modifi-cation of HUP to GUP.

In recent decades, numerous studies on the effects of GUP in various classical
and quantum mechanical systems have been performed (see e.g., GUPwork1 ;
GUPwork2 ; GUPwork3 ; GUPwork4 ). Uncertainty in momentum is also bounded from
below and it is proposed that its minimum is non-zero, a proposal which
modifies HUP to the extended uncertainty principle (EUP) EUPwork . In the
presence of EUP and GUP, the general form of modified HUP is proposed as

$\displaystyle\Delta x\Delta p\geq\frac{\hbar}{2}\left(1+\alpha(\Delta
x)^{2}+\eta(\Delta p)^{2}+\gamma\right),$ (1)

in which $\alpha$, $\eta$ and $\gamma$ are positive deformation parameters
generalform . It should be noted that there is another formulation of GUP and
EUP GUPEUP2 , and also the extended forms of HUP, like GUP, may break the
fundamental symmetries such as Lorentz invariance and CPT symGUP .

On the other hand, it is known that heavy ions can be accelerated to very high
kinetic energies constituting an ensemble of ideal gas with relativistic
velocities in large particle accelerators heavyions . In such high energy
regimes, minimal length effects may appear and could have their own influences
on the statistics of ideal gases. Therefore, particle accelerators could
provide a setting to examine the phenomena related to short-distance physics
gravityac . Based on minimum observable length, the quantum gravity
implications on the statistical properties of ideal gases have been
investigated in many studies, see, e.g., phase1 ; phase2 ; phase3 ;
Hossenfelder2013  and references therein. In the framework of GUP, $i$)
deformed density of states and an improved definition of the statistical
entropy have been introduced in densGUP , $ii$) Maxwell-Boltzmann statistics
is investigated in maxwellGUP , and $iii$) employing Maxwell-Boltzmann
statistics, the thermodynamics of relativistic ideal gas has also analyzed in
thermalGUP . In the same manner, there have been some works on the deformation
of statistical concepts in the EUP framework statEUP .

Jüttner distribution is a generalization of Maxwell-Boltzmann statistics to
the relativistic regimes, which appears in high energy physics. Since quantum
gravity is a high energy physics scenario, its statistical effects may be more
meaningful in the framework of Jüttner distribution function compared to the
Maxwell-Boltzmann distribution jut . Here, our main aim is to study the
effects of GUP and EUP on Maxwell-Boltzmann and Jüttner distributions and
density of states in energy space. To achieve this goal, we begin by providing
an introductory note on Maxwell-Boltzmann and Jüttner distributions. We then
address a way to find these functions by starting from the partition function
of the system in the next section. The effects of GUP and EUP on these
statistics are also studied in the subsequent sections, respectively. The last
section is devoted to a summary on the work.

##  II The Maxwell-Boltzmann and Jüttner distribution functions

We begin by considering an ideal gas composed of non-interacting particles and
set the units so that $\omega_{0}=2\pi\hbar=1$, where $\omega_{0}$ denotes the
fundamental volume of each cell in the two-dimensional phase-space. This value
of $\omega_{0}$ originates from the well-known commutation relation between
the canonical coordinates $x$ and $p$, and indeed, it is the direct result of
HUP pathria . Therefore, any changes in HUP can affect $\omega_{0}$.

### Non-relativistic gas

Let us consider a $3$-dimensional classical gas consisting of $N$ identical
non-interacting particles of mass $m$ with $E=\frac{m}{2}v^{2}$, where $E$ and
$v$ denote the energy and velocity of each particle, respectively. At
temperature $T$, the Maxwell-Boltzmann (MB) distribution function is given by

$\displaystyle f_{{\rm MB}}(v,\beta)=Z_{{\rm MB}}\ \exp\left(-\frac{\beta
mv^{2}}{2}\right),$ (2)

where $Z_{{\rm MB}}$ is a normalization constant, and $\beta\equiv{1}/{K_{\rm
B}T}$ with $K_{\rm B}$ being the Boltzmann constant. In terms of $E$ we have

$\displaystyle f_{{\rm MB}}(E,\beta)$ $\displaystyle=$ $\displaystyle 4\pi
v^{2}(E)\ f_{{\rm MB}}(v(E),\beta)\frac{dv}{dE}$ (3) $\displaystyle=$
$\displaystyle Z_{{\rm MBE}}\ E^{\frac{1}{2}}\exp(-\beta E),$

in which $Z_{{\rm MBE}}$ is a new normalization constant and $E^{\frac{1}{2}}$
denotes the density of states with energy $E$. The normalization constants can
be calculated using the normalization constraint

$\displaystyle\int_{0}^{\infty}f_{{\rm
MB}}(v,\beta)d^{3}v=\int_{0}^{\infty}f_{{\rm MB}}(E,\beta)dE=1$ (4)

The extremum of $f_{{\rm MB}}(E,\beta)$ is located at
$E=\frac{1}{2\beta}\equiv E_{{\rm MB}}^{{\rm ext}}$ or equally
$v=\frac{1}{\sqrt{\beta m}}\equiv v^{{\rm ext}}_{{\rm MB}}$. One can also
evaluate the partition function of the mentioned gas (with Hamiltonian
$H=\frac{p^{2}}{2m}$) as

$\displaystyle Q_{N}=\frac{1}{N!}\underbrace{\int_{...}\int}_{3N}\exp(-\beta
H)\ d^{3N}x\ d^{3N}p=\frac{\left(Q_{1}^{\rm NR}\right)^{N}}{N!},$ (5)

where

$\displaystyle Q^{{\rm NR}}_{1}=\int\exp(-\beta H)\ d^{3}x\
d^{3}p=V\left(\frac{2\pi m}{\beta}\right)^{\frac{3}{2}},$ (6)

denotes the single particle partition function of a non-relativistic gas and
$V$ refers to the total volume of the system. In this manner, the
corresponding thermal energy per particle ($U$) takes the form

$\displaystyle U^{{\rm NR}}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}E\ f_{{\rm MB}}(E,\beta)dE=-\frac{\partial\ln
Q^{{\rm NR}}_{1}}{\partial\beta}=\frac{3}{2\beta}$ (7) $\displaystyle=$
$\displaystyle 3E_{{\rm MB}}^{{\rm ext}}.$

Although the use of Eq. (7) dates back to before the discovery of the special
relativity theory by Einstein, the ultra-relativistic expression of $E$
produces interesting results in this framework pathria .

### Relativistic gas

In the relativistic situations, where $E=\sqrt{p^{2}c^{2}+m^{2}c^{4}}$ in
which $c$ denotes the light velocity and $m$ is the rest mass, one can employ
Eq. (5) to get

$\displaystyle Q_{1}^{\rm R}$ $\displaystyle=$ $\displaystyle Q_{1}^{\rm
NR}\Psi(\sigma),$ $\displaystyle\Psi(\sigma)$ $\displaystyle=$
$\displaystyle\frac{i^{\frac{3}{2}}m^{3}c^{6}\sqrt{\frac{\pi}{2}}\
H_{2}^{(1)}(i\sigma)}{(i\sigma)^{\frac{5}{2}}},$ (8)

as the partition function of single particle pauli ; pauli1 ; jut . Finally,
we obtain the thermal energy per particle as

$\displaystyle U^{{\rm
R}}=\frac{1}{\beta}\left[1-i\sigma\frac{H_{2}^{\prime(1)}(i\sigma)}{H_{2}^{(1)}(i\sigma)}\right]=-\frac{\partial}{\partial\beta}\ln
Q_{1}^{\rm R}.$ (9)

In the above equations, $\sigma=\beta mc^{2}$, $H^{(j)}_{n}(i\sigma)$ is the
$n$-th order Hankel function of the $j$-th kind, and prime denotes derivative
with respect to the argument of function. The above results were first
reported by in $1911$ Jüttner jut , who attempt to calculate the energy of a
relativistic ideal gas using the conventional theory of relativistic
statistical mechanics. According to Jüttner’s results, a comprehensive study
of a $3$-dimensional relativistic system requires the Jüttner distribution
($f_{\rm J}$)

$\displaystyle f_{\rm J}(\gamma,\beta)=Z_{\rm
J}(\gamma^{2}-1)^{\frac{1}{2}}\gamma\exp\big{(}-\beta m\gamma\big{)},$ (10)

instead of MB distribution ($f_{{\rm MB}}$) jut1 ; ange ; prl ; mon ; ghod ;
njp . Jüttner distribution is indeed the relativistic extension of generalized
isotropic MB distribution when $E(p)=m\gamma(p)c^{2}$. Here, $Z_{\rm J}$ is
the normalization constant and $\gamma=\frac{1}{\sqrt{1-v^{2}}}$ refers to the
Lorentz factor, where the units have been set so that $c=1$. In terms of $E$,
simple calculations give us Jüttner distribution as

$\displaystyle f_{\rm J}(E,\beta)=Z_{{\rm JE}}a_{\rm J}(E)\exp\big{(}-\beta
E\big{)},$ (11)

where $a_{\rm J}(E)=E\sqrt{E^{2}-m^{2}}$ denotes the density of states in
energy representation, and in terms of $v$ one finds

$\displaystyle f_{\rm J}(v,\beta)=Z_{{\rm
JV}}\left(\frac{1}{\sqrt{1-v^{2}}}\right)^{5}\exp\left(-\beta\frac{m}{\sqrt{1-v^{2}}}\right),$
(12)

in which $Z_{{\rm JE}}$ and $Z_{{\rm JV}}$ are new normalization constants
ange . These constants can be evaluated using the normalization condition

$\displaystyle\int_{1}^{\infty}f_{\rm J}(\gamma,\beta)\ d\gamma$
$\displaystyle=$ $\displaystyle\int_{m}^{\infty}f_{\rm J}(E,\beta)dE=1$ (13)
$\displaystyle=$ $\displaystyle\int_{0}^{1}f_{\rm
J}(\gamma(v),\beta)\frac{d\gamma}{dv}dv$ $\displaystyle=$
$\displaystyle\int_{0}^{1}f_{\rm J}(v,\beta)\ d^{3}v,$

which clearly states that

$\displaystyle f_{\rm J}(v,\beta)=\frac{1}{4\pi v^{2}}f_{\rm
J}(\gamma(v),\beta)\frac{d\gamma}{dv}.$ (14)

It is finally useful to note that the extremum of $f_{\rm J}(v,\beta)$ is
located at $v=\sqrt{1-(\frac{\beta m}{5})^{2}}\equiv v^{{\rm ext}}_{\rm J}$
leading to $E^{{\rm ext}}_{\rm J}=\frac{5}{\beta}$, a solution which is valid
only when $\beta m<5$. There are also other proposals for Jüttner function
($f_{\rm J}(\gamma,\beta)$) jut1 ; ange ; prl ; mon ; ghod ; njp , but the
standard form Eq. (11) is considered in this paper is confirmed by some
previous studies prl ; mon ; ghod . The corresponding thermal energy per
particle (i.e., $\langle\gamma\rangle$) (or equally, the ratio $\frac{U}{N}$
in Eq. (12)) can also be obtained by using $f_{\rm J}(\gamma,\beta)$, as

$\displaystyle U^{{\rm R}}\equiv m\langle\gamma\rangle=\int_{m}^{\infty}E\
f_{\rm J}(E,\beta)\ dE=-\frac{\partial}{\partial\beta}\ln Q_{1}^{\rm R}.$ (15)

Although Eqs. (7) and (15) are simple examples, they confirm that the mean
value of energy (or equally, thermal energy) can be calculated by using either
the partition function or the distribution function. Moreover, employing these
equations, one can find the distribution functions whenever the partition
function is known. Indeed, if the phase-space geometry is deformed, then the
partition function will also be modified. Therefore, one can find the
corresponding modified MB and Jüttner distributions by directly using Eqs. (7)
and (15) for the non-relativistic cases, repetitively.

##  III Generalized uncertainty principle, partition and distribution
functions

In the units of $\hbar=c=1$, $[x_{k},p_{l}]=i\delta_{kl}$ is the standard
commutation relation between the canonical coordinates $x$ and $p$. This
relation leads to HUP in the framework of quantum mechanics and is the
backbone of calculating $\omega_{0}$ pathria ; maxwellGUP . Thus, the volume
element $d^{3}x\ d^{3}p$ changes whenever different coordinates (commutation
relations) are used phase1 ; phase2 ; phase3 ; maxwellGUP . For GUP, we have
phase2 ; prdgup

$\displaystyle\Delta X\Delta P\geq\frac{1}{2}[1+\eta(\Delta P)^{2}+...],$ (16)

where $\eta$ denotes the GUP parameter and is based on modified commutation
relations

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\big{[}X_{k},P_{l}\big{]}$
$\displaystyle=$ $\displaystyle i\bigg{(}\delta_{kl}(1+\eta
P^{2})+\eta^{\prime}P_{k}P_{l}\bigg{)},$
$\displaystyle\big{[}P_{k},P_{l}\big{]}$ $\displaystyle=$ $\displaystyle 0,$
$\displaystyle\\!\\!\\!\\!\\!\big{[}X_{k},X_{l}\big{]}$ $\displaystyle=$
$\displaystyle i\frac{2\eta-\eta^{\prime}+(2\eta+\eta^{\prime})\eta
P^{2}}{1+\eta P^{2}}$ (17) $\displaystyle(P_{k}X_{l}-P_{l}X_{k}),$

where $k,l=1,2,3$ for a $3$-dimensional space phase2 ; prdgup . $P$ and $X$
are generalized coordinates which are not necessarily equal to the canonical
coordinates $p$ and $x$. In this manner, assuming $\eta^{\prime}=0$ and $\eta$
is independent of $\hbar$, one finds

$\displaystyle d^{3}x\ d^{3}p\rightarrow\frac{d^{3}X\ d^{3}P}{(1+\eta
P^{2})^{3}},$ (18)

which must be considered as the volume element in $X-P$ space instead of
$d^{3}x\ d^{3}p$ phase1 ; phase2 ; phase3 ; maxwellGUP . This means that the
density of states in the $X-P$ phase space is affected by the factor of
$(1+\eta P^{2})$ phase1 . In this situation, the single particle partition
function can also be found as

$\displaystyle Q_{1}^{{\rm GUP}}=\int\exp(-\beta H(P,X))\frac{d^{3}X\
d^{3}P}{(1+\eta P^{2})^{3}},$ (19)

where $H(P,X)$ denotes Hamiltonian in generalized coordinates phase1 ; phase2
; thermalGUP ; statEUP . The corresponding thermal energy ($U^{{\rm GUP}}$)
can be calculated using the relation

$\displaystyle U^{{\rm GUP}}=-\frac{\partial}{\partial\beta}\ln Q_{1}^{{\rm
GUP}},$ (20)

along with Eq. (19), which finally gives

$\displaystyle U^{{\rm GUP}}=\frac{\int H(P)\ \exp(-\beta
H(P))\frac{d^{3}P}{(1+\eta P^{2})^{3}}}{\int\exp(-\beta
H(P))\frac{d^{3}P}{(1+\eta P^{2})^{3}}}.$ (21)

In obtaining this equation, the fact that $H(\equiv E)$ is independent of $X$
has been used which cancels integration over $d^{3}X$. Indeed, the density of
states in phase-space is changed under the shadow of GUP maxwellGUP ; phase2 ;
phase1 , a result which affects the distribution function.

For a single free particle with $H=\frac{P^{2}}{2m}$, the ideal gas law is
still valid, and therefore

$\displaystyle Q_{1}^{{\rm NR,GUP}}=Q^{{\rm NR}}_{1}\
I\left(\frac{2m\eta}{\beta},3\right),$ (22)

while the explicit form of the function
$I\big{(}\frac{2m\eta}{\beta},3\big{)}$ can be followed in phase3  and
$Q_{1}^{{\rm NR}}$ is introduced in Eq. (6). The effects of GUP are stored in
$I\big{(}\frac{2m\eta}{\beta},3\big{)}$, and at the limit of $\eta\rightarrow
0$, one gets the ordinary single partition function of a free particle.

Correspondingly, the partition function of a single free relativistic particle
can also be evaluated using $H^{2}=P^{2}+m^{2}$ in Eq. (19). By doing so one
finds

$\displaystyle Q_{1}^{{\rm
R,GUP}}=\int\exp\left(-\beta\sqrt{P^{2}+m^{2}}\right)\frac{d^{3}X\
d^{3}P}{(1+\eta P^{2})^{3}},$ (23)

for which the solution reads

$\displaystyle Q_{1}^{{\rm R,GUP}}=Q^{{\rm
R}}_{1}\left(1-\eta\frac{15}{2}\frac{1}{\beta m}\right),$ (24)

when $m\gg\frac{1}{\beta}$ thermalGUP ; statEUP .

### Maxwell-Boltzmann statistics

Bearing the recipe which led to the expression for $f_{{\rm MB}}(E,\beta)$ in
mind, one can get the modified MB distribution in the $X-P$ space as

$\displaystyle f_{{\rm MB}}^{{\rm GUP}}(E,\beta)$ $\displaystyle=$
$\displaystyle 4\pi P^{{}^{2}}(E)\frac{\exp(-\beta E)}{(1+\eta
P^{{}^{2}}(E))^{3}}\frac{dP}{dE}$ (25) $\displaystyle=$ $\displaystyle Z_{{\rm
MBE}}^{{\rm GUP}}\ \frac{E^{\frac{1}{2}}\exp(-\beta E)}{(1+2\eta mE)^{3}},$

where $Z_{{\rm MBE}}^{{\rm GUP}}$ denotes the normalization constant in the
presence of GUP. The thermal energy then reads

$\displaystyle U^{{\rm GUP}}=\int_{0}^{\infty}E\ f_{{\rm MB}}^{{\rm
GUP}}(E,\beta)\ dE.$ (26)

One can also find the normalization constant $Z_{{\rm MBE}}^{{\rm GUP}}$ as

$\displaystyle Z_{{\rm MBE}}^{{\rm GUP}}$ $\displaystyle=$
$\displaystyle\bigg{(}\int_{0}^{\infty}\frac{E^{\frac{1}{2}}\exp(-\beta
E)}{(1+2\eta mE)^{3}}dE\bigg{)}^{-1},$ (27)

which is equal to $Z_{{\rm MBE}}$ in the limit where $\eta\rightarrow 0$.
Obviously, the MB distribution $f_{MB}(E,\beta)$ is recovered through Eq. (25)
at the appropriate limit of $\eta=0$. For the density of states in HUP
framework we have $a_{{\rm MB}}(E)=\sqrt{E}$. This relation is modified in the
presence of GUP effects and thus, the density of states will take the
following form

$\displaystyle a_{{\rm MB}}^{{\rm GUP}}(E)=\frac{\sqrt{E}}{(1+2\eta mE)^{3}},$
(28)

which is in agreement with the results of maxwellGUP . The extremum of
$f_{{\rm MB}}^{{\rm GUP}}(E,\beta)$ is also located at

$\displaystyle\mathcal{E}_{{\rm MB}}^{{\rm ext}}=$
$\displaystyle\frac{1+10m\eta E_{{\rm MB}}^{{\rm
ext}}}{4m\eta}\left(\sqrt{1+\frac{8m\eta E_{{\rm MB}}^{{\rm ext}}}{(1+10m\eta
E_{MB}^{ex})^{2}}}-1\right),$

which clearly indicates $\mathcal{E}_{{\rm MB}}^{{\rm ext}}\rightarrow E_{{\rm
MB}}^{{\rm ext}}$ whenever $\eta\rightarrow 0$. In Fig. (1), the effects of
GUP on the distribution function in MB statistics are shown where the
temperature is considered to be constant ($\beta=1$).

Figure 1: MB distribution versus energy for $\eta=0.5,1,1.5$. The ordinary MB
distribution is denoted by the solid curve. Here, we have set the units so
that $\hbar=c=k_{B}=1$.

### Jüttner statistics

In the relativistic situation, where $H=\sqrt{P^{2}+m^{2}}(\equiv E)$,
following the above recipe, we get the modified Jüttner distribution as

$\displaystyle f_{\rm J}^{{\rm GUP}}(E,\beta)=Z_{{\rm JE}}^{{\rm GUP}}\
\frac{E\sqrt{E^{2}-m^{2}}\exp\big{(}-\beta
E\big{)}}{\big{(}1+\eta[E^{2}-m^{2}]\big{)}^{3}},$ (30)

which recovers $f_{\rm J}(E,\beta)$ in the limit where $\eta\rightarrow 0$.
Here, $Z_{{\rm JE}}^{{\rm GUP}}$ is also the normalization constant evaluated
by utilizing the normalization constraint $\int_{m}^{\infty}f_{J}^{{\rm
GUP}}(E,\beta)\ dE=1$. We also find

$\displaystyle a^{{\rm GUP}}_{\rm
J}(E)=\frac{E\sqrt{E^{2}-m^{2}}}{\big{(}1+\eta[E^{2}-m^{2}]\big{)}^{3}},$ (31)

as the density of states in Jüttner statistics in the presence of GUP. Figure
(2) shows the behavior of $f_{\rm J}^{{\rm GUP}}(E,\beta)$ for some positive
values of $\eta$ parameter.

Figure 2: Behavior of Jüttner distribution against energy for
$\eta=0.5,1,1.5$. The solid curve belongs to the ordinary Jüttner distribution
and we have set the units so that $\hbar=c=k_{B}=1$.

##  IV Extended uncertainty principle, partition and distribution functions

The modified Heisenberg algebra in the EUP framework can be recast into the
following form

$\displaystyle[X_{i},P_{j}]=i(\delta_{ij}+\alpha X_{i}X_{j}),$ (32)

where $\alpha$ is a small positive parameter known as the EUP parameter. In
the limit of $\alpha\rightarrow 0$, the canonical commutation relation of the
standard quantum mechanics is recovered. Based on the commutation relation
(32), the HUP is modified by

$\displaystyle(\Delta X_{i})(\Delta P_{i})\geq\frac{1}{2}[1+\alpha(\Delta
X_{i})^{2}+...],$ (33)

which leads to a non-zero minimum uncertainty in momentum as $(\Delta
P_{i})_{\rm min}=\sqrt{\alpha}$. Here, we apply the coordinate representation
of the operators $X_{i}$ and $P_{i}$ expressed as

$\displaystyle X_{i}$ $\displaystyle=$ $\displaystyle x_{i},$ $\displaystyle
P_{i}$ $\displaystyle=$ $\displaystyle(\delta_{ij}+\alpha x_{i}x_{j})p_{j},$
(34)

where $x_{i}$ and $p_{j}$ satisfy the standard commutation relation of
ordinary quantum mechanics. This representation yields the following
commutation relation for the momentum operator

$\displaystyle[P_{i},P_{j}]=i\alpha(x_{i}p_{j}-p_{i}x_{j}).$ (35)

In the $X-P$ space, the modified volume element

$\displaystyle\frac{d^{3}Xd^{3}P}{(1+\alpha X^{2})^{3}},$ (36)

should be considered instead of $d^{3}xd^{3}p$ statEUP . We then proceed to
consider the consequences of such a modification in calculating the partition
and distribution functions. For a single particle, the partition function in
$X-P$ space can be found as

$\displaystyle Q_{1}^{{\rm EUP}}=\int\exp(-\beta H(P,X))\frac{d^{3}X\
d^{3}P}{(1+\alpha X^{2})^{3}},$ (37)

whence we get the corresponding thermal energy as

$\displaystyle U^{{\rm EUP}}=-\frac{\partial}{\partial\beta}\ln Q_{1}^{{\rm
EUP}},$ (38)

The above expression can also be combined with Eq. (37) to give Eq.(21). For
the free non-relativistic and relativistic particles, one finds

$\displaystyle Q_{1}^{{\rm NR,EUP}}=V_{{\rm eff}}(\alpha,r)\left(\frac{2\pi
m}{\beta}\right)^{\frac{3}{2}}=\frac{V_{{\rm eff}}(\alpha,r)}{V}Q_{1}^{{\rm
NR}},$ (39)

and

$\displaystyle Q_{1}^{{\rm R,EUP}}=\frac{V_{{\rm eff}}(\alpha,r)}{V}Q_{1}^{\rm
R}.$ (40)

respectively, where we have defined $V_{\rm
eff}(\alpha,r)=\int_{0}^{r}\frac{d^{3}X}{(1+\alpha X^{2})^{3}}$ as the
effective volume, and in the limit of $\alpha\rightarrow 0$, the usual volume
$V$ is recovered. Since $V_{\rm eff}(\alpha,r)$ is independent of $\beta$, the
thermal energy related to EUP is the same as what we obtained in Eqs. (7) and
(9), respectively. Consequently, EUP does not affect the Maxwell-Boltzmann and
Jüttner distribution functions, because the corresponding effective volume has
no dependence on $\beta$.

##  V Conclusion

The Jüttner function is the relativistic version of MB distribution and is
proper for studying relativistic (high energy) systems. On the other hand, the
minimal length comes into play in the realms of high energy physics. Hence,
compared with MB distribution, the study of its ef-fects on Jüttner
distribution would be more meaningful. Thus, our attempt in the present work
was to address an algorithm with the help of which, one can get the
distribution function, starting from the partition function. Motivated then by
the abovementioned arguments, we studied the effects of GUP and EUP (two
aspects of quantum gravity) on Jüttner distribution and the corresponding
density of states in energy space. We also addressed the consequence of
applying our approach to the MB distribution in order to find the density of
states Eq. (28) which is in agreement with previous reports maxwellGUP ;
phase1 , a result which confirms our approach. The results of our study are
summarized in Tables (1) and (2) for the non-relativistic and relativistic
regimes, respectively.

It is obvious from Figs. (1) and (2), that the effects of the existence of a
non-zero minimal length ($\eta\neq 0$) on distribution functions become more
sensible as energy increases. This means that the probability of achieving
high energy states when $\eta\neq 0$ is smaller than the $\eta=0$ case. It is
also worth mentioning that though there exist some proposals to test
observable effects of the minimal length Gustavo2019 , the Planck scale is
currently far beyond our reach. Since, by comparing the Planck energy
($\approx 10^{16}$TeV) Plenergy  to the energy achieved in the Large Hadron
Collider ($\approx 10$TeV) LHCenergy , or the Planck length ($\approx
10^{-35}$m) to the uncertainty within the position of the LIGO mirrors
($\approx 10^{-18}$m) LIGOlength  or the Planck time ($\approx 10^{-44}$s) to
the shortest light pulse produced in laboratory ($\approx 10^{-17}$s) Kun2012
, we observe that we are at best 15 orders of magnitude away from achieving
the Planck scale. In this regard, more future developments within these
experimental setups are expected in order to seek for the footprints of GUP
effects in nature.

Finally, regarding the results reported in Tsallis3  and Tsallis4  the
usefulness of Tsallis distribution function at high energy physics is
expected. In line with these results, some researchers study the possibility
of describing the distribution of transverse momentum in Large Hadron Collider
and Relativistic Heavy Ion Collider employing the Tsallis distribution,
expressed as Tsallis1 ; Tsallis11 ; Tsallis12 ; Tsallis2R

$\displaystyle f_{T}(q,\beta)=Z_{T}[1-(1-q)\beta E]^{\frac{1}{1-q}}.$ (41)

Here $Z_{T}$ and $q$ denote the normalization constant and non-extensivity
parameter, respectively. Although, utilizing our approach to investigate the
effects of GUP and EUP on (41) is straightforward, such a study needs more
careful analysis owing to the issues raised by Tsallis3  which states a
criterion on the domains of validity of Maxwell-Blotzmann, Jüttner and Tsallis
distributions as a special high energy phenomenon. Therefore, it can be
considered as an attractive topic for future studies.

Table 1: Non-relativistic ideal gas ($2\pi\hbar=1$) | HUP | GUP | EUP
---|---|---|---
The volume of phase space element | $1$ | $(1+\eta P^{2})^{3}$ | $(1+\alpha X^{2})^{3}$
Density of States | $\sqrt{E}$ | $\frac{\sqrt{E}}{(1+2m\eta E)^{3}}$ | $\sqrt{E}$
Single Partition Function | $V(\frac{2\pi m}{\beta})^{\frac{3}{2}}$ | $V(\frac{2\pi m}{\beta})^{\frac{3}{2}}I(\frac{2m\eta}{\beta},3)$ | $V_{\rm eff}(\alpha,r)(\frac{2\pi m}{\beta})^{\frac{3}{2}}$
Table 2: Relativistic ideal gas ($2\pi\hbar=1$) | HUP | GUP | EUP
---|---|---|---
The volume of phase space element | $1$ | $(1+\eta P^{2})^{3}$ | $(1+\alpha X^{2})^{3}$
Density of States | $E\sqrt{E^{2}-m^{2}}$ | $\frac{E\sqrt{E^{2}-m^{2}}}{(1+\eta(E^{2}-m^{2}))^{3}}$ | $E\sqrt{E^{2}-m^{2}}$
Single Partition Function | $V(\frac{2\pi m}{\beta})^{\frac{3}{2}}\Psi(\sigma)$ | $V(\frac{2\pi m}{\beta})^{\frac{3}{2}}\Psi(\sigma)(1-\eta\frac{15}{2}\frac{1}{\beta m}),\ \textmd{when}\ m\gg\frac{1}{\beta}$ | $V_{\rm eff}(\alpha,r)(\frac{2\pi m}{\beta})^{\frac{3}{2}}\Psi(\sigma)$

## References

  * (1) M. Maggiore, Phys. Rev. D 49, (1994) 5182; K. Konishi, G. Paffuti, P. Provero, Phys. Lett. B 234, (1990) 276; F. Scardigli, Phys. Lett. B, 452, (1999) 39.
  * (2) M. Maggiore, Phys. Lett. B 319, (1993) 83.
  * (3) R. Casadio, F. Scardigli, Phys. Lett. B 807, (2020) 135558; E. C. Vagenas, A. F. Ali, M. Hemeda, H. Alshal, EPJC, 79, (2019) 398; H. Shababi, W. S. Chung, Phys. Lett. B 770, (2017) 445; H. Shababi, P. Pedram, W. S. Chung, Int. J. Mod. Phys. A 31, (2016) 1650101; P. Pedram, Phys. Rev. D 91, (2015) 063517; H. S. Snyder, Phys. Rev. 71, (1947) 38; C. N. Yang, Phys. Rev. 72, (1947) 874.
  * (4) Z. W. Feng, H. L. Li, X. T. Zu, S. Z. Yang, Eur. Phys. J. C 76, (2016) 212; A. N. Tawfik, A. M. Diab, Rep. Prog. Phys. 78, (2015) 126001; M. M. Khalil, AHEP 2014, (2014) 619498; S. Miraboutalebi, L. Farhang Matin, Can. J. Phys. 93, (2014) 574; A. F. Ali, M. Moussa, AHEP 2014, (2014) 629148; M. A. Motlaq, P. Pedram, J. Stat. Mech. (2014) P08002; S. Das, E. C. Vagenas, Phys. Rev. Lett. 101, (2008) 221301; S. K. Rama, Phys. Lett. B 519, (2001) 103; A. Kempf, J. Phys. A: Math. Gen. 30, (1997) 2093.
  * (5) D. Park, Phys. Rev. D 101, (2020) 106013; P. Bosso, S. Das, I. Pikovski, M. R. Vanner, Phys. Rev. A 96, (2017) 023849; I. Pikovski, M. R. Vanner, M. Aspelmeyer, M. S. Kim, C. Brukner, Nature Phys. 8, (2012) 393; S. Das, R. B. Mann, Phys. Lett. B 704, (2011) 596.
  * (6) G. G. Luciano, L. Petruzziello, Eur. Phys. J. C 79, (2019) 283; G. Gecim, Y. Sucu, Phys. Lett. B 773, (2017) 391; V. Husain, S. S. Seahra, E. J. Webster, Phys. Rev. D 88, (2013) 024014; W. Chemissany, S. Das, A. F. Ali, E. C. Vagenas, JCAP 12 (2011) 017; M. Sprenger, M. Bleicher, P. Nicolini, Class. Quant. Grav 28, (2011) 235019; T. Zhu, J. R. Ren, M. F. Li, Phys. Lett. B 674, (2009) 204.
  * (7) J. R. Mureika, Phys. Lett. B 789, (2019) 88; W. S. Chung, H. Hassanabadi, Mod. Phys. Lett. A 32, (2017) 1750138; S. Mignemi, Mod. Phys. Lett. A 25, (2010) 1697; C. Bambi, F.R. Urban, Class. Quantum Gravity 25, (2008) 095006.
  * (8) A. Kempf, J. Math. Phys. 38, (1997) 1347; H. Hinrichsen, A. Kempf, J. Math. Phys. 37, (1996) 2121.
  * (9) M. P. Dabrowski, F. Wagner, Eur. Phys. J. C 79, (2019) 716.
  * (10) G. Lambiase, F. Scardigli, Phys. Rev. D 97, (2018) 075003.
  * (11) D. Walker, J. H. Fremlin, Nature, 171, (1953) 189.
  * (12) T. Kalaydzhyan, Phy. Lett. B 750, (2015) 112; G. A. Camelia, Living Rev. Relativ. 16, (2013) 5.
  * (13) L. Chang, D. Minic, N. Okamura, T. Takeuchi, Phys. Rev. D 65, (2002) 125028.
  * (14) T. V. Fityo, Phys. Lett. A 372, (2008) 5872.
  * (15) P. Wang, H. Yang, X. Zhang, JHEP 2010, 43 (2010).
  * (16) S. Hossenfelder, Living Rev. Relativ. 16, (2013) 1-90.
  * (17) A. E. Shalyt-Margolin, A. Y. Tregubovich, Mod. Phys. Lett. A 19, (2004) 71; A. E. Shalyt-Margolin, J. G. Suarez, Int. J. Mod. Phys. D 12, (2003) 1265.
  * (18) B. Vakili, M. A. Gorji, Journal Statistical Mechanics: Theory and Experiment, 2012 (2012) P10013.
  * (19) M. Mirtorabi, S. Miraboutalebi, A. A. Masoudi, L. Farhang Matin, Physica A, 506 (2018) 602.
  * (20) W. S. Chung, H. Hassanabadi, Int. J. Mod. Phys. A 34, (2019) 1950041; K. Nozari, A. Etemadi. Phys. Rev. D 85, (2012) 104029.
  * (21) F. Jüttner, Ann. Phys. 34, (1911) 856.
  * (22) R. K. Pathria, P. D. Beale, “Statistical Mechanics,”3rd ed.; Elsevier: Burlington, MA, USA, 2011.
  * (23) W. Pauli, “The Theory of Relativity,” Pergamon, Pergamon Press: London, UK, 1958.
  * (24) L. P. Horwitz, W. C. Schieve, C. Piron, Ann. Phys. 137, (1981) 306.
  * (25) S. R. de Groot, W. A. van Leeuwen, Ch. G. van Weert, “Relativistic Kinetic Theory-Principles and Applications,” North-Holland: Amsterdam, The Netherlands 1980.
  * (26) G. Livadiotis, Ann. Geophys., 34, (2016) 1145.
  * (27) D. Cubero, J. Casado-Pascual, J. Dunkel, P. Talkner, P. Hänggi, Phys. Rev. Lett. 99, (2007) 170601.
  * (28) A. Montakhab, M. Ghodrat, M. Barati, Phys. Rev. E 79, (2009) 031124.
  * (29) M. Ghodrat, A. Montakhab, Phys. Rev. E 82, (2010) 011110.
  * (30) J. Dunkel, P. Talkner, P. Hünggi, New J. Phys. 9, (2007) 144.
  * (31) A. Kempf, G. Mangano, R. B. Mann, Phys. Rev. D 52, (1995) 1108.
  * (32) A. E. Rastegin, Entropy 2018, 20(5), 354.
  * (33) A. Tawfik and A. Diab, Int. J. Mod. Phys. D 23, 1430025 (2014).
  * (34) C. O’Luanaigh, CERN, Available online:https://home.cern/news/news/accelerators/first-successful-beam-record-energy-65-tev (accessed on 10 April 2015).
  * (35) B. P. Abbott et al. Phys. Rev. Lett. 116, 061102 (2016).
  * (36) K. Zhao, Q. Zhang, M. Chini, Y. Wu, X. Wang, and Z. Chang, Optics Letters 37, (2012) 3891.
  * (37) B. Walton, J. Rafelski, Phys. Rev. Lett 84, (2000) 31 .
  * (38) U. Tirnakli, E. Borges, E. Sci. Rep. 6, 23644 (2016).
  * (39) A. S. Parvan, Eur. Phys. J. A 56, (2020) 106.
  * (40) F. H. Liua, Y. Q. Gao, B. C. Li, Eur. Phys. J. A 50, (2014) 123.
  * (41) R. F. Si, H. L. Li, F. H. Liu, AHEP 2018, 2018, 1-12.
  * (42) H. Zheng, Lilin Zhu and A. Bonasera, Phys. Rev. D 92, 074009 (2015).