Abstract
Recent studies reveal that Allee effect may play important roles in the growth of tumor. We present one of the first mathematical models of avascular tumor that incorporates the weak Allee effect. The model considers the densities of tumor cells in three stages: proliferating cells, quiescent cells, and necrotic cells. We investigate how Allee effect impacts the growth of the avascular tumor. We also investigate the effect of apoptosis of proliferating cells and necrosis of quiescent cells. The system is numerically solved in 2D using different sets of parameters. We show that Allee effect and apoptosis play important roles in the growth of tumor and the formation of necrotic core.
1. Introduction
Genetic mutation in DNA modifies the normal cell proliferation rate and causes a normal cell to become cancerous. A cancer cell then searches for nutrient such as oxygen and glucose from nearby tissues and proliferates rapidly. This results in a clump of cancer cells known as tumor. The growth of a tumor depends on many factors such as its location, cell type, and also the availability of nutrients. There are three stages of tumor growth, namely, avascular stage, angiogenesis stage, and metastasis stage ([1]). At avascular stage, the tumor often grows in an almost spherical form (Figure 1). At this stage, the tumor is harmless and may reach a saturation size due to the limited nutrient supply to the core of the tumor. However, nutrient-deficient tumor cells produce signaling proteins such as vascular endothelial growth factor (VEGF). When enough such signaling proteins are produced, the nearby blood vessels will provide additional nutrient through neovasculature to sustain the growth of the tumor, leading to the second stage called angiogenesis. Once the tumor makes connection with nearby blood vessels, mobile tumor cells can migrate to the other parts of the body and develop additional tumors at different sites, resulting in the metastasis stage ([2]).
(a)
(b)
In experiments, tumor cell lines may be grown in many different ways, for example, as a two dimensional monolayer on a substrate, as a three dimensional spheroid suspended in a liquid or gel, and as a mouse xenograft. Each of these methods yields different information about the cell line. For example, it is known that the response of the tumor culture to treatment differs between monolayer and spheroid cultures ([3]), and mathematical models exhibit this differential response also ([4]). Glioma cells grown as spheroids produce more lactic acid, use more glucose, and exhibit more markers of hypoxia than those grown in monolayer.
There are multiple methods for growing spheroids ([5]). Tumor cell lines that are adhesive enough to be grown as in vitro spheroids always exhibit a necrotic core ([6]). It is commonly believed that the necrotic core is formed when the oxygen drops to a critical level, although it is also known that a feedback mechanism due to the production of tumor necrosis factor alpha causes apoptosis in the proliferating cell population ([7, 8]). Both mechanisms are supported by the observation that when the tumor spheroids grown in vitro have a diameter reaching 500 , a necrotic core will form ([9]), and their growth eventually ceases ([10]). Both mechanisms are also supported by assorted simple mathematical models ([11]). In support of factors besides hypoxia, some researchers observed that necrotic core can also form in relatively high oxygen level ([12, 13]). In ([13]), Muller-Klieser pointed out that necrotic core is not necessarily a consequence of hypoxia or anoxia but emerges as an accumulation of apoptosis. Apoptotic fragments may accumulate in the tumor, leading to the formation of a necrotic core. They observe that the diameter of the necrotic core expands at approximately the same rate as the spheroid volume increases. This results in a relatively constant thickness of the viable cell rim with increasing spheroid diameter ([12]). In this article, we consider both hypoxia-induced apoptosis and necrosis as mechanisms that lead to the formation of central necrotic core, but leave the role of tumor necrosis factor alpha for a future study, although it is known that glioma cells are sensitive to TNF-alpha ([14]). See Figure 2 for a comparison between apoptosis and necrosis.
In recent years, there has been a growing interest in viewing tumors and its microenvironment as an ecosystem ([15, 16]). A cancer ecosystem contains multiple strains of tumor cells as well as epithelial cells, nutrients, growth factors, etc. This perspective emphasizes their interactions and consequences ([17]) instead of treating cancer as a collection of muted cells. Our knowledge of ecological dynamics within tumors is limited, especially in small tumors prior to detection or after treatment ([16]). However, even under the assumption that a single cell line is grown in a biologically inert substrate, ecological phenomena may occur, in particular, the Allee effect. First described in the 1930s, Warder Allee was able to show that, contrary to intuition, goldfish grow more rapidly in a tank with more goldfish ([18]). Since then, ecologists and mathematicians have studied this effect, now termed Allee effect, in many different situations. A strong Allee effect describes a population that can grow at intermediate population density but declines when the number of organisms is either too small or too large. Such populations are then likely to collapse and become extinct if their population size falls below a certain threshold. Clearly, any cancer cell lines that can propagate indefinitely in vitro are highly unlikely to be subject to a strong Allee effect. However, a weak Allee effect describes a population that grows with a reduced per capita growth rate at a low population size, although the growth rate is still positive. For example, low-grade cell line cultures, in vivo and in vitro, have low chances of persistence and low reproducibility. On the contrary, high-grade cell lines have a higher chance to establish ([19, 20]). The possible existence of a weak Allee effect in tumor cells has been described in a very recent study ([21]). The proliferation rate of tumor cells at low density was measured via a series of in vitro experiments with two glioblastoma cell lines ([22]), where it was found that when the cell density is low, the growth rate increases with the population density, while decreasing at larger density. This study, which suggests a weak Allee effect in glioblastoma, is echoed in some cellular automata models of spheroid growth, in which low density affects slow proliferation ([23]). In cultures of cancer cells weak Allee effect likely arises due to autocrine growth factors, diffusive signaling molecules produced and secreted by cells that enhance growth and proliferation of other cells ([22]).
The study of tumor growth using mathematical models has been in focus for the past four decades or so. Starting from early diffusion and differential equation models by Thomlinson and Gray, which was then extended by Burton and Greenspan, tumor models have expanded to include ordinary differential equation (ODE) models, partial differential equation (PDE) models, cellular automata (CA) models, and the more sophisticated hybrid models that combines PDE models and CA models ([24–26]). For a complete review on mathematical modeling of avascular tumors, we refer the readers to the review by Araujo and McElwain ([27]) or the more recent one by ([28]). Here we describe briefly a handful of models closely related to the one presented in this study.
A malignant tumor invasion model with strong Allee effect was recently studied in ([29]), focusing on properties of travelling wave solutions. A PDE system comprising both live-cell density, dead-cell density, and nutrient density was developed in ([30]). However, they did not distinguish between proliferating and quiescent cells, nor did they consider different types of cell death (via necrosis or apoptosis). Later they built upon this model by incorporating two necrotic depletion mechanisms: leakage of the dead cellular material by diffusion or consumption of cellular material for the construction of new cells ([31]).
Sherratt and Chaplain proposed an interesting model that accounts for three different states of tumor cells ([32]). In their model, the necrotic core comes directly from the necrosis of quiescent cells and there is no weak Allee effect. Their model also does not consider the fact that quiescent cells consume nutrients as well, only at a slower rate. Including cell movement due to contact inhibition, they considered the following reaction diffusion model in 1D:
The existence of travelling wave solutions for some special cases was studied in Zhu and Ou [33] using the Banach fixed point theorem and perturbation methods. Travelling wave solutions to similar systems for contact inhibition were also studied in Bertsch et al. [34]. The impact of Allee effect on the spreading of species has been well studied in ecological modeling. It is generally believed that Allee effect can adversely affect the spreading speed. The Allee effect can also explain the time lag that is often noticed between the time when a species first establishes and the time it starts to spread rapidly ([35]). Similarly, an Allee effect caused a long time lag between the resection and recurrence in a one dimensional spatial model glioblastoma ([22]).
The early study of the travelling wave solution of the following reaction diffusion equation with Allee effect can be found in ([36, 37]). represents the Allee threshold below which the growth rate of the species becomes negative. The main results can be summarized as follows: (a) For , there exists a unique wave front solution with speed . For , the front propagates to the region where the species is absent, i.e., the invasion of the species; for , the front propagates to the region where the species is at its carrying capacity, i.e., the retreat of the species; when , the front is stationary.
(b) For , there exists a maximal wave speed such that the travelling wave solution exists when . Moreover, for (derived by linearizing about zero) and when .
A weak Allee effect in a highly necrotic expanding spheroid would necessarily apply to any spatial region where the density of live cells is quite low. This is a situation where mathematical modeling can reveal unknown or counter-intuitive principles that may have been overlooked before. It can test theories by simplifying the underlying mechanisms ([38]). Our model builds on prior ones by extending the spatial domain to two dimensions, including diffusion of nutrient, motion of cells, the weak Allee effect, the quiescent cell compartment, and two mechanisms of cell death: apoptosis and necrosis. Our hypothesis is that, for expanding spheroids in substrate suitable for migration, a weak Allee effect will slow the rate of expansion of the proliferating edge, as measured by the speed of the wavefront. To capture the spatial characteristics of central necrosis, a ring of quiescent cells, and a proliferating rim developing in response to nutrient delivery and depletion we propose a system of PDEs in two spatial dimensions, and with weak Allee effect, for the purpose of computationally testing this hypothesis.
2. Mathematical Model
Consider the following model: subject to the initial conditions and boundary conditions
Here represents the density of proliferating cells, the density of quiescent cells, and the density of necrotic cells. denotes the concentration of oxygen.
The functions and will be decreasing functions of with . is an increasing function of . will be increasing with . is the apoptotic rate at which proliferating cells die. For simplicity, we assume is a constant. The exact forms of these switching functions remain unknown. Many authors have considered the simplest case, for example, a linear switching rate between proliferating cells and quiescent cells. is the rate at which proliferating cells become quiescent. We use a Michaelis-Menten-type expression where is the critical oxygen level and is positive constant which governs the sharpness of change near the critical concentration. is the switching rate when the oxygen is 0. In the numerical simulations that follow, we take as the standard value. A typical profile of is presented in Figure 3.
For the switching rate from quiescent cells back to proliferating cells, we take where is switching rate when is high. In the simulations that follow, we take as the standard value. See Figure 4 for a typical profile of .
For the apoptotic rate, we let where is the apoptosis rate of proliferating cells.
For the necrosis rate of quiescent cells, we adopt the following form where is the maximum necrosis rate of quiescent cell. is the critical oxygen level for necrosis. See Figure 5.
We let with .
As , tends to a Heaviside function. These two choices of switching function do not significantly alter the outcome as their profiles can become very similar when is large.
The oxygen diffuses according to a diffusion coefficient . The degradation rate of oxygen is denoted by and it is utilized by proliferation cells at a rate . The rate that is being utilized by quiescent cells is denoted by . Throughout this paper, is assumed to be , within the experimentally determined range Freyer et al. [39].
Given that the diffusion rate of oxygen throughout the tumor spheroid is much faster comparing to the scale of growth, we adopt the standard quasi steady-state assumption that where , , and are dimensionless parameters. For simplicity, we assume .
Our model depends crucially on the biological parameters. These parameters vary greatly depending on tumor type, localization, progression, etc. Since our goal is to investigate the impact of different parameters on the growth of avascular tumors, parameters are chosen within the experimentally observable range rather than a quantitative description of any specific tumor. Most of the parameters involved are available in literature ([40, 41]). Our system is already given in nondimensional form for computational reasons. We take characteristic scale of time and length second and cm, respectively. These nondimensional parameters are summarized in Table 1.
3. Impact of Allee Effect on Travelling Wave Solution and Wave Speed
In Petrovskii et al. [42], the authors considered how Allee effect impacts different scenarios of biological invasion in a predator-prey system. They noticed that as a result of the interplay between the Allee effect and the predation, species extinction and unbounded spatial spread are not the only possible system dynamics. For certain parameter values, initial species distribution can evolve to quasi-stationary patches. The position of these patches does not change with time while their shape is either stationary or oscillatory. They pointed out that invasive species can be held localized due to certain interspecific and intraspecific interactions such as the interplay between the Allee effect and predation.
We consider travelling wave solutions of the following form: where and is the wave speed. Substituting these solutions into the original system yields the following ODEs As usual, we linearize the first equation of this system about and assume , , and to leading orders ([37]); we have Thus for nontrivial solutions, we must have Thus positive solution exists if and only if where Thus the minimal speed (if exists) for our model must satisfy .
It remains unclear if linear determinacy holds for this problem. By linear determinacy, we mean .
4. Numerical Results and Discussion
All our numerical simulations employ the zero-flux boundary conditions with a system size 210210 discretized through and , with . Our simulations use forward difference in time and central difference in space with a small time step and space step . The concentration at mesh point at time is denoted by and they are given by with and defined by
Here , , and are given by
The initial conditions are set as follows:
4.1. Allee Effect
It is now well known that the growth curves for multicellular spheroids are characterized by doubling time which decreases with increasing growth until a stable saturation was attained ([7]). However, the original model proposed in ([32]) does not support this behavior. The solutions exhibit travelling wave solutions and the tumor continues to grow outward. In this section, we show that Allee effect plays an important role in the growth of avascular tumor. Essentially, as the density of proliferating cells drops, the weak Allee effect further slows down its growth, which overall slows down the growth of the tumor.
In Figure 6, we compare typical profiles of proliferating cells, quiescent cells, and necrotic cells at . Starting with the same initial data, both models exhibit three layers of different cells, with a necrotic core inside. With Allee effect, the tumor size is significantly reduced.
(a) Without Allee effect
(b) With Allee effect
(c) Heat map of P without Allee effect
(d) Heat map of P with Allee effect
Figure 7 shows individual profiles for proliferating cells, quiescent cells, and necrotic cells without Allee effect at time .
(a) Proliferating cells
(b) Quiescent cells
(c) Necrotic cells
In Figure 8, we show how the strength of Allee effect impacts the size of the tumor. Figure 8(a) shows the density of proliferating cells at without Allee effect. Figure 8(b) shows the density of proliferating cells at with weak Allee effect (). When the strength of Allee effect is stronger (), there is a significant reduction in tumor size as shown in Figure 8(c).
(a) No Allee effect
(b) Weak Allee effect
(c) Stronger Allee effect
Figure 9 shows the impact of Allee effect on the profile of necrotic core.
(a) No Allee effect
(b) Weak Allee effect
(c) Stronger Allee effect
In Figure 10, we study the impact of Allee effect on the size of the tumor at different time. We locate the peak of the proliferating cell density and use its position to estimate the radius of the tumor. The results are shown in Figure 10. We note that the radius of tumor grows linearly without Allee effect. It is significantly reduced when weak Allee effect is present. When the strength of Allee effect is stronger, however, the radius grows exponentially at the very beginning, then quickly becomes linearly, and eventually reaches a plateau.
4.2. Impact of Apoptotic Cell Death Rate
This section considers the effect of cell death from both the quiescent cells (necrosis) and proliferating cells (apoptosis). Comparing to small apoptosis rate (), we increase this rate to . The profiles of different tumor cells at different times are shown in Figure 11. Higher level of apoptosis leads to fewer quiescent cells and more necrotic cells. The rim of proliferating cells still remains on the boundary. A closer look at the concentration of nutrient reveals that a higher level of apoptosis leads to higher nutrient concentration at the center of the spheroid. This is due to the fact that a higher level of apoptosis increases the amount of necrotic cells at the center, thus reducing the amount of viable cells.
(a) Large apoptotic cell death at
(b) Large apoptotic cell death at
(c) Large apoptotic cell death at
(d) Large apoptotic cell death at
4.3. Impact of Proliferation Rate
In this section, we simulate the system with a smaller proliferation rate and a larger proliferation rate . A larger proliferation rate leads to a higher density in proliferating cells which will also lead to more quiescent cells and necrotic cells. The result is shown in Figure 12. For a better illustration, please see the heat map in Figure 13.
(a) Small proliferation rate
(b) Standard proliferation rate
(c) Large proliferation rate
(a) Proliferating cells
(b) Quiescent cells
(c) Necrotic cells
(d) Proliferating cells
(e) Quiescent cells
(f) Necrotic cells
4.4. Impact of Oxygen Degradation Rate
To study the effects of oxygen on the growth of the avascular tumor, we compare tumor cell density profiles at different oxygen degradation rate at and at time . The spreading speed of the tumor increases as the degradation rate gets higher. At , the radius of the tumor is estimated to be 27 when , comparing to 30 when (Figure 14). As the degradation rate of oxygen increases, we also notice a significant reduction on necrotic region. For a better illustration, please see the heat map of all the cells in Figure 15.
(a) Large oxygen degradation rate
(b) Small oxygen degradation rate
(a) Proliferating cells
(b) Quiescent cells
(c) Necrotic cells
(d) Proliferating cells
(e) Quiescent cells
(f) Necrotic cells
4.5. Impact of Critical Oxygen Level
In this section we observe a very important relationship between the growth rate of the tumor and the critical oxygen level for proliferation. To study the impact of critical oxygen level on the growth of the tumor, we calculate the tumor radius at for different critical values. The result is shown in Figure 16. The simulations show that the tumor growth speed decreases dramatically as the critical oxygen level increases.
4.6. Impact of Crowding Factor
This section compares the growth of the tumor layers for different crowding factors. A larger crowding factor , i.e., increasing the cell crowding, lowers the availability of oxygen. As a result, the tumor has a smaller density of proliferating cell and quiescent cell. The necrotic core of the tumor, however, becomes denser. Figure 17 shows the tumor profiles at for and , respectively.
(a) Small crowding factor
(b) Large crowding factor
5. Summary and Discussions
This work is the first 2D mathematical model that incorporates Allee effect in the growth of avascular tumor. In the present paper we study the impact of various parameters on tumor progression. This work further expands our understanding of the process of avascular tumor growth by considering weak Allee effect and apoptosis as an additional mechanism of cell death. The extension to 2 spatial dimensions allows a more precise description of the processes of avascular tumor growth. One of our main findings is that weak Allee effect may lead to a plateau size. Another finding is that apoptosis increases the size of necrotic core and reduces the amount of viable cells. Additional numerical findings can be summarized as follows: (a) Large proliferation rate increases the overall size of tumor. (b) Large degradation rate of oxygen significantly reduces the size of necrotic region. (c) Increasing the critical oxygen concentration reduces the growth rate of the tumor. We point out that our aim in this work is not to propose completely new mathematical models, but to investigate two new mechanisms (Allee effect and apoptosis) through numerical simulations of Owen-Sherratt model. It is surely more realistic to consider the spheroid in 3D. But the main ingredient of the current research relies on the assumption that the nutrient reaches the core of the avascular tumor much faster than the reaction time of proliferating cells, etc. For 3D spheroid, this assumption no longer holds.
The impact of weak Allee effect on the growth of avascular tumors is investigated in this article. The understanding of this early stage is relevant to progression, relapse, and metastasis. A model that considers weak Allee effect may help to improve the predictions of tumor growth and relapse dynamics.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Disclosure
The result in this article has been presented at 38th Southeastern-Atlantic Regional Conference on Differential Equations.
Conflicts of Interest
The authors have no conflicts of interest.