Next Article in Journal
Immunoassay Urine Drug Testing among Patients Receiving Opioids at a Safety-Net Palliative Medicine Clinic
Previous Article in Journal
Enhancing Prostate Cancer Diagnosis with a Novel Artificial Intelligence-Based Web Application: Synergizing Deep Learning Models, Multimodal Data, and Insights from Usability Study with Pathologists
Previous Article in Special Issue
Exploring Cell Migration Mechanisms in Cancer: From Wound Healing Assays to Cellular Automata Models
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Hybrid Cellular Automata Modeling Reveals the Effects of Glucose Gradients on Tumour Spheroid Growth

1
Dipartimento di Ingegneria Chimica, dei Materiali e della Produzione Industriale, Università degli Studi di Napoli Federico II, 80125 Naples, Italy
2
CEINGE-Biotecnologie Avanzate Franco Salvatore, Via G. Salvatore 436, 80131 Naples, Italy
3
Mathematics in Medicine Program, Houston Methodist Research Institute, Houston, TX 77030, USA
4
Department of Physiology and Biophysics, Weill Cornell Medical College, New York, NY 10065, USA
5
Neal Cancer Center, Houston Methodist Research Institute, Houston, TX 77030, USA
6
Department of Imaging Physics, University of Texas M.D. Anderson Cancer Center, Houston, TX 77030, USA
7
Physiology, Biophysics, and Systems Biology Program, Graduate School of Medical Sciences, Weill Cornell Medicine, New York, NY 10065, USA
*
Authors to whom correspondence should be addressed.
These authors contributed equally to this work.
Cancers 2023, 15(23), 5660; https://doi.org/10.3390/cancers15235660
Submission received: 12 October 2023 / Revised: 15 November 2023 / Accepted: 28 November 2023 / Published: 30 November 2023
(This article belongs to the Special Issue Experimental and Modeling Efforts to Target Metabolism in Cancer)

Abstract

:

Simple Summary

In recent years, mathematical models have revolutionized cancer research, illuminating the complex dynamics of tumor growth and aiding drug development. These models, reflecting biological and physical processes, are increasingly used in clinical practice, offering precise patient-specific predictions. Our work introduces an innovative in silico model to simulate tumor growth and invasiveness. The automated hybrid cell, replicating key tumor cell features, enables exploration of 3D tumor spheroid evolution. Sensitivity analyses reveal that tumor growth is primarily influenced by cell replication speed and adhesion, while invasiveness relies on chemotaxis. These insights shed light on tumor development mechanisms, guiding effective strategies against tumor progression. Our model serves as a valuable tool for advancing cancer biology research and potential therapeutic interventions.

Abstract

Purpose: In recent years, mathematical models have become instrumental in cancer research, offering insights into tumor growth dynamics, and guiding the development of pharmacological strategies. These models, encompassing diverse biological and physical processes, are increasingly used in clinical settings, showing remarkable predictive precision for individual patient outcomes and therapeutic responses. Methods: Motivated by these advancements, our study introduces an innovative in silico model for simulating tumor growth and invasiveness. The automated hybrid cell emulates critical tumor cell characteristics, including rapid proliferation, heightened motility, reduced cell adhesion, and increased responsiveness to chemotactic signals. This model explores the potential evolution of 3D tumor spheroids by manipulating biological parameters and microenvironment factors, focusing on nutrient availability. Results: Our comprehensive global and local sensitivity analysis reveals that tumor growth primarily depends on cell duplication speed and cell-to-cell adhesion, rather than external chemical gradients. Conversely, tumor invasiveness is predominantly driven by chemotaxis. These insights illuminate tumor development mechanisms, providing vital guidance for effective strategies against tumor progression. Our proposed model is a valuable tool for advancing cancer biology research and exploring potential therapeutic interventions.

1. Introduction

According to the World Health Organization, cancer remains one of the leading causes of death in advanced countries. In 2020 alone, there were 19.3 million new cases and 10 million deaths reported worldwide [1,2]. However, there has been a consistent decrease in cancer-related mortality, primarily due to successful prevention efforts and advancements in treatment options [3]. It is noteworthy that approximately 66.7% of cancer deaths are related to a process known as metastasis [4], which refers to the formation of secondary tumors in parts of the body different from the site of cancer origin. The process involves cancer cells crossing endothelial walls, circulating through in the bloodstream, and eventually extravasating from capillaries crossing the basement membrane to colonize the new site [5].
A key feature of the abovementioned steps is the ability of cancer cells to migrate directionally in response to external stimuli, known as chemotaxis, which plays a fundamental role in enabling tumor invasiveness [6]. Chemotaxis is defined as the directional movement of an organism in response to the concentration gradient of a given chemical species. In the context of cell chemotaxis, it refers to the sensitivity of cells to specific gradients, quantified as the ratio between the concentration gradient and the local concentration value ( S G = C / C ) [7,8]. Chemokines (e.g., interleukin) and growth factors (e.g., FBS) are examples of chemical species that drive chemotaxis in cancer cells [7]. Metabolites, such as glucose or oxygen, which are also present in many growth factors, as well as catabolite gradients, can also induce chemotaxis.
In the case of cancer tissues, cell over-proliferation induces a lack of nutrients and an accumulation of catabolites. Consequently, tumor masses lose compactness and, in order to keep a high surface to volume ratio, tend to invade the surrounding tissues. This morphologic instability [9], known as diffusional instability [10,11], can be predicted by mathematical models. The investigation of such a complex phenomenon requires the use of adequate models able to include both biological and transport phenomena aspects. Different approaches have been followed in cancer research, ranging from simple 2D cell cultures, complex 3D scaffolds in vitro, murine in vivo models, to clinical studies.
In cancer research, in vitro approaches, which involve the study of biological systems in artificial laboratory conditions, have been widely employed to investigate the underlying mechanisms of cell migration, proliferation, and other processes related to tumorigenesis. Most of the published data regarding known cell-based processes are derived from experiments performed in two-dimensional (2D) conditions. Two-dimensional cell cultures, growing on solid substrates, such as plastic, do not fully reproduce the complex three-dimensional (3D) architecture and complexity of living tissues. Important aspects, such as cell–cell interactions, tissue phenotypes, the role of cell density and extracellular matrix (ECM) [12,13,14,15], proliferation regulators [16], and metabolic functions [17] are often missing or not adequately represented in 2D cultures.
In order to mimic the native in vivo scenario [18,19,20] where cancer growth occurs, 3D models have been used in cancer research as a compromise between 2D cell cultures and whole-animal systems. However, 3D tumor microenvironments in murine models may not accurately represent the human scenario, making it challenging to control and interpret experimental outcomes.
Recently, the tumor spheroid [8], a tightly bound aggregate of cells, is gaining popularity as a 3D model due to of its ability to strikingly mirror the 3D cellular context. Being characterized by naturally physiological and chemical gradients, cellular spheroids consist of actively proliferating cells on the outside with quiescent cells in the inner rim, and a central nutrient-deprived zone (necrotic core) [21], mimicking the natural scenario in a tumor in situ.
Each of the abovementioned approaches has limitations that can be taken into account by adequately coupling experimental investigations with in silico mathematical models. One of the main challenges is to capture the complex and dynamic interactions between cancer cells and their microenvironment. In vitro experimental conditions may introduce artifacts that can potentially compromise the validity of conclusions drawn from these studies if not properly accounted for. In silico approaches, based on mathematical modeling and computer simulations, have the potential to overcome these limitations by focusing on the representation of the interactions of cells with the microenvironment, which are difficult to mimic experimentally [22,23]. Furthermore, in silico models, once developed and adequately validated, can be used to conduct extensive virtual experimental campaigns at a significantly lower cost compared to in vivo and in vitro research. This approach is crucial for the fundamental understanding of mechanisms of cancer growth and can facilitate the simulation of specific treatments on individual patients, or the study of the effects of rare mutations or genetic variations, implementing precision and personalized medicine [24] protocols.
For this reason, in silico approaches [25] are attracting more attention. The earliest attempts to mathematically model tumor growth and invasiveness date back to 1967, and since then, the number and accuracy of these models have continued to grow [14]. These models can be classified into continuum, discrete, and hybrid models [26], and they are employed to simulate and study a wide variety of phenomena, such as the resistance of tumors to different drug treatments [27,28,29], the interaction between tumor cells and their microenvironment [13,15], and the immune system [25,30].
Hybrid agent-based models, such as hybrid cellular automata (HCA), have proven to be effective in investigating complex tumor systems by combining deterministic reaction–diffusion partial differential equations with the representation of the cells as single and autonomous entities (cellular automata, CA [25]) governed by deterministic or stochastic laws. These models consider the interactions between individual cells [26,27,28,29,31] and can provide greater insights into the mechanisms underlying cancer growth and progression, supporting the design of new therapeutic strategies [30,32].
HCA models are particularly useful in many applications in cancer research, where a multi-scales model is required to depict the dynamics of cancer development over time, including the evolution of cell phenotypes and genotypes [25]. These models account for the interaction of individual cells with parameters characterizing the surrounding environment, such as concentration fields or changes in pH [33].
In this study, a computational model based on HCA is developed to simulate the growth and invasiveness of spheroids under different gradients of nutrients. Specifically, chemotaxis is taken into account by estimating the glucose concentration profiles surrounding a cancer cell aggregate (tumor spheroid). The model aims to investigate the roles of cell migration, duplication, and other parameters, such as cell–cell adhesion and sensitivity to chemotactic gradients, in the phenomenon of cancer invasion. A global sensitivity analysis (GSA) and a local sensitivity analysis (LSA) are performed to evaluate the model’s sensitivity to the variations in the abovementioned cellular parameters.
The goal of this study is to develop a model that can capture the complexity of the cancer invasion mechanism, as envisioned by the diffusional instability model [9], while maintaining computational feasibility. This will provide a more comprehensive understanding of cancer biology and aid in the design of new therapeutic strategies.

2. Materials and Methods

In the following section, the structure of the hybrid cellular automata model, referred as 2D-HCA, used for simulating the in vitro growth of avascular tumor spheroids, is presented in detail. The model description is organized into several subsections. Firstly, the model domain (Section 2.1.1) and numerical methods used to calculate continuous functions (in our case, glucose concentration) (Section 2.1.2) are presented. The HCA model is described by presenting the cell phenotypes (Section 2.1.3) and rules governing cellular dynamics (Section 2.1.4). Statistical and sensitivity analyses are described in an independent subsection (Section 2.2).
By structuring the model description in this manner, the paper provides a comprehensive overview of the 2D-HCA model and its application in simulating the in vitro growth of avascular tumor spheroids.

2.1. Model Development

The 2D domain used in the model represents a layer of the extracellular matrix (ECM) that contains a single tumor spheroid. This domain is discretized into a squared grid, referred to as a lattice. Each element of the lattice (lattice point, LP) can be occupied by a single cell or remain empty, representing the presence of ECM only. The entire grid evolves through a series of discrete time steps, following a set of pre-defined rules. The evolution of each cell is governed by the chemical stimuli in the LP occupied, and by the state of neighboring LPs. In this model, the chemical stimuli are represented by the concentration profile of a chemoattractant, specifically glucose. The model can be viewed as a superposition of two identical square grids, where one represents the cellular layer (see Section 2.1.3) and the other is used to calculate the glucose concentration field (see Section 2.1.2). The model can be easily generalized to account for more chemical species, or other types of stimuli, such as pressure fields, simply adding further layers.
In Figure 1a, a schematic representation (not in scale) of the cellular layer superimposed onto the glucose layer of our model is presented. In the cellular layer, the LPs are colored either in orange or blue, corresponding to the positions occupied or not by the cells. In each LP, the glucose concentration is calculated and schematically represented in Figure 1a in a color scale. As qualitatively reported in Figure 1a, the glucose concentration at the grid margins is higher compared to the center of the domain. The resulting symmetrical radial concentration gradient is achieved due to an isotropic source of glucose from each of the four edges of the domain, and a consumption in the center where the cell spheroid is located.

2.1.1. Domain Building

In this study, the domain considered for the model was a squared layer with dimensions of 2 mm × 2 mm. This domain was discretized into a grid consisting of 100 × 100 LPs. Each LP represented a cell with an approximate diameter of 20 μm. Two scenarios were examined in this study, defined as isotropic and gradient. Both scenarios shared the same set of parameters, except for the initial and boundary conditions in the glucose layer.
In the isotropic case, the initial glucose concentration was set to 5.5 mM over the entire domain (initial condition), mimicking a physiological level of glucose in the ECM [34]; the concentration was fixed at the boundary of the domain (5.5 mM). As a result, the tumor spheroid placed in the center of the domain was subjected to an isotropic chemical stimulus, analogous to what was qualitatively reported in Figure 1a. In the gradient case, we defined different boundary conditions at the left and right edges of the domain, imposing fixed glucose concentrations of 8 [34,35] and 3 mM [36,37], respectively. As a consequence, the initial condition was a linear concentration gradient C ( x , y ) x = 2.5   m M / m , with an average concentration (and at the centre of the domain) still equal to 5.5 mM, as in the isotropic case. The initial concentration was constant along the y-axis.
The tumor spheroid was initially represented as an aggregate of a few cells, occupying a sub-domain with a radius of 50 μm (21 cells) in the center of the lattice representing the ECM, which, for the sake of simplicity, was assumed to be composed of collagen. Spheroid growth was simulated for 48 h (simulated time). During the experiment, the spheroid evolved while individual cells proliferated and migrated invading the ECM under the chemotactic stimuli of glucose concentration.

2.1.2. Glucose Layer

In the model, glucose was defined as a function of the spatial variable x = (x,y) and time t. The dynamics of the glucose concentration field C(x,t) in time and space is determined by solving the classical Fickian reaction–diffusion equation:
C x _ , t t = D 2 C x _ , t k x _ , t with   C x _ , 0 = C 0
where D = 7 10 10 m 2 / s [38] is the diffusion coefficient of glucose in the collagen-based ECM and k(x,t) is the consumption rate of glucose by the cells. The computational cost to obtain the numerical solution of Equation (1) can thus be high, given the considerable disparity between the time scales of cell division (hours to days) and glucose diffusion (seconds); cellular proliferation was treated as an adiabatic perturbation in the chemical field [29]. Thus, using the adiabatic perturbation approximation, Equation (1) is approximated as a pseudo-stationary problem [39]:
0 = D 2 C x _ , t k x _ , t
which is numerically solved by the method of simultaneous over-relaxation with the Chebyshev acceleration [29]. This method requires the discretization of the differential equation in terms of finite differences, where y and x represent the row and column indices of the elements on the grid, and Δ is the size of a single LP (i.e., each of the two edges along x and y):
( C x + 1 , y C x , y ) + ( C x 1 , y C x , y ) + ( C x , y + 1 C x , y ) + ( C x , y 1 C x , y ) 2 k x , y D = 0
where Δ is the size of a single LP (in our case 20 μm). Equation (3) can be rearranged by defining its residual ξx,y as:
ξ x , y = C x + 1 , y + C x 1 , y + C x , y + 1 + C x , y 1 4 + 2 k x , y D C x , y
and calculating an approximate solution for the concentration field in the subsequent simulation time step t + dt as:
C x , y t + d t = C x , y t + ξ x , y t 4 + 2 k x , y D
Glucose consumption by cells (in LPs where cells are present) is calculated according to a Michaelis–Menten kinetics:
k x , y = V m a x C x , y K m + C x , y I x , y
where V m a x = 0.05   m o l / ( m 3 s ) [40] is the maximum consumption rate, K m = 2   m o l / m 3 [40] is the Michaelis–Menten constant, and I is a Boolean indicator that takes on the value of 1 or 0 for lattice points that may or may not contain a cell. Because, in our model, two different cell phenotypes (active and starved, see subsequent section) able to consume glucose can occupy LPs, the indicator I in Equation (6) is calculated as the sum of two Dirac delta functions, one for each cell phenotype:
I x , y = δ x x a c t , y y a c t + δ x x s t a r v , y y s t a r v
where (xact, yact) and (xstar, ystarv) are the coordinates of the LPs where active and starving cells are located. The use of the Dirac delta function allows for the representation of a point-like object, in this case, the presence of a cell in a specific LP.

2.1.3. Cellular Layer and Cell Phenotypes

In the cellular layer of our model, the dynamics of individual cells were governed by specific rules that governed cellular dynamics and phenotyping, including migration, proliferation, starvation, and death.
Cell time evolution was governed by the nutrient condition. In the cellular layer, three different cellular phenotypes were distinguished: active, starving (or quiescent), and necrotic cells. The occurrence of various phenotypes depended on the cellular microenvironment and affected the metabolic activity of the cells. Two critical glucose concentrations defined the threshold for starving and necrotic cells (Cstarv > Cnec). A cell was active if the local glucose concentration Cx,y (i.e., concentration in the LP occupied) was higher than both these predefined thresholds, Cx,y > Cstarv > Cnec. An active cell consumed glucose according to Equation (6) and its metabolism included the possibility to migrate and eventually duplicate upon the completion of its cell cycle.
If Cnec < Cx,y < Cstarv, the cells did not have enough nutrients to duplicate, and thus starved. The cells still consumed residual glucose and could migrate, looking for more nutrient-rich LPs. As time progressed, if Cx,y increased back to levels higher than Cstarv, a starved cell could become active again, unless it remained in the starved state for too long a period, going in apoptosis (see subsequent section).
If the nutrient concentration decreased further, Cx,y < Cnec, the cell underwent necrosis. This process is irreversible, independent of any future change in the glucose concentration, and from that time on, the necrotic cell cannot proliferate, migrate, or consume glucose.

2.1.4. Rules Governing Cellular Dynamics

The dynamic evolution of cells is highly dependent on nutrient availability, particularly glucose concentration. In response to nutrient availability, cells can experience different fates, including necrosis and apoptosis (cell death).
Each non-necrotic cell in the cellular layer can undergo migration or (if active) proliferation. The two mechanisms are regulated by two independent characteristic times, defined as the migration time Tm and duplication time Td.
(a)
Cell migration
As the simulation time progresses, at every time interval Tm, a cell can migrate. The probability of the cell to migrate is affected by cell–cell interactions if the cell under evaluation is attached to a cell cluster and not isolated (detached). If a migration event occurs, the cell occupies one of its eight neighboring LPs, provided the target location is empty (Figure 1b); alternatively, if migration does not occur, the cell remains in the same position until the subsequent time interval, Tm.
When a cell tries to migrate, it must overcome cell adhesion defined by a probability Padh, defined as:
P a d h = N c e l l a t t K + N c e l l a t t if   attached P a d h = 0 if   attached
where N c e l l a t t is the number of attached cells confined with the cell of interest ( 0 < N c e l l a t t < 8 ), while K  K R is a parameter quantifying the role of cell–cell adhesion. For K = 0, cell–cell adhesion is supposed to be strong enough to inhibit the possibility of any movement to the cells, while high values of K are related to a weak cell–cell interaction and the higher motility of the cells. When an attached cell attempts to migrate, a random number 0 < σ < 1 is generated and compared to P a d h . If σ < P a d h , the migration fails, until the subsequent time interval, Tm. If σ > P a d h , migration is allowed and a change in the cell adhesion state can be induced; in our model, for simplicity, we assume that cell detachment is irreversible, i.e., once detached, a cell is not allowed to re-attach. If the cell of interest is already detached P a d h = 0 , it shows the typical behavior of isolated cells and is expected to migrate every Tm.
If migration is allowed, the cell has to choose a motility direction, which is defined according to a biased random walk approach. Each of the nine possible positions P i ( i = 1 , 2 , , 9 ) that a cell can take (Figure 1b), including the position already occupied (P5), is associated with a numeric interval R i = a i 1 ,   a i . The 9 intervals have different sizes ( R i = a i a i 1 ) , are contiguous and non overlapping, and span the entire range 0–1 ( 0 < R i < 1 ,   i = 1 9 R i = 1 ) (Figure 1c). A further random number 0 λ 1 is generated and compared to the R i intervals; if λ R i , the cell migrates towards direction P i , if the corresponding LP is empty. In particular, if λ R 5 , i.e., the LP already occupied by the cell, the cell does not migrate. If the cell, upon arrival at its new location, is not contiguous with any attached cells, it undergoes a transition into the detached state, provided it is not already in that state.
To evaluate the range of R i and define the probability of migration in each direction, a score S i is evaluated for each of the 9 candidate positions as:
S i = e α A i
where α R is a parameter describing the glucose chemotactic sensitivity of the cell, and A i is the local glucose specific gradient, defined as A i = ( C i C 5 ) / C 5 . S i values are further normalized s ¯ i = S i j = 1 9 S j . Finally, the ranges R i = a i ,   a i + 1 are calculated as a i = j = 1 i 1 s ¯ j (Figure 1c).
(b)
Cell proliferation
As the simulation time progresses, at every time interval Td, active cells proliferate, generating a new daughter cell, identical to its progenitor. The new cell is randomly located in one of the free positions among the eight LPs surrounding the progenitor. If no empty spot is available, the new cell’s location is chosen by identifying the direction where the minimum number of cells separate the progenitor from the edge of the cluster. All the cells along the selected direction shift one position away from the progenitor cell, and the new-born cell occupies the vacancy. It is worth mentioning our model considered the memory of the history of each cell to consider the possible cell-dependent variables, such as the random mutations of cell parameters, even if, in this work, this feature was not used. In each time step, all the cells were considered the same, with the only differences being among the active, starved, and necrotic phenotypes, while cell dynamic evolution was dependent on the nutrient availability only. The code implemented in our model also allowed for more complex interactions. Another simplification of our model was that only active cells attached to the spheroid could duplicate, while detached cells were expected to enter irreversibly in a migration state, unless nutrient availability induced their deaths.
(c)
Cell death
In our model, if a cell remained in starved state for a time longer than the apoptosis time (Tapop), it could spontaneously die. This biological event, known as apoptosis [41], is a mechanism of defense of the cells to prevent the propagation of lesions to the future generation. In our model, apoptosis corresponded to the degradation of the cell, which was dissolved in the ECM, and left the previously occupied LP empty. It is worth mentioning that, in the case where the glucose concentration reduces further to values below Cnec, while the cell is already in a starved state, the cell does not enter a state of apoptosis, but evolves into a necrotic state, where it remains indefinitely. Necrotic cells, according to our model, do not proliferate nor migrate, but continue indefinitely to occupy the same LP. The flowchart of the whole cell dynamic algorithm is presented in Figure 2. The model implementation of cell dynamics is also summarized as a sequence of operations in the figure caption.

2.1.5. Analysis

In this study, we focused on analyzing the effects of four key input parameters on the dynamic evolution of cell spheroids in our model: T d , T m , K , and α . By investigating the influence of these parameters, we aimed to gain insights into the roles of cell adhesion, proliferation, and migration in tumor growth and invasion processes.
In Table 1, the output variables computed in our model are briefly summarized and they include the number of cells attached within the spheroid ( N c e l l c o r e ), the number of cells that migrated away from the spheroid ( N c e l l m i g ), and the total number of cells ( N c e l l T O T ). Additionally, the percentages of adhered and migrated cells relative to the total number were calculated as N % c o r e and N % m i g , respectively. The ratio of adhered to migrated cells,   ϕ , was also determined by calculating the ratio of N c e l l c o r e / N c e l l m i g , and considered as a measure of the invasiveness of the tumor. The lower the ϕ value, the higher the tendency of the cancer cells to invade and colonize the surrounding ECM.
Furthermore, the spatial domain where the spheroid was located was divided into four sectors according to the diagonals of the square domain and named as North (N), South (S), East (E), and West (W) (Figure S1). The number of migrated cells located in each of these four quadrants were counted and represented as N c e l l N , N c e l l S , N c e l l E , and N c e l l W . The corresponding percentages of migrated cells in each quadrant relative to the total number of migrated cells ( N c e l l m i g ) were also calculated.

2.2. Statistical Analysis

In this study, we employed statistical analysis techniques to achieve a greater understanding of our model by conducting both global and local sensitivity analysis. GSA enabled us to assess the collective impact of perturbing all input parameters simultaneously, providing insights into how variations in multiple parameters influenced the system as a whole, while LSA allowed us to investigate the effects of perturbing individual parameters, providing detailed information on the sensitivity of the system to specific parameter changes. This approach allowed us to identify the parameters that had a greater impact on specific aspect of the model, and to identify the controlling mechanism that drove the entire dynamic of the system. In our model, we investigated, in particular, the effect of varying the input parameters related to cell duplication time, cell migration, chemotactic motility, and cell–cell adhesion: Td, Tm, K, and α.

2.2.1. Average and Variance Convergence

Since stochastic events and the values of randomly generated numbers governed our HCA model, each simulation run was unique, even if the values of the parameters, initial conditions, and boundary conditions were kept the same. For this reason, the outcomes of one single simulation for a given set of the four input parameters (Td, Tm, K, and α) were not enough to guarantee the statistical significance of the results. Therefore, for each set of input parameters, the simulation was replicated n times, defining the actual outcomes as the arithmetic average of the n iteration. To identify a value of n reliable from the statistical point of view, a convergence analysis was performed to study the effect of the number of simulation replicates on the model outputs (see Supplementary Materials, Figure S2). Convergence was verified on 4 different sets of model parameters, for the gradient experimental condition, running up to 100 simulations for each set. For n ≅ 10, all the model outcomes became constant and independent on n (data reported in Figure S2) for each of the 4 sets of parameters investigated. Therefore, in the results presented below, each simulation was reiterated and mediated for n = 40.

2.2.2. Global and Sensitivity Analyses

To explore the impact of simultaneous parameter perturbations on the key model outputs, a GSA was conducted, according to the protocols reported in the literature [42,43,44]. In brief, Latin hypercube sampling (LHS) was used to sample the multidimensional parameter space generating a set of 500 combinations of the 4 key parameters of interest (α ∈ [0, 30], K ∈ [0, 4], Td ∈ [10, 60 h] [45,46], and Tm ∈ [10, 120 min] [47]).
As there is no literature available to guide the choice of ranges for α and K, the selection was based on informed judgment and the biological context. For the α parameter, which was an index of cell chemotaxis, the lower limit of 0 represented a situation where cells did not sense concentration gradients, while the positive range reflected a scenario where the cells were attracted to increasing glucose concentrations. The upper limit of 30 was chosen to provide a degree of certainty that most cells would localize in regions that were rich in nutrients. For the K parameter, which was an index of the cell–cell adhesion strength, the lower limit of 0 corresponded to the strongest cell adhesion that inhibited any cell movement, while the upper limit of 4 was arbitrarily chosen to represent an adhesion strength weak enough for all cells to detach from the spheroid in a time frame of 48 h.
For each of the 500 combinations of the parameters, a simulation of n = 10 iterations was run. A multivariate linear regression analysis (MLRA) was performed using the built-in MATLAB function mvregress on a whole set of samples to obtain the linear regression coefficients for each parameter. The MLRA technique is based on the idea to express the output of interest Y i as a linear function of the input parameters X j , as shown in Equation (9):
Y i = ε i + β 1 X 1 + β 2 X 2 + β 3 X 3 + = ε i + j β j X j
ε i is an error term and β j is the regression coefficient of the input parameter. In our case, there were four input parameters ,   X j (Td, Tm, K, and α), and nine outputs, Y i , which were N c e l l c o r e , N c e l l m i g , N c e l l E , N c e l l S , N c e l l W   N % c o r e , N % E , N % S , and N % W . This subset of outputs, Y i , was chosen to be linearly independent, as required by the MLRA. The entire procedure (LHS and MLRA) was repeated 10 times, on 10 different sets of 500 combinations of the input parameter X j value combinations. As result, four regression coefficient distributions were obtained. Each distribution (one for each input parameter, X j ) represented the measure of the sensitivity index (SI) of the respective parameter.
Finally, a one-way ANOVA test was performed, followed by Tukey’s test (using the built-in MATLAB functions anova1 and multcompare) on the SI distributions. This procedure allowed us to rank the parameters according to their relative significance in affecting the model outputs. The GSA was conducted for both the isotropic and gradient cases.

2.2.3. Local Sensitivity Analysis

Each input parameter was perturbed independently, while keeping the others constant at their respective baseline values. The parameter ranges chosen were the following: chemotactic index α [−30, 30] (basal value α = 6), cell–cell adhesion parameter K [0, 4] (basal value K = 0.2), doubling time Td [10, 60 h] (basal value Td = 18 h), and migration time Tm [10, 120 min] (basal value Tm = 30 min). For α, the range was set between −30 and 30 to account for a chemo-repellent effect. The baseline value of 6 was arbitrarily chosen based on trial and error in order to represent a moderate chemotactic force. For K, the range was still between 0 and 4, with a baseline value of 0.2. Furthermore, the baseline value was selected through trial and error to ensure a balanced level of adhesion strength. For each parameter, 150 values were sampled uniformly and randomly within the specified range. Also, in this case, the LSA was conducted for both the isotropic and gradient cases.

3. Results

In this section, we present the results of the numerical model simulations, which aim to study the evolution of a tumor spheroid under two different glucose concentration profiles, defined as isotopic and gradient. The results are organized into three main paragraphs. In the first paragraph, we present the results of the simulation of the tumor spheroid’s evolution assuming baseline values for the input parameters of the model. In the second paragraph, we report the result of a GSA used to identify the key parameters that have the most significant impact on the tumor’s evolution. Finally, in the third paragraph, we report the results of an LSA of the parameters to further understand the relationship between the parameters and tumor evolution.

3.1. Baseline Case

In this paragraph, we present the results of the simulation of the tumor spheroid’s evolution imposing baseline values onto the four input parameters (i.e., α = 6, K = 0.2, Td = 18 h, and Tm = 30 min) both under isotropic and gradient conditions.
Figure 3 displays snapshots of the temporal evolution of the cellular and glucose layers at five time points (0, 12, 24, 36, and 48 h) for both isotropic (Figure 3a,b) and gradient (Figure 3c,d) conditions. The cellular layer (Figure 3a,c) displays the extracellular matrix (ECM) in dark blue, with the spheroid in the center with attached and detached cells, respectively, reported in orange and light blue. The glucose layer (Figure 3b,d) displays the concentration field of glucose in a color scale, ranging from blue (low concentration) to red (high concentration).
The general behavior observed was that the main tumor mass grew over time and consumed nutrients, inducing the formation of concentration gradients from the surrounding environment to the cell-populated area. This chemical stress stimulated cell motion and led to the cells’ detachment and migration toward regions less populated where a higher concentration of nutrients was available, in agreement with the predictions of the diffusional instability theory [9], and its experimental verification [38]. The two conditions investigated presented relevant differences. Under isotropic conditions, the cells moved radially away from the core, without any apparent preferential direction (Figure 3a). Under gradient conditions, as a consequence of the anisotropy in the stimulus, the detached cells tended to travel preferentially toward the source of nutrients at the left edge of the domain (Figure 3c).
It is worth mentioning that, in the base value conditions we investigated, given the limited size of the spheroid and the short time frame simulated, the glucose concentration in the nutrient layer, determined by diffusion and consumption, always remained higher than C s t a r v = 0.1   m o l / m 3 . Therefore, all cells remained proliferative, and neither a necrotic core nor a starvation rim were formed in the cellular layer. When running a simulation on a larger domain (5 mm × 5 mm) and for a longer time (20 days), the appearance of necrotic cells was observed, as expected (the results are reported in the Supplementary Figure S3).
To quantitatively investigate the observed phenomenon and evaluate the numerical measures of tumor evolution over time, we calculated the total number of cells N c e l l T O T and the number of cells that adhered to the spheroid N c e l l c o r e as a measure of tumor growth. The ratio of N c e l l c o r e to cells that migrated away from the spheroid N c e l l m i g , defined as , allowed us to compare the proliferation and migration activity, and it was a measure of the invasiveness of the tumour. We independently counted the number of cells that migrated in the four sectors (Figure S1) N c e l l N ,   N c e l l S ,   N c e l l W , a n d   N c e l l E and the respective percentages N % N ,   N % S ,   N % W , a n d   N % E . The comparison of these four values was a measurement of anisotropy in cancer invasion and a directional response to chemotactic stimuli (diffusional instability).
Figure 4 presents the time-course plots of isotropic and gradient conditions (Figure 4a–d and e–g, in the top and bottom parts of the figure, respectively). The solid lines report the average output (n = 40) values and the ribbons represent the associated standard errors. The orange data report the number (n = 40) of cells attached to the spheroid N c e l l c o r e , while the light-blue data represent detached cells N c e l l m i g . Figure 4a,e report the total cell count over time, while the graphs in Figure 4b,f report the percentages of the core and migrating cells. Figure 4c,g shows the time-course plot of ratio ϕ.
The results show that the tumor core ( N c e l l c o r e ) has an initial rapid growth rate over time, with an exponential trend for the first 24–36 h; following this initial rapid expansion, N c e l l c o r e continues to grow almost linearly (Figure 4a,e). The number of cells detaching from the main body and migrating ( N c e l l m i g ) increases linearly over time, starting from the initial value of 0, for the entire temporal window investigated (Figure 4a,e). Accordingly, the percentage of migrated cells ( N % m i g ) increases while the percentage of adhering cells ( N % c o r e ) decreases as time progresses. After an initial period of instability, both the percentages of adherent and migrated cells reach a steady-state plateau, with the core cells ( N % c o r e ) representing ~80% of the total number of cells in the domain (Figure 4b,f). This steady state can be easily seen by the constant value of the ratio between the numbers of adherent and migrated cells (ϕ = 4) after 24 h (Figure 4c,g). This result is due to the fact that both the size of the spheroid and the number of migrating cells increase linearly over time, with different rates, for long periods of time. The system, in other terms, reaches a steady state where the generation of cells due to proliferation and the flow of cells leaving the spheroid are balanced.
The abovementioned behavior is comparable in terms of total cells in the core and migrating cells, among the two conditions investigated here, i.e., isotropic and gradient, while a relevant difference can be observed in the direction of migration of the invading cells.
Figure 4d,h report the percentages of migrated cells in the four different directions, defined according to the source of nutrients (see Figure S1). In isotropic conditions, the detached cells do not show any preferential migratory direction and they uniformly distribute into the four sectors of the domain ( N % N N % E N % S N % W = 25 % ) (Figure 4d). On the other hand, in the gradient case, a significant difference (p<0.005) in the percentage of detached cells that migrate in the four directions is observed. As shown in Figure 4h, the nutrient gradient induces the cells to move in a biased random walkway, such that a higher percentage of the population migrates towards the West (higher glucose concentration). At the end of the simulation (t = 48 h) under anisotropic conditions, 42 ± 12% of cells migrate towards the direction of an increasing concentration gradient, and only 13 ± 8% of the cells move against the gradient, heading towards the East. The remaining migrating cells distribute almost uniformly between the remaining North and South sectors, with 21 ± 8% and 24 ± 10% of cells present in each sector, respectively.

3.2. Global and Sensitivity Analysis (GSA)

To perform the analysis, we estimated the SI values of each input parameter ( X j ) and ranked them according to their importance.
In Figure 5, the regression coefficients obtained from the MLRA for the isotropic and gradient cases are plotted. Each bar represents the regression coefficient mean, which is a measure of the SI parameter. The higher the SI, the greater the sensitivity of the parameter in affecting a given model output. To determine the ranking of the relative sensitivity of the model parameters, we conducted a one-way ANOVA test and post hoc analysis (Tukey’s test) after checking that the regression coefficients followed a normal distribution, a crucial requirement for the two tests. The bars marked with an asterisk indicated the statistically significant parameters for each model output, obtained from the MLRA.
From the GSA, we observed that N c e l l c o r e was independent on α (an analogous result can be obtained by measuring the core area), but it was significantly affected by T d   (Figure 5a,j). The number of migrated cells N c e l l m i g independent on α and T m under both isotropic and gradient conditions (see Figure 5b,k). Looking at the migration direction-related outputs N c e l l E , N c e l l S , N c e l l W , N % E , N % S , and N % W , relevant differences were visible when comparing the two conditions investigated. In the gradient case, all the migration direction-related outputs depended on all four parameters and, in particular, on α (Figure 5l–r). This was expected since α was related to chemotaxis. On the other hand, in the isotropic case, N c e l l E , N c e l l S , and N c e l l W were only dependent on Td and K (Figure 5c–e), with the former being the most relevant, while N % E , N % S , and N % W were only affected by T m (Figure 5g–i). From these observations, we can see that α only influences the outputs related to the directionality of cell migration, solely in the presence of a gradient, while it has no effect in isotropic conditions; parameters K and Td affect both the size of the spheroid and the number of cells detaching from it, in both isotropic and gradient conditions; and Tm has a reduced effect on spheroid growth compared to K and Td, but it significantly influences the directionality of cell migration, only in the presence of a gradient.

3.3. Local Sensitivity Analysis (LSA)

In order to define the empirical functional relationship between the parameters and model outputs, we conducted an LSA. The results of the LSA provide a greater insight into the non-linear relationships between the parameters and tumor evolution, and highlight the regions of the parameter space where the response variables are most sensitive to changes.
We analyzed the effects of the key input parameters α , K , T d , and T m on all model outputs (see Table 1). However, for the sake of brevity, we limited our discussion to only the most significant input–output relationships suggested by the results we obtained from the GSA.

3.3.1. Chemotaxis Sensitivity Index α

As described in the Materials and Methods Section, parameter α was an indicator of cell sensitivity to chemical stimuli, represented by glucose in this work. The higher the absolute value of α , the greater the tendency of the cells to migrate following a given concentration gradient.
Figure 6a,b depict the fraction of migrating cells invading from the spheroid along the four sectors (North, South, West, and East; see Figure S1) of the domain, which are subjected to different chemical stimuli, due to the geometry of our domain. The cells were counted after 48 h of a real-time simulation and reported on the y-axis as the fraction of cells in each direction; the x-axis showed the α values. The solid lines represent the mean fraction of cells, while the ribbons represent the associated standard deviations calculated for the n = 40 iterations (see Section 2.2.1).
Under isotropic conditions (Figure 6a), the fraction of cells migrating in the four sectors was independent on the α parameter, with the cells uniformly distributed throughout the entire domain ( N % N N % E N % S N % W 25 % ) in agreement with the isotropic stimulus imposed.
In contrast, under anisotropic conditions (Figure 6b), the distribution of migrating cells in the four sectors strongly depended on α. We also investigated the negative values of α to study the potential impact of a chemorepellent (such as a catabolite or toxic drug) on the migration of cancerous cells. Three cases are distinguishable in the graph: α = 0 , α > 0 , and α < 0 .
For α = 0 , the cells were insensitive to the concentration gradient and the system restored a pseudo-isotropic condition, with migrating cells uniformly distributed ( N % N N % E N % S N % W 25 % ) , as in the case of Figure 6a.
For α > 0 , the cells showed a preferential migration towards the source of the chemoattractant, which was West in our domain. As α increased, N % W (purple) increased linearly to reach saturation levels for high values of α ( α > 20 ), with the fraction of cells plateauing at around 85%. The fraction of cells in the remaining sectors decreased accordingly. It can be seen that the N % N (light blue) and N % S (orange) curves almost overlap, due to the North/South symmetry of our setup, while the East direction ( N % E , yellow) decreases more rapidly, being opposite to the source of chemoattractant.
A completely inverted situation was observed for α < 0 , where the cells moved away from the source of the chemical and showed preferential migration towards decreasing gradients, with the directionalities of the cells inverted. The trends of the curves associated with cells migrating towards the west and East were inverted, but the overall trend of the chart was symmetrical.

3.3.2. Cell–Cell Adhesion Parameter K

As described in the Materials and Methods Section, parameter K models cell–cell adhesion; the higher the value of K , the greater the tendency of an attached cell to move to one of the adjacent positions. This implies a greater chance that a cell can detach from the spheroid and invade the surrounding environment. Low values of K are related to a limited tendency to invade.
As described in the GSA, parameter K mostly affected the number of cells connected within the spheroid ( N c e l l c o r e ), the number of detached cells that migrated away ( N c e l l m i g ), and their percentages ( N % c o r e and N % m i g ). The general result, as expected, was that N % c o r e monotonically decreased as time progressed and K increased.
The trends of N c e l l c o r e and N c e l l m i g were evaluated at four different time points (12, 24, 36, and 48 h) as a function of parameter K , and are shown in Figure 6c,d, respectively. We reported the data only for the gradient case; the results obtained were independent of the concentration field, as confirmed by the GSA. As in the previous paragraph, solid lines represent the mean values of the cell numbers, while ribbons represent the associated standard deviations.
For simplicity, we distinguished four scenarios: K = 0 , 0 < K < 0.5 , K = 0.5 , and K > 0.5 . When K = 0 , cell adhesion is indissoluble and, thus, N c e l l c o r e monotonically increases over time, while N c e l l m i g remains at 0. For values of 0 < K < 0.5 , cell adhesion is still quite strong and N c e l l c o r e can still increase over time, although at a decreasing rate as K increases. This growth is slowed down by the progressive detachment of migratory cells, which increases over time and as K increases.
For the critical value of K = 0.5 , a steady-state equilibrium was established (Figure 6c, inset). The number of cells adhering to the spheroid remained constant over time and was approximately equal to the initial number of cells (21 cells, t = 0). On the other hand, the number of migratory cells evaluated at 48 h (red) reached its absolute maximum for this value of K. This indicates that there is a balance between the “flow” of cells leaving the tumor and the generation of new cells within the tumor.
For values of K > 0.5 , cell adhesion weakened further, and the generation of new cells could not compensate for the outflow of migratory cells. The number of adhering cells monotonically decreased over time until the spheroid completely disintegrated at 48 h. Similarly, N c e l l m i g decreased from the peak reached for K = 0.5 until it reached a plateau, as there was no longer a tumor core able to generate new cells (Figure 6d). The higher the value of K (above 1), the sooner the tumor core disappears.

3.3.3. Doubling Time Td

Based on the GSA results, Td mostly impacts the number of adherents ( N c e l l c o r e ) and migrating cells ( N c e l l m i g ) . Figure 6e,f display the trends of N c e l l c o r e and N c e l l m i g , evaluated at four different time points (12, 24, 36, and 48 h) as Td variates. Again, we present the results only for the gradient case. Upon increasing T d , the N c e l l c o r e count (Figure 6e) rapidly decreases, and a similar, but less marked, trend is followed by N c e l l m i g (Figure 6f), reconfirming the results of the GSA.

3.3.4. Migration Time Tm

Once a cell is detached, T m represents the time required by the cell to travel a distance equal to its diameter, thereby indicating the motility of cells. Alternatively, when a cell is attached, its capacity to move is limited, and T m represents the inverse of the frequency at which the attached cell attempts to move.
As described in the GSA, parameter T m mostly affects the fraction of cells comprising the spheroid ( N % c o r e ) both in isotropic and anisotropic cases, and only in the gradient case is the fraction of migrating cells located in the four sectors ( N % N ,   N % E ,   N % S ,   N % W ) .
Figure 6g depicts N % c o r e evaluated at four time points (12, 24, 36, and 48 h) in gradient conditions. Decreasing T m , the core reduces in size because a larger number of cells migrate toward the surroundings (Figure 6g) due to an increased tendency of migration.
Figure 6h shows N % N (light blue), N % S (orange), N % W (purple), and N % E (yellow) evaluated at 48 h as a function of Tm. In the gradient case, the effect of   T m is observed by an increase in the cell population percentage migrating West (a higher glucose concentration) when   T m decreases (Figure 6h). However, when   T m becomes sufficiently large, the four percentages ( N % N , N % S , N % E , N % W ) seem to converge to a ~25% value, which can be attributed to the reduced frequency of all cells to migrate.

4. Discussion

The HCA model developed in this study successfully simulated the temporal and spatial evolutions of a tumor spheroid under different chemical stimuli. The model incorporated the key biological mechanisms driving tumor progression, such as cell proliferation, migration, apoptosis, and interactions, with the extracellular matrix and neighboring cells. In this study, the growth and invasiveness of an avascular tumor spheroid were simulated in two different chemical fields, isotropic and anisotropic, for 48 h to study its evolution over time. For a fixed set of input parameters, and independent from the type of chemical field, the spheroid core area increased over time, and a fraction of cells detached and invaded the tumor’s surroundings. The effect of an external glucose gradient was to convert cell migration from a random walk to a chemotaxis-driven biased random walk.
Through GSA and LSA, we were able to understand how the selected input parameters affected the system, including how the parameters interacted with each other and how they acted individually. The GSA revealed that the impact of α was negligible in isotropic conditions but was a fundamental factor in governing tumor invasion under the gradient case; parameters Td and K were critical for the evolution of the core area and the number of cells detaching from the spheroid core, in both chemical fields, while Tm had a similar but less effective influence on tumor growth and invasiveness.
The LSA quantitatively assessed the effect of parameter perturbation, showing that the value of α influenced the migration direction of cells towards or away from a positive gradient, while α equal to zero resulted in cells becoming insensitive to the gradient. The parameter K, representing cell–cell adhesion, played a crucial role in tumor evolution, with K = 0.5 marking the border between scenarios where adhesion promoted proliferation, or the flux of migrating cells dissolved the tumor over time. Perturbing Td led to a drastic decrease in the numbers of core and invading cells, while perturbing Tm had the opposite effect, accelerating tumor growth.
Despite its simplicity, the HCA model proved to be highly effective in simulating tumor growth and invasion over time, capturing essential features of tumor biology and reproducing realistic spatiotemporal patterns under different experimental conditions. The model’s versatility allowed for easy modification and adaptation to include additional biological, chemical, and physical processes relevant to tumor growth and invasion. The model could be calibrated and validated using the experimental data obtained from in vitro or in vivo assays, thus providing a more accurate representation of the biological system under investigation.
Moreover, the model could be used to design and predict the outcomes of new experiment setups, thereby reducing the need for extensive and costly experiments. By incorporating new biological parameters extracted from simple and rapid experiments, the model could simulate more complex scenarios beyond the limitations of laboratory-based assays. Future improvements can involve the integration of complex signaling pathways, immune responses, hypoxia, nutrient transport, and mechanical stresses, such as oncotic pressure or laminar flows affecting the tumor surface.
In conclusion, the HCA model represents a promising tool for investigating the mechanisms of tumor growth and invasion, as well as for guiding the design and assessment of novel therapeutic strategies. Its ability to capture the complexity of tumor biology and its adaptability to various experimental settings makes it a valuable asset for cancer research.

5. Conclusions

This work introduced a model based on the HCA (hybrid cellular automaton) approach, capable of simulating tumor growth and cancer cell invasiveness. The model’s performance was tested, simulating the evolution of a 3D tumoral model, spheroid, at the early stages of its development, and under isotropic and anisotropic glucose concentration fields. Furthermore, this study investigated how perturbations in selected cellular parameters, including chemotactic sensitivity, cellular adhesion, doubling time, and cell motility, affected tumor dynamics.
These parameters were ranked according to their influence on model outcomes by performing a GSA. The results demonstrate that the automata-based model accurately describes the initial exponential-like growth of tumors and depicts cells adapting their migration mechanisms in response to external chemical stimuli. The chemotactic index was identified as crucial for cellular migration in the presence of a chemoattractant gradient, but it was irrelevant for the tumor development in isotropic chemical fields. Cellular adhesion plays an essential role in tumor growth and metastasis. The study revealed an optimal value of parameter K, representing cellular adhesion, where the tumor produced the highest number of migrating cells while maintaining a constant size over time. Additionally, the investigation highlighted that a hypothetical cancer cell line characterized by high motility and a short doubling time generated cancers with rapid growth outcomes and high invasiveness.
To ascertain the reliability and utility of our model, its validation through in vitro experiments was essential. For this reason, future developments of this research should primarily focus on validating the model with the empirical data derived from in vitro tumor spheroid growth studies. This crucial step is intended to refine the model’s practicality and expand our comprehension of tumor dynamics. Specifically, in our forthcoming studies, we will analyze the morphological responses of spheroids, focusing on the changes in area and cell detachment when subjected to chemical stimuli. This will include an in-depth investigation of various cell lines, each with unique attributes, such as doubling times, motility values, and cell–cell adhesion forces. These experimental validations will utilize the chemotaxis experiments previously conducted by our research group [38], providing a standardized and replicable framework for our investigations. Supporting this, the preliminary data presented in the Supplementary Video S1 demonstrate a promising qualitative concordance between the in vitro and in silico outcomes. This alignment pertains to both the morphological response and cell behavior over time, providing a qualitative comparison of the response to the chemoattractant. These experimental validations are anticipated to enhance the model’s accuracy and its capacity to elucidate complex aspects of tumor behavior.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/cancers15235660/s1, Figure S1: Schematic illustration of the division into four sectors of the cellular layer; Figure S2: Convergence analysis; Figure S3: Simulating advanced tumor morphology and temporal changes under isotropic conditions. Video S1: A qualitative comparison of the morphological response to the chemoattractant in both in silico and in vitro experiments, left and right, respectively.

Author Contributions

Conceptualization, S.C. and P.D.; methodology, L.M., M.J.P. and P.D.; formal analysis, P.D.; investigation, L.M.; data curation, R.F. and L.M.; writing—original draft preparation, L.M. and R.F.; writing—review and editing, P.D., Z.W., V.C. and S.C.; supervision, S.C. and P.D. All authors have read and agreed to the published version of the manuscript.

Funding

This work was partially supported by funding from the Cockrell Foundation (PD). This study was conducted under the umbrella of the International Academic Affiliation Agreement between Houston Methodist Academic Institute (Houston, TX, USA) and the University of Naples Federico II (Napoli, Italy).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Sung, H.; Ferlay, J.; Siegel, R.L.; Laversanne, M.; Soerjomataram, I.; Jemal, A.; Bray, F. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J. Clin. 2021, 71, 209–249. [Google Scholar] [CrossRef]
  2. World Health Organization. Global Health Estimates 2020: Deaths by Cause, Age, Sex, by Country and by Region, 2000–2019; World Health Organization: Geneva, Switzerland, 2020. [Google Scholar]
  3. Siegel, R.L.; Miller, K.D.; Fuchs, H.E.; Jemal, A. Cancer statistics, 2022. CA Cancer J. Clin. 2022, 72, 7–33. [Google Scholar] [CrossRef]
  4. Dillekås, H.; Rogers, M.S.; Straume, O. Are 90% of deaths from cancer caused by metastases? Cancer Med. 2019, 8, 5574–5576. [Google Scholar] [CrossRef]
  5. Fares, J.; Fares, M.Y.; Khachfe, H.H.; Salhab, H.A.; Fares, Y. Molecular principles of metastasis: A hallmark of cancer revisited. Signal Transduct. Target. Ther. 2020, 5, 28. [Google Scholar] [CrossRef]
  6. Roussos, E.T.; Condeelis, J.S.; Patsialou, A. Chemotaxis in cancer. Nat. Rev. Cancer 2011, 11, 573–587. [Google Scholar] [CrossRef]
  7. Vasaturo, A.; Caserta, S.; Russo, I.; Preziosi, V.; Ciacci, C.; Guido, S. A novel chemotaxis assay in 3-D collagen gels by time-lapse microscopy. PLoS ONE 2012, 7, e52251. [Google Scholar] [CrossRef]
  8. Caserta, S.; Campello, S.; Tomaiuolo, G.; Sabetta, L.; Guido, S. A methodology to study chemotaxis in 3-D collagen gels. AIChE J. 2013, 59, 4025–4035. [Google Scholar] [CrossRef]
  9. Cristini, V.; Frieboes, H.B.; Gatenby, R.; Caserta, S.; Ferrari, M.; Sinek, J. Morphologic instability and cancer invasion. Clin. Cancer Res. 2005, 11, 6772–6779. [Google Scholar] [CrossRef]
  10. Bearer, E.L.; Lowengrub, J.S.; Frieboes, H.B.; Chuang, Y.-L.; Jin, F.; Wise, S.M.; Ferrari, M.; Agus, D.B.; Cristini, V. Multiparameter computational modeling of tumor invasion. Cancer Res. 2009, 69, 4493–4501. [Google Scholar] [CrossRef]
  11. Cristini, V.; Lowengrub, J. Multiscale Modeling of Cancer: An Integrated Experimental and Mathematical Modeling Approach; Cambridge University Press: Cambridge, UK, 2010. [Google Scholar]
  12. Sixt, M.; Lämmermann, T. In vitro analysis of chemotactic leukocyte migration in 3D environments. In Cell Migration; Springer: Berlin/Heidelberg, Germany, 2011; pp. 149–165. [Google Scholar]
  13. Lämmermann, T.; Bader, B.L.; Monkley, S.J.; Worbs, T.; Wedlich-Söldner, R.; Hirsch, K.; Keller, M.; Förster, R.; Critchley, D.R.; Fässler, R. Rapid leukocyte migration by integrin-independent flowing and squeezing. Nature 2008, 453, 51–55. [Google Scholar] [CrossRef]
  14. Friedl, P.; Bröcker, E.B. The biology of cell locomotion within three-dimensional extracellular matrix. Cell. Mol. Life Sci. CMLS 2000, 57, 41–64. [Google Scholar] [CrossRef] [PubMed]
  15. Sant, S.; Johnston, P.A. The production of 3D tumor spheroids for cancer drug discovery. Drug Discov. Today Technol. 2017, 23, 27–36. [Google Scholar] [CrossRef] [PubMed]
  16. Ng, K.W.; Leong, D.T.W.; Hutmacher, D.W. The challenge to measure cell proliferation in two and three dimensions. Tissue Eng. 2005, 11, 182–191. [Google Scholar] [CrossRef]
  17. Rhodes, N.P.; Srivastava, J.K.; Smith, R.F.; Longinotti, C. Metabolic and histological analysis of mesenchymal stem cells grown in 3-D hyaluronan-based scaffolds. J. Mater. Sci. Mater. Med. 2004, 15, 391–395. [Google Scholar] [CrossRef]
  18. Nyga, A.; Cheema, U.; Loizidou, M. 3D tumour models: Novel in vitro approaches to cancer studies. J. Cell Commun. Signal. 2011, 5, 239–248. [Google Scholar] [CrossRef] [PubMed]
  19. Yamada, K.M.; Cukierman, E. Modeling tissue morphogenesis and cancer in 3D. Cell 2007, 130, 601–610. [Google Scholar] [CrossRef] [PubMed]
  20. Hoarau-Véchot, J.; Rafii, A.; Touboul, C.; Pasquier, J. Halfway between 2D and animal models: Are 3D cultures the ideal tool to study cancer-microenvironment interactions? Int. J. Mol. Sci. 2018, 19, 181. [Google Scholar] [CrossRef]
  21. Li, C.K.N. The glucose distribution in 9L rat brain multicell tumor spheroids and its effect on cell necrosis. Cancer 1982, 50, 2066–2073. [Google Scholar] [CrossRef]
  22. Kam, Y.; Rejniak, K.A.; Anderson, A.R.A. Cellular modeling of cancer invasion: Integration of in silico and in vitro approaches. J. Cell. Physiol. 2012, 227, 431–438. [Google Scholar] [CrossRef]
  23. Stillman, N.R.; Kovacevic, M.; Balaz, I.; Hauert, S. In silico modelling of cancer nanomedicine, across scales and transport barriers. NPJ Comput. Mater. 2020, 6, 92. [Google Scholar] [CrossRef]
  24. Karolak, A.; Rejniak, K.A. Micropharmacology: An in silico approach for assessing drug efficacy within a tumor tissue. Bull. Math. Biol. 2019, 81, 3623–3641. [Google Scholar] [CrossRef]
  25. Dogra, P.; Butner, J.D.; Chuang, Y.-l.; Caserta, S.; Goel, S.; Brinker, C.J.; Cristini, V.; Wang, Z. Mathematical modeling in cancer nanomedicine: A review. Biomed. Microdevices 2019, 21, 40. [Google Scholar] [CrossRef]
  26. Burdett, E.; Kasper, F.K.; Mikos, A.G.; Ludwig, J.A. Engineering tumors: A tissue engineering perspective in cancer biology. Tissue Eng. Part B Rev. 2010, 16, 351–359. [Google Scholar] [CrossRef] [PubMed]
  27. Norton, K.-A.; Gong, C.; Jamalian, S.; Popel, A.S. Multiscale agent-based and hybrid modeling of the tumor immune microenvironment. Processes 2019, 7, 37. [Google Scholar] [CrossRef] [PubMed]
  28. Kapałczyńska, M.; Kolenda, T.; Przybyła, W.; Zajączkowska, M.; Teresiak, A.; Filas, V.; Ibbs, M.; Bliźniak, R.; Łuczewski, Ł.; Lamperska, K. 2D and 3D cell cultures–a comparison of different types of cancer cell cultures. Arch. Med. Sci. AMS 2018, 14, 910. [Google Scholar] [CrossRef]
  29. Patel, A.A.; Gawlinski, E.T.; Lemieux, S.K.; Gatenby, R.A. A cellular automaton model of early tumor growth and invasion: The effects of native tissue vascularity and increased anaerobic tumor metabolism. J. Theor. Biol. 2001, 213, 315–331. [Google Scholar] [CrossRef] [PubMed]
  30. Wang, Z.; Zhang, L.; Sagotsky, J.; Deisboeck, T.S. Simulating non-small cell lung cancer with a multiscale agent-based model. Theor. Biol. Med. Model. 2007, 4, 50. [Google Scholar] [CrossRef]
  31. Costa, E.C.; Moreira, A.F.; de Melo-Diogo, D.; Gaspar, V.M.; Carvalho, M.P.; Correia, I.J. 3D tumor spheroids: An overview on the tools and techniques used for their analysis. Biotechnol. Adv. 2016, 34, 1427–1441. [Google Scholar] [CrossRef]
  32. Chang, S.L.; Cavnar, S.P.; Takayama, S.; Luker, G.D.; Linderman, J.J. Cell, isoform, and environment factors shape gradients and modulate chemotaxis. PLoS ONE 2015, 10, e0123450. [Google Scholar] [CrossRef]
  33. Gatenby, R.A.; Gawlinski, E.T. A reaction-diffusion model of cancer invasion. Cancer Res. 1996, 56, 5745–5753. [Google Scholar]
  34. Masur, K.; Vetter, C.; Hinz, A.; Tomas, N.; Henrich, H.; Niggemann, B.; Zänker, K.S. Diabetogenic glucose and insulin concentrations modulate transcriptom and protein levels involved in tumour cell migration, adhesion and proliferation. Br. J. Cancer 2011, 104, 345–352. [Google Scholar] [CrossRef] [PubMed]
  35. Lamers, M.L.; Almeida, M.E.S.; Vicente-Manzanares, M.; Horwitz, A.F.; Santos, M.F. High glucose-mediated oxidative stress impairs cell migration. PLoS ONE 2011, 6, e22865. [Google Scholar] [CrossRef] [PubMed]
  36. Palorini, R.; Votta, G.; Pirola, Y.; De Vitto, H.; De Palma, S.; Airoldi, C.; Vasso, M.; Ricciardiello, F.; Lombardi, P.P.; Cirulli, C. Protein kinase A activation promotes cancer cell resistance to glucose starvation and anoikis. PLoS Genet. 2016, 12, e1005931. [Google Scholar] [CrossRef] [PubMed]
  37. Leithner, K.; Hrzenjak, A.; Trötzmüller, M.; Moustafa, T.; Köfeler, H.C.; Wohlkoenig, C.; Stacher, E.; Lindenmann, J.; Harris, A.L.; Olschewski, A. PCK2 activation mediates an adaptive response to glucose depletion in lung cancer. Oncogene 2015, 34, 1044–1050. [Google Scholar] [CrossRef] [PubMed]
  38. Ferraro, R.; Ascione, F.; Dogra, P.; Cristini, V.; Guido, S.; Caserta, S. Diffusion-induced anisotropic cancer invasion: A novel experimental method based on tumor spheroids. AIChE J. 2022, 68, e17678. [Google Scholar] [CrossRef]
  39. Levenspiel, O. Chemical Reaction Engineering; John Wiley & Sons: Hoboken, NJ, USA, 1998. [Google Scholar]
  40. Schornack, P.A.; Gillies, R.J. Contributions of cell metabolism and H+ diffusion to the acidic pH of tumors. Neoplasia 2003, 5, 135–145. [Google Scholar] [CrossRef]
  41. Morana, O.; Wood, W.; Gregory, C.D. The apoptosis paradox in cancer. Int. J. Mol. Sci. 2022, 23, 1328. [Google Scholar] [CrossRef]
  42. Dogra, P.; Butner, J.D.; Ramírez, J.R.; Chuang, Y.-l.; Noureddine, A.; Brinker, C.J.; Cristini, V.; Wang, Z. A mathematical model to predict nanomedicine pharmacokinetics and tumor delivery. Comput. Struct. Biotechnol. J. 2020, 18, 518–531. [Google Scholar] [CrossRef]
  43. Saltelli, A. Global sensitivity analysis: An introduction. In Proceedings of the 4th International Conference on Sensitivity Analysis of Model Output (SAMO 2004), Santa Fe, NM, USA, 8–11 March 2004. [Google Scholar]
  44. Wang, Z.; Deisboeck, T.S.; Cristini, V. Development of a sampling-based global sensitivity analysis workflow for multiscale computational cancer models. IET Syst. Biol. 2014, 8, 191–197. [Google Scholar] [CrossRef]
  45. Pollock, K.; Stroemer, P.; Patel, S.; Stevanato, L.; Hope, A.; Miljan, E.; Dong, Z.; Hodges, H.; Price, J.; Sinden, J.D. A conditionally immortal clonal stem cell line from human cortical neuroepithelium for the treatment of ischemic stroke. Exp. Neurol. 2006, 199, 143–155. [Google Scholar] [CrossRef]
  46. Herlyn, M.; Steplewski, Z.; Herlyn, D.; Koprowski, H. Colorectal carcinoma-specific antigen: Detection by means of monoclonal antibodies. Proc. Natl. Acad. Sci. USA 1979, 76, 1438–1442. [Google Scholar] [CrossRef] [PubMed]
  47. Milo, R.; Phillips, R. Cell Biology by the Numbers; Garland Science: New York, NY, USA, 2015. [Google Scholar]
Figure 1. Model schematic showing cellular and glucose field layers. (a) Typical initial configuration of the cellular layer (upper grid) with cell-free LPs (dark blue) surrounding cell-occupied LPs (orange). Not in scale. In the glucose layer (bottom grid) glucose concentration is calculated in a pseudo-stationary condition in each lattice point (LP), according to Equation (6). In the scheme, a typical concentration profile in a cell spheroid is reported using a color scale. (b) Eight possible migration directions for a representative cell (located in P5); (c) schematic representation of the migration direction probabilities (Ri values).
Figure 1. Model schematic showing cellular and glucose field layers. (a) Typical initial configuration of the cellular layer (upper grid) with cell-free LPs (dark blue) surrounding cell-occupied LPs (orange). Not in scale. In the glucose layer (bottom grid) glucose concentration is calculated in a pseudo-stationary condition in each lattice point (LP), according to Equation (6). In the scheme, a typical concentration profile in a cell spheroid is reported using a color scale. (b) Eight possible migration directions for a representative cell (located in P5); (c) schematic representation of the migration direction probabilities (Ri values).
Cancers 15 05660 g001
Figure 2. Cell dynamic simulation cycle. The CA model iteratively scans the cellular layer, where each LP represents a cell. If the LP is occupied by a cell, the glucose concentration Cx,y is sampled from the glucose layer; if Cx,y < Cnec, the cell is classified as necrotic; if Cnec < Cx,y < Cstarv, the cell is classified as starved; if Cx,y > Cstarv, the cell is classified as active. Starved cells after Tapop undergo apoptosis and release the occupied LP. Active cells after Td proliferate. Active cells after Tm can migrate, according to the migration rules.
Figure 2. Cell dynamic simulation cycle. The CA model iteratively scans the cellular layer, where each LP represents a cell. If the LP is occupied by a cell, the glucose concentration Cx,y is sampled from the glucose layer; if Cx,y < Cnec, the cell is classified as necrotic; if Cnec < Cx,y < Cstarv, the cell is classified as starved; if Cx,y > Cstarv, the cell is classified as active. Starved cells after Tapop undergo apoptosis and release the occupied LP. Active cells after Td proliferate. Active cells after Tm can migrate, according to the migration rules.
Cancers 15 05660 g002
Figure 3. Baseline numerical solution. Snapshots of cellular layer (a,c) and glucose layer (b,d) at 0, 12, 24, 36, and 48 h following seeding of cancerous cells in the center of the simulation domain at 0 h, under isotropic (a,b) and gradient (c,d) conditions. Cellular layers display cells attached to the core (orange) and cells detached from the core (light blue), which migrate to the collagen matrix (dark blue).
Figure 3. Baseline numerical solution. Snapshots of cellular layer (a,c) and glucose layer (b,d) at 0, 12, 24, 36, and 48 h following seeding of cancerous cells in the center of the simulation domain at 0 h, under isotropic (a,b) and gradient (c,d) conditions. Cellular layers display cells attached to the core (orange) and cells detached from the core (light blue), which migrate to the collagen matrix (dark blue).
Cancers 15 05660 g003
Figure 4. Quantification of baseline solution. Results of the simulation assuming basal conditions with input parameters α = 6; K = 0.2, Td = 18 h, Tm = 30 min in isotropic (ad) and gradient cases (eh); (a,e) core (orange) and migrated cell numbers over time (light blue); (b,f) core (orange) and migrated cell percentages (light blue); (c,g) migrated cell number ratio ϕ; (d,h) percentage of cells migrated in the four directions: North (light blue), East (yellow), South (orange), and West (purple) over time. The solid lines represent the means (n = 40) and the ribbons represents the standard errors.
Figure 4. Quantification of baseline solution. Results of the simulation assuming basal conditions with input parameters α = 6; K = 0.2, Td = 18 h, Tm = 30 min in isotropic (ad) and gradient cases (eh); (a,e) core (orange) and migrated cell numbers over time (light blue); (b,f) core (orange) and migrated cell percentages (light blue); (c,g) migrated cell number ratio ϕ; (d,h) percentage of cells migrated in the four directions: North (light blue), East (yellow), South (orange), and West (purple) over time. The solid lines represent the means (n = 40) and the ribbons represents the standard errors.
Cancers 15 05660 g004
Figure 5. MLRA regression coefficients (SI) in isotropic (ai) and gradient (jr) cases. The columns refer to the key parameters α , K , T d , and T m . The rows refer to the outputs of interest: (a,j) core cell number; (b,k) migrated cell number; (c,l) number of cells migrating East; (d,m) number of cells migrating South; (e,n) number of cells migrating West; (f,o) core cell percentage; (g,p) percentage of cells migrating East; (h,q) percentage of cells migrating South; (i,r) percentage of cells migrating West. The asterisk indicates the significance parameters (p < 0.05) for the related outputs.
Figure 5. MLRA regression coefficients (SI) in isotropic (ai) and gradient (jr) cases. The columns refer to the key parameters α , K , T d , and T m . The rows refer to the outputs of interest: (a,j) core cell number; (b,k) migrated cell number; (c,l) number of cells migrating East; (d,m) number of cells migrating South; (e,n) number of cells migrating West; (f,o) core cell percentage; (g,p) percentage of cells migrating East; (h,q) percentage of cells migrating South; (i,r) percentage of cells migrating West. The asterisk indicates the significance parameters (p < 0.05) for the related outputs.
Cancers 15 05660 g005
Figure 6. Most relevant LSA results according to the GSA. Results of LSA perturbing input parameter α: (a) percentages of cells migrating North (light blue), East (yellow), South (orange), and West (purple) at 48 h in the isotropic case; (b) percentages of cells migrating North (light blue), East (yellow), South (orange), and West (purple) at 48 h in the gradient case. Results of LSA perturbing input parameter K: (c) core cell numbers at 12 (green), 24 (blue), 36 (yellow), and 48 h (red) in gradient case; (d) migrating cell numbers at 12 (green), 24 (blue), 36 (yellow), and 48 h (red) in gradient case. Results of LSA perturbing input parameter Td: (e) core cell numbers at 12 (green), 24 (blue), 36 (yellow), and 48 h (red); (f) migrating cell numbers at 12 (green), 24 (blue), 36 (yellow), and 48 h (red). Results of LSA perturbing input parameter Tm: (g) core percentages at 12 (green), 24 (blue), 36 (yellow), and 48 h (red); (h) percentages of cells migrating North (light blue), East (yellow), South (orange), and West (purple) at 48 h. The solid lines represent the means (n = 40) and the ribbons represent the standard deviations.
Figure 6. Most relevant LSA results according to the GSA. Results of LSA perturbing input parameter α: (a) percentages of cells migrating North (light blue), East (yellow), South (orange), and West (purple) at 48 h in the isotropic case; (b) percentages of cells migrating North (light blue), East (yellow), South (orange), and West (purple) at 48 h in the gradient case. Results of LSA perturbing input parameter K: (c) core cell numbers at 12 (green), 24 (blue), 36 (yellow), and 48 h (red) in gradient case; (d) migrating cell numbers at 12 (green), 24 (blue), 36 (yellow), and 48 h (red) in gradient case. Results of LSA perturbing input parameter Td: (e) core cell numbers at 12 (green), 24 (blue), 36 (yellow), and 48 h (red); (f) migrating cell numbers at 12 (green), 24 (blue), 36 (yellow), and 48 h (red). Results of LSA perturbing input parameter Tm: (g) core percentages at 12 (green), 24 (blue), 36 (yellow), and 48 h (red); (h) percentages of cells migrating North (light blue), East (yellow), South (orange), and West (purple) at 48 h. The solid lines represent the means (n = 40) and the ribbons represent the standard deviations.
Cancers 15 05660 g006
Table 1. Model outputs.
Table 1. Model outputs.
Output VariableFormulaDescription
N c e l l c o r e Number of cells adhered to the spheroid
N c e l l m i g Number of cells that migrated away from the spheroid
N c e l l T O T N c e l l c o r e + N c e l l m i g Total number of cells
N % c o r e ( N c e l l c o r e / N c e l l T O T ) 100 Percentage of cells adhered to the spheroid
N % m i g ( N c e l l m i g / N c e l l T O T ) 100 Percentage of cells that migrated away from the spheroid
ϕ N c e l l c o r e / N c e l l m i g Ratio of adhered cells to migrated cells
N c e l l N Number of migrated cells located in the North quadrant
N c e l l S Number of migrated cells located in the South quadrant
N c e l l E Number of migrated cells located in the East quadrant
N c e l l W Number of migrated cells located in the West quadrant
N % N ( N c e l l N / N c e l l m i g ) 100 Percentage of migrated cells located in the North quadrant
N % S ( N c e l l S / N c e l l m i g ) 100 Percentage of migrated cells located in the South quadrant
N % E ( N c e l l E / N c e l l m i g ) 100 Percentage of migrated cells located in the East quadrant
N % W ( N c e l l W / N c e l l m i g ) 100 Percentage of migrated cells located in the West quadrant
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Messina, L.; Ferraro, R.; Peláez, M.J.; Wang, Z.; Cristini, V.; Dogra, P.; Caserta, S. Hybrid Cellular Automata Modeling Reveals the Effects of Glucose Gradients on Tumour Spheroid Growth. Cancers 2023, 15, 5660. https://doi.org/10.3390/cancers15235660

AMA Style

Messina L, Ferraro R, Peláez MJ, Wang Z, Cristini V, Dogra P, Caserta S. Hybrid Cellular Automata Modeling Reveals the Effects of Glucose Gradients on Tumour Spheroid Growth. Cancers. 2023; 15(23):5660. https://doi.org/10.3390/cancers15235660

Chicago/Turabian Style

Messina, Luca, Rosalia Ferraro, Maria J. Peláez, Zhihui Wang, Vittorio Cristini, Prashant Dogra, and Sergio Caserta. 2023. "Hybrid Cellular Automata Modeling Reveals the Effects of Glucose Gradients on Tumour Spheroid Growth" Cancers 15, no. 23: 5660. https://doi.org/10.3390/cancers15235660

APA Style

Messina, L., Ferraro, R., Peláez, M. J., Wang, Z., Cristini, V., Dogra, P., & Caserta, S. (2023). Hybrid Cellular Automata Modeling Reveals the Effects of Glucose Gradients on Tumour Spheroid Growth. Cancers, 15(23), 5660. https://doi.org/10.3390/cancers15235660

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop