Next Article in Journal
Enumerating Subtrees of Flower and Sunflower Networks
Next Article in Special Issue
A Generalised Time-Dependent Mathematical Formulation for Magnetoelectrically Coupled Soft Solids at Finite Strains
Previous Article in Journal
Singularities for Timelike Developable Surfaces in Minkowski 3-Space
Previous Article in Special Issue
A Low-Cost Algorithm for Uncertainty Quantification Simulations of Steady-State Flows: Application to Ocular Hemodynamics
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Mathematical Modeling of Spherical Shell-Type Pattern of Tumor Invasion

Department of Mechanical Engineering, University of Victoria, Victoria, BC V8P 5C2, Canada
*
Author to whom correspondence should be addressed.
Symmetry 2023, 15(2), 283; https://doi.org/10.3390/sym15020283
Submission received: 24 December 2022 / Revised: 11 January 2023 / Accepted: 16 January 2023 / Published: 19 January 2023
(This article belongs to the Special Issue Symmetry in Finite Element Modeling and Mechanics)

Abstract

:
Cancer cell migration, as the principal element of tumor invasion, involves different cellular mechanisms. Various modes of cell migration including single and collective motions contribute to the invasion patterns. The competition between adhesive cell–cell and cell–matrix forces is a key factor that determines such patterns. In this paper, we study a distinct shell-type mode of tumor invasion observed in brain and breast tumors. In this mode, cells at the outer layer of the tumor collectively move away from the core and form a shell-type shape. Both the core and the shell sustain a sharp interface between cells and the surrounding matrix. To model the preserved interface, we adopted a Cahn–Hilliard-type free energy relation with the contribution of the interfacial stress. This nonconvex form of free energy allows for cells to remain together and preserve the tumor core via adhesive cell–cell forces while separating the core from the surrounding matrix across a continuous sharp interface. In addition, the motion of the shell was modeled using the chemotactic migration of cells in response to the gradient of nutrients. The associated fluxes of cells were implemented in a general form of balance law. A non-Michaelis–Menten kinetics model was adopted for the proliferation rate of cells. The flux of nutrients was also modeled using a simple diffusion equation. The comparison between the model predictions and experimental observations indicates the ability of the model to manifest the salient features of the invasion pattern.

1. Introduction

Cancer cell migration is the primary component of tumor invasion and metastasis, which are the main cause of death in patients [1]. Cell migration starts when cells attain a malignant morphology enabling them to detach the tumor body and move along the extracellular matrix (ECM) [2]. Cancer cell migration in three-dimensional tissue involves different cellular mechanisms; hence, various invasion patterns have been observed [3,4]. In general, cancer cells show two main migration patterns in histological sections: single and collective. In single-cell migration, e.g., in breast cancer, cells break intercellular contacts and acquire mobility via either mesenchymal or amoeboid morphology. Details of this invasion mode, the corresponding cell morphology characteristics, and the cellular mechanisms were reviewed in [5]. In collective cell migration, cells maintain intercellular adhesion and form an invading front sheet of cells that extend into the tissue and may or may not maintain their contact with the trailing cells and the primary tumor [6,7]. This group of cohesive cells maintains polarity, and expresses higher level of integrins and matrix metalloproteases (MMPs). Experimental studies using microfluidic platforms showed that the plasticity of tumor invasion modes depends on the tissue context, biochemical signals, and mechanical properties [8,9]. Collective migration is evidenced by human cancer pathology as the primary pattern of invasion in most epithelial cancers [10]. The tumor microenvironment (TME) plays an important role in the invasion pattern [11]. For example, collagen IV stimulates collective migration, while fibronectin triggers single-cell migration [12,13]. Collective migration can lead to a distinct pattern of invasion in which a cluster of cells start to break away from the tumor and form a precursor layer of moving cells. For instance, Koh et al. [14] showed that the ECM microenvironment can significantly impact patient-derived glioblastoma mutiforme (pdGBM) invasion and dissemination, as shown in Figure 1a. Another study, by Ling et al. [15], showed that altering cell–cell adhesion and contact in collagen-embedded MCF10AT1 cell-derived spheroids cocultured with adipose stromal cells (ASCs) isolated from human breast adipose tissue can lead to a collective invasion, shown in Figure 1b. These studies suggested that collective cell migration may give rise to a distinctive pattern of invasion in various types of tumors.
The diversity of mechanisms involved in tumor invasion requires different areas of knowledge to contribute to the development of treatment strategies [16,17]. Mathematical modeling can provide both qualitative and quantitative information to better understand the complex behavior of biological tissues, particularly tumors [18,19,20]. That information helps in interpreting the complicated interactions between tumors and their microenvironment. Tumor growth and invasion were mathematically studied from different perspectives: continuum, discrete, and hybrid (continuum-discrete) [21]. Continuum models, which are typically built upon partial differential equations (PDEs) to determine the concentrations of cells and nutrients, showed acceptable agreement with experimental observations [22,23,24,25]. These models correlate the growth and invasion with kinetic interactions between the key components, such as the rates of nutrient/oxygen consumption, cell proliferation, and diffusion for transport phenomena [26,27]. The ability of tumors to maintain a compact aggregation of cells was modeled by introducing adhesive cell–cell forces [28]. Further studies investigated other factors, e.g., level of PH, to predict the growth and invasion of tumors [29,30]. Other aspects of tumor progression, such as anisotropic growth and mechanical stress distribution, were also investigated in the framework of continuum mechanics [31,32,33,34]. On the other hand, discrete mathematical modeling incorporates signaling pathways, and inter- and intracellular interactions in tumor progression. These models are capable of predicting tumor progression at the subcellular level. For instance, a lattice Monte Carlo model was utilized to predict the dynamics of avascular tumor growth, and the effects of microenvironmental conditions and growth promoters/inhibitors on tumor cell survival [35]. Another example of discrete models is cellular automata (CA), which predict the collective behavior of self-organizing systems, such as tumors, by focusing on the interactions between the components [36].
To relate the cellular mechanisms with tissue-level responses, hybrid mathematical models integrate continuum and discrete approaches, and allow for descriptions of discrete cellular interactions together with macroscopic variables [21]. The effects of cellular pathways and interactions on the clinical-size morphology of tumors can be studied using hybrid models [37].
The results of these studies provide the prediction of various growth and invasion patterns, including the standard dissemination of cells, anisotropic growth, and asymmetric invasion, but are missing the exposition of the spherical-shell-type invasion pattern observed in breast and brain tumors.
Ttumors exhibit different patterns of invasion, among which the spherical shell-type pattern is a distinct mode that is less studied. In this work, we study this unique mode of invasion in which cancer cells collectively migrate away from the tumor, forming the precursor layer of a moving shell. The interplay between cell–cell and cell–matrix adhesion allows for the shell and the tumor to develop and preserve sharp interfaces. To capture this, we adopted a nonconvex free energy that allowed for phase separation. The motion of the shell was also modeled using chemotactic migration in response to the nutrient gradient. The role of various parameters in the invasion pattern, such as cell and nutrient diffusion, cell proliferation, and nutrient uptake, were investigated. The model was able to predict the formation and movement of the spherical shell pattern. The novelty of the presented work is in incorporating the interplay between cell–cell and cell–matrix adhesion using a nonconvex free energy in which the capillary stress allowed for a continuous transition between low and high concentrations across the interface.

2. Model Formulation

Different forces may induce cell migration within a tumor microenvironment. A typical form of cell migration is due to the nonuniform distribution of cells within tumor aggregation, which governs uniform distribution via random motions. This form of cell migration is correlated to the gradient of concentration, and the corresponding cell flux may be obtained using a standard Fick’s law relationship. In contrast, adhesive forces due to cell–cell and/or cell–matrix bindings prevent cell diffusion, hence maintaining the structure of a tumor. Another stimulus that may regulate cell migration either collectively or individually is the gradient of chemoattractants such as nutrients and growth factors. In this form of motion, cells break their cell–cell bindings, e.g., E-cadherin junctions, change their morphology, and move in response to the gradients of chemoattractants, so-called chemotaxis [38,39]. The competition among diffusive, adhesive, and chemotactic forces determines the ultimate pattern of cell migration.
To model the contribution of different forces, we adopted a general form of balance equations for cells and nutrient concentrations,
C t + div J C = ζ ,
N t + div J N = η ,
where C and N are the concentrations of cells and nutrients, i.e., mass per volume (kg·m 3 ), J C and J N are the cell and nutrient fluxes (kg·m 2 · s 1 ), η is the rate of nutrient consumption (kg·m 3 · s 1 ), and ζ is the rate of cell proliferation (kg·m 3 · s 1 ). The flux of nutrients in the tumor microenvironment can be described with a standard Fick’s law equation. However, cell flux is mainly generated according to the interplay between diffusive and adhesive forces, and the directional motion due to chemotaxis. As observed in experimental studies, when cell–cell adhesive forces are dominant, tumors are able to maintain their compact structure, while highly invasive tumors could exhibit a long trace of cell migration [40]. In addition, the composition of the ECM surrounding the tumor is an important stimulus that may switch the invasion mode [14].
To account for the overall flux of cells, we considered a nonconvex form of free energy and a standard chemotactic model. The total flux of cells and nutrients can be expressed as follows:
J C = J f + J d ,
J N = D N grad N ,
where J f is the flux due to the free energy, J d is the chemotactic flux and D N is nutrient diffusivity coefficient. In the next section, we study both fluxes and present a full formulation for Equation (1).

2.1. Free Energy and Interfacial Flux

The pattern of invasion in Figure 1 shows that the tumor’s sharp interface with the surrounding gel was sustained. The continuous transition from high to low concentrations across the interface required a distributed surface tension. The surface tension of the diffusive layer between the regions with different concentrations generated a form of capillary stress. This implies that the properties must depend on the interface contribution, i.e., Ψ = ( grad C ¯ · grad C ¯ ) / 2 , as proposed by Anderson et al. [41]. Hence, the free energy could take the form f = f ^ ( C ¯ ) + f ˇ ( Ψ ) . We adopted a nonconvex form of free energy, similar to Cahn–Hilliard interfacial free energy [42], which allowed for phase separation:
f = C ¯ + C ¯ 2 / 10 C ¯ 3 / 6 + C ¯ 4 / 12 + β Ψ ,
in which C ¯ is the dimensionless cell concentration, and β is a coefficient related to the length of transition that should be determined on the basis of the nature of the interaction between the cells and the ECM. Other coefficients were selected such that the corresponding flux retained two nontrivial roots as the equilibrium states in the domain 0 C ¯ 1 . This form of free energy allows for cells to remain together and form the tumor body due to cell–cell adhesive forces with a continuous transition across the interface that could form a sharp edge and separate the tumor from the surrounding matrix. According to [41], the cell flux takes the following form:
J f = D f grad f C ¯ div f Ψ grad C ¯ ,
where D f is a diffusion coefficient. We considered the isothermal process and neglected the dependency on temperature. The substitution of (3) into (4) yields
J f = D f grad C ¯ 5 C ¯ 2 2 + C ¯ 3 3 β div grad C ¯ .
The proposed form of free energy allows for an equilibrium state of separated phases, as the term f ^ / C ¯ in (4) had two extrema, and the corresponding flux function J ^ = D f grad ( f ^ / C ¯ ) had two nontrivial roots, as depicted in Figure 2a.
Cell distribution at equilibrium in the absence of nutrients is shown in Figure 2b. The initial distribution of cells, that is, the initial tumor consisting of a concentrated core and a sharp interface with the matrix, reached an equilibrium state of high C h and low C l concentrations coexisting with a continuous transition. The values of C h and C l were determined on the basis of the coefficients in Equation (3).
Next, we describe the role of nutrients in cell proliferation and chemotactic flux.

2.2. Chemotactic Flux and Cell Proliferation

Cellular migration mechanisms depend on cell morphology and external stimuli. Invasive tumor cells become morphologically polarized and develop membrane protrusions, allowing for them to migrate forward, the so-called mesenchymal motion, when cells move via traction and adhesion to the ECM [43]. In addition, cells squeeze through the pores in the ECM, i.e., amoeboid motion, which is relatively fast and does not require strong cell–matrix adhesive forces [44]. Nutrients also play a significant role in inducing directional motion, i.e., chemotaxis [45]. Observations suggest that the gradient of nutrients is the characteristic property of nutrient distribution that can drive chemotaxis. Therefore, we adopted a simple form of the chemotactic flux:
J d = μ C grad N ,
where μ is the chemotactic coefficient, and the flux was towards larger values of nutrient concentration. In addition, the rate of cell proliferation depends on the local concentrations of nutrients, as they are the main components of cellular metabolism and essential for cell growth. Cells nonlinearly consume nutrients, i.e., the higher the concentration of the supplied nutrient is, the higher the consumption rate and hence the higher the rate of proliferation, which saturates at a certain nutrient level. In addition to the nutrients, cells require intercellular proteins and enzymes for their growth. The local concentrations of these components are proportional to the concentration of cells. The cellular metabolisms in which nutrients and intercellular components are involved take place through cooperative enzymatic reactions [46]. Therefore, we adopted a non-Michaelis–Menten kinetics model for the proliferation rate ζ that included two substrates, nutrient N and cells C concentrations:
ζ = V m a x 1 1 + ( k 1 / C ) c 1 1 1 + ( k 2 / N ) c 2 ,
where V m a x is the maximal (saturated) rate, k 1 , k 2 are reaction constants, and c 1 , and c 2 are the Hill coefficients. We took c 1 = c 2 = c for the sake of simplicity. The nutrient consumption rate is also concentration-dependent. For simplicity, we took the linear form η = η 0 C N , where η 0 is the consumption coefficient.

2.3. Problem Statement

In vitro studies of tumors often involve using microfluidic chips comprising a main chamber for the tumor-embedded hydrogel and side channels to supply the nutrients, as shown in Figure 3. In such experiments, tumors are either mixed with the hydrogel and injected into the chamber or placed within the hydrogel already filling the chamber. In either case, the initial condition of tumors is essentially a compact aggregation of cells preserved by adhesive cell–cell forces that start responding to the matrix and the nutrient when placed inside the chamber. To simulate such a condition, we set the initial condition of cells to be a uniform axisymmetric spherical distribution that accounted for the tumor, followed by a sharp interface with the surrounding matrix. For nutrients, we took uniform initial concentration. The radius of the chamber R is typically much larger than that of tumor R 0 ; hence, we considered zero flux at the boundary, shown by J c = 0 in the figure. Nutrients are constantly and symmetrically provided via the microfluidic channels; therefore, we set a constant concentration of nutrient supply at the boundary, N 0 = N s . The simulation was performed over a spherical domain, as shown in the figure.

2.4. Governing Equations

Balance Equations (1) can be expressed using dimensionless parameters. Scaling the distance with the domain radius R , time with total time τ , and cell and nutrient concentrations with initial cell concentration C 0 and the nutrient source concentration N s , respectively, yields
C ¯ t ¯ = α div grad f ¯ ( C ¯ ) β ¯ div ( grad C ¯ ) γ C ¯ grad N ¯ + V ¯ ( 1 + ( k ¯ 1 / C ¯ ) c ) ( 1 + ( k ¯ 2 / N ¯ ) c ) ,
N ¯ t ¯ = ψ div grad N ¯ η ¯ C ¯ N ¯ ,
with dimensionless parameters introduced as C ¯ = C / C 0 , N ¯ = N / N s , t ¯ = t / τ , r ¯ = r / R , α = τ D f / R 2 , β ¯ = β R 2 , γ = N s μ / D f , V ¯ = τ V m a x / C 0 , k ¯ 1 = k 1 / C 0 , k ¯ 2 = k 2 / N s , f ¯ = f / C 0 , ψ = τ D n / R 2 and η ¯ = τ η 0 C 0 . Using spherical coordinates and assuming axisymmetric distributions, cell and nutrient concentrations were only a function of t and r. The system of equations presented in (8) was subjected to the initial and boundary conditions followed by Section 2.3. We took the uniform initial concentrations of cells C 0 = const . within R and C 0 = 0 elsewhere, and nutrients N 0 = N s :
C ( r , t = 0 ) = C 0 , 0 < r < R 0 0 , elsewhere and N ( r , t = 0 ) = N s ,
where R 0 is the initial radius of the tumor. At the boundary, we took the cell flux to be zero and applied a fixed nutrient source concentration N s = const . ; hence, boundary conditions were
C ( r , t ) r | r = R = 0 ,
N ( r = R , t ) = N s .
The symmetry condition at the center of the tumor applied to both cells and nutrients:
C ( r , t ) r | r = 0 = N ( r , t ) r | r = 0 = 0 .
The nonlinear system of Equations (8) was numerically solved with the finite differences method. Considering radially symmetric solutions, a uniform spatial discretization of 1000 grid points was applied, and the Euler method was implemented for forward integration in time using 5000 time steps. The convergence of the solution was verified with a sequence of grid refinements. We took the following baseline values for the parameters, C 0 = 1 × 10 2 kg/cm 3 , N s = 1 × 10 3 kg/cm 3 , R = 1 × 10 1 cm, R 0 = 5 × 10 2 cm, D f = 1 × 10 8 cm 2 /s, μ = 1 × 10 7 cm 2 /s, D n = 1 × 10 6 cm 2 /s, k 1 = 5 × 10 3 , k 2 = 5 × 10 4 and τ = 24 h.
In the following section, we study the role of dimensionless parameters in tumor invasion pattern. First, we obtained the steady-state solution to Equation (8a) in the absence of nutrients. Next, the sustainability of the tumor’s sharp interface, the distribution of the migrated cells, and the shape and the formation of the shell are demonstrated in response to the variation of the dimensionless parameters β ¯ , γ , V ¯ , η ¯ , ψ and c by solving the full system of Equation (8a,b) subjected to (9)–(11).

3. Results

To better understand the role of free energy and the interfacial stress proposed in Section 2.1, we first look into the steady-state axisymmetric spherical form of Equation (8a):
· f ¯ ( C ¯ ) β ¯ · ( C ¯ ) = 0 ,
where ( · ) = ( · ) / r and · ( · ) = ( 2 / r ) ( · ) / r + 2 ( · ) / r 2 . This equation, in the absence of nutrients, represents the steady state cell distribution when the adhesive and diffusive forces are balanced. This equilibrium was determined with parameter β ¯ , as shown in Figure 4.
Dimensionless parameter β ¯ controls the transitional length between the high- and low-concentration regions. Higher values of β ¯ led to a smooth edge where the balance of diffusive and adhesive forces occurred at a longer transitional length. However, lower values of β ¯ allowed for the tumor to sustain a sharp continuous interface. The physical interpretation of β ¯ could be explained by considering the role of ECM in cell migration. An ECM composition that allows for strong cell–matrix adhesion and facilitates the mesenchymal motion corresponds to higher values of β ¯ . In this case, cells migrated through the matrix and exhibited a typical finger-type invasion pattern where the tumor interface was generally indistinct. However, lower values of β ¯ are associated with the ECM composition that regulates weak adhesive bindings with invasive cells. In this case, invasion is either inhibited or limited to the amoeboid motion of cells through the matrix pores. In addition, amoeboid-mesenchymal plasticity comprises both migration patterns that could be regulated by a specific ECM composition, and allows for a switch between the aforementioned modes.
Next, we study the role of dimensionless parameters γ , V ¯ , η ¯ , ψ , and c for the baseline values presented in Section 2.4, depicted in Figure 5.
Dimensionless cell diffusion number γ = N s μ / / D represents the ratio of chemotactic flux to the diffusive flux, for particular values of nutrient source concentration N s . Higher values of γ correspond to higher chemotactic effects that derive cells migrating toward the gradient of nutrients and initiate the formation of a peak in the cell concentration at the precursor cells, as shown in Figure 5a. Increasing γ gave rise to a higher level of chemotactic invasion with a less distinct interface between the tumor core and the matrix, while the leading cluster of cells still preserved the interface.
A similar effect related to V ¯ , i.e., dimensionless cell growth number, can be seen in Figure 5b, where a higher growth number V ¯ increased cell concentration and led to the formation of a peak in cell concentration. A comparison between the effects of γ and V ¯ indicates that a variation in γ only changed the concentration of cells outside the tumor core, while variation in V ¯ changed cell concentration everywhere. Dimensionless nutrient consumption number η ¯ had an opposite effect to that of growth number γ . Higher consumption led to a uniform decrease in cell concentration and eliminated the leading interface, as shown in Figure 5c.
Dimensionless nutrient diffusion number ψ plays a dual role in cell migration, as shown in Figure 5d. Higher values of ψ led to the formation of a more distinct tumor edge separating the tumor core and migrated cells, but removing the peak and the leading edge. Higher diffusion reduced the nutrient gradient; hence, less chemotactic migration occurred within the invasion zone.
Hill coefficient c determines the formation of a peak in the cell concentration at the leading edge. As depicted in Figure 5e, c alters the sensitivity of migrating cells to the change in cell and nutrient concentrations. Higher values of c reduce the concentration of cells in the invasion zone due to a significant change in the proliferation rate while leading to the formation of a notable peak at the leading interface. The peak observed in the axisymmetric radial solution corresponded a spherical-shell-type cluster of cells in the three-dimensional domain.
Lastly, the dynamic progression of the invasion and the nutrient distribution are shown in Figure 6. As captured by the model, a cluster of cells at the outer layer of the tumor broke away and formed a shell-type layer of cells moving up the gradient of nutrients while preserving a notable interface, as the experimental data in Figure 1 show. This can be explained as the interplay between adhesive cell–cell and cell–matrix forces that led to the separation of a cell cluster in the form of a shell layer diffusing into the surrounding matrix. As discussed in the previous section, we modeled this phenomenon by adopting nonconvex free energy and accounting for the interfacial properties, which allowed for the equilibrium states of high and low concentrations. We also adopted a simple reaction–diffusion model for the chemotactic migration of cells. The formation of the shell was predicted as a peak in the cell concentration over a small interval followed by a drop denoting the interface. This peak corresponded to the cell aggregation seen as a ring and marked by red arrows in Figure 1. As shown in the figure, the tumor and the invading cells preserved their interfaces with the surrounding gel while the invasion proceeded. The same behavior was predicted by the model, as shown in Figure 6, where the shell-type cluster of cells preserved a sharp continuous interface and broke away. The gradient of nutrients due to the consumption by the precursor cells evoked the chemotactic migratory response, hence inducing the radial motion of the shell over time.
The proposed model presents a mathematical framework that translates biological phenomena, such as cell–cell/–matrix adhesion, chemotaxis, proliferation and consumption rates, into the context of free energy and capillary stress. Using this model, we were able to interpret the interplay between various biological factors that influenced tumor invasion and their effects on inducing a distinctive pattern of invasion observed in different tumors. In this invasion mode, a cluster of cells migrated into the surrounding matrix, with a preserved interface exhibiting a spherical-shell-type layer moving away from the tumor core. This pattern may be perceived as a ring in experimental studies due to the projection of tumor images onto a two-dimensional plane in microscopic imaging. The results presented here are limited to the axisymmetric invasion of tumors with a homogeneous cell population. Anisotropic tumor invasion may occur due to different stimuli, such as heterogeneity in a cell population, e.g., hypoxic and normoxic cell phenotypes, and fiber-embedded tumor ECM.

4. Conclusions

In this work, we developed a continuum model to study the distinctive mode of tumor invasion resulting from the interplay between cell–cell and cell–matrix adhesion. A nonconvex free energy is integrated with a standard chemotactic model to predict the collective migration of cells in the form of a shell-type invasion. The adopted Cahn–Hilliard-type free energy with the contribution of the interfacial stress allowed for the cells to remain together and preserve the tumor core via adhesive cell–cell forces while separating the shell from the core and the surrounding matrix across continuous sharp interfaces. Cells moved in response to the gradient of nutrients, i.e., chemotaxis, while preserving the shell layer. This indicates that the adhesive cell–cell force was strong enough to maintain the structure of the shell throughout the migration. Due to nutrient consumption, the chemotactic effect was reduced as the shell passed along, leaving behind a region with fewer cells. The effects of different parameters, such as diffusion, growth, and consumption numbers, on cell distribution were studied. Higher values of cell diffusion and growth numbers increased the growth of cells moving outwards from the tumor core, yielding a peak in cell concentration with a distinct interface. Higher values of nutrient consumption and diffusion numbers reduced cell concentration and eliminated the peak. In addition, the formation of the peak was controlled by the Hill coefficient, c, which determined the sensitivity of the proliferation rate to the gradient of nutrients, giving rise to a larger peak at the leading interface. The proposed model could demonstrate the principal features of the invasion pattern, including the formation, shape, and movement of the shell, and the continuous interface separating the tumor from the surrounding matrix. Nevertheless, empirical studies are required to calibrate the model parameters for quantitative predictions. Future research directions include experimental studies on tumor invasion patterns using tumor-on-a-chip platforms that enable testing different hydrogel compositions, and help in understanding the role of cell/matrix adhesion. This further helps in the calibration of model parameters to enhance quantitative predictions. In addition, the model can be integrated with a discrete mathematical framework to include cellular interactions. This improves the accuracy of predictions in the case of anisotropic invasion.

Author Contributions

Conceptualization, M.A., H.S. and B.N.; methodology, M.A., H.S. and B.N.; formal analysis, M.A., H.S. and B.N.; software, M.A.; writing—original draft preparation, M.A.; writing—review and editing, M.A., H.S. and B.N.; supervision, B.N. and H.S.; funding acquisition, B.N. All authors have read and agreed to the published version of the manuscript.

Funding

Funding was provided by the Natural Science and Engineering Research Council (NSERC) of Canada.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data will be available upon request.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Van Zijl, F.; Krupitza, G.; Mikulits, W. Initial steps of metastasis: Cell invasion and endothelial transmigration. Mutat. Res. Mutat. Res. 2011, 728, 23–34. [Google Scholar] [CrossRef] [PubMed]
  2. Spano, D.; Heck, C.; De Antonellis, P.; Christofori, G.; Zollo, M. Molecular networks that regulate cancer metastasis. In Proceedings of the Seminars in Cancer Biology; Elsevier: Amsterdam, The Netherlands, 2012; Volume 22, pp. 234–249. [Google Scholar]
  3. Friedl, P. Prespecification and plasticity: Shifting mechanisms of cell migration. Curr. Opin. Cell Biol. 2004, 16, 14–23. [Google Scholar] [CrossRef]
  4. Liu, S.; Li, Q.; Chen, K.; Zhang, Q.; Li, G.; Zhuo, L.; Zhai, B.; Sui, X.; Hu, X.; Xie, T. The emerging molecular mechanism of m6A modulators in tumorigenesis and cancer progression. Biomed. Pharmacother. 2020, 127, 110098. [Google Scholar] [CrossRef]
  5. Yilmaz, M.; Christofori, G.; Lehembre, F. Distinct mechanisms of tumor invasion and metastasis. Trends Mol. Med. 2007, 13, 535–541. [Google Scholar] [CrossRef] [PubMed]
  6. Hegerfeldt, Y.; Tusch, M.; Bröcker, E.B.; Friedl, P. Collective cell movement in primary melanoma explants: Plasticity of cell-cell interaction, β1-integrin function, and migration strategies. Cancer Res. 2002, 62, 2125–2130. [Google Scholar] [PubMed]
  7. Weijer, C.J. Collective cell migration in development. J. Cell Sci. 2009, 122, 3215–3223. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  8. Yang, W.; Luo, C.; Lai, L.; Ouyang, Q. A novel microfluidic platform for studying mammalian cell chemotaxis in different oxygen environments under zero-flow conditions. Biomicrofluidics 2015, 9, 44121. [Google Scholar] [CrossRef] [Green Version]
  9. Wolf, K.; Wu, Y.I.; Liu, Y.; Geiger, J.; Tam, E.; Overall, C.; Stack, M.S.; Friedl, P. Multi-step pericellular proteolysis controls the transition from individual to collective cancer cell invasion. Nat. Cell Biol. 2007, 9, 893–904. [Google Scholar] [CrossRef]
  10. Bronsert, P.; Enderle-Ammour, K.; Bader, M.; Timme, S.; Kuehs, M.; Csanadi, A.; Kayser, G.; Kohler, I.; Bausch, D.; Hoeppner, J.; et al. Cancer cell invasion and EMT marker expression: A three-dimensional study of the human cancer–host interface. J. Pathol. 2014, 234, 410–422. [Google Scholar] [CrossRef]
  11. Connor, Y.; Tekleab, Y.; Tekleab, S.; Nandakumar, S.; Bharat, D.; Sengupta, S. A mathematical model of tumor-endothelial interactions in a 3D co-culture. Sci. Rep. 2019, 9, 1–14. [Google Scholar] [CrossRef]
  12. Rey-Barroso, J.; Calovi, D.S.; Combe, M.; German, Y.; Moreau, M.; Canivet, A.; Wang, X.; Sire, C.; Theraulaz, G.; Dupré, L. Switching between individual and collective motility in B lymphocytes is controlled by cell-matrix adhesion and inter-cellular interactions. Sci. Rep. 2018, 8, 1–16. [Google Scholar]
  13. Plou, J.; Juste-Lanas, Y.; Olivares, V.; Del Amo, C.; Borau, C.; García-Aznar, J. From individual to collective 3D cancer dissemination: Roles of collagen concentration and TGF-β. Sci. Rep. 2018, 8, 1–14. [Google Scholar] [CrossRef]
  14. Koh, I.; Cha, J.; Park, J.; Choi, J.; Kang, S.G.; Kim, P. The mode and dynamics of glioblastoma cell invasion into a decellularized tissue-derived extracellular matrix-based three-dimensional tumor model. Sci. Rep. 2018, 8, 1–12. [Google Scholar] [CrossRef] [Green Version]
  15. Ling, L.; Mulligan, J.A.; Ouyang, Y.; Shimpi, A.A.; Williams, R.M.; Beeghly, G.F.; Hopkins, B.D.; Spector, J.A.; Adie, S.G.; Fischbach, C. Obesity-Associated Adipose Stromal Cells Promote Breast Cancer Invasion through Direct Cell Contact and ECM Remodeling. Adv. Funct. Mater. 2020, 30, 1910650. [Google Scholar] [CrossRef]
  16. Alizadeh, A.A.; Aranda, V.; Bardelli, A.; Blanpain, C.; Bock, C.; Borowski, C.; Caldas, C.; Califano, A.; Doherty, M.; Elsner, M.; et al. Toward understanding and exploiting tumor heterogeneity. Nat. Med. 2015, 21, 846–853. [Google Scholar] [CrossRef]
  17. Gao, Y.; Zhang, H.; Lirussi, F.; Garrido, C.; Ye, X.Y.; Xie, T. Dual inhibitors of histone deacetylases and other cancer-related targets: A pharmacological perspective. Biochem. Pharmacol. 2020, 182, 114224. [Google Scholar] [CrossRef] [PubMed]
  18. Yin, A.; Moes, D.J.A.; van Hasselt, J.G.; Swen, J.J.; Guchelaar, H.J. A review of mathematical models for tumor dynamics and treatment resistance evolution of solid tumors. CPT Pharmacometrics Syst. Pharmacol. 2019, 8, 720–737. [Google Scholar] [CrossRef] [Green Version]
  19. Chaudhary, R.K.; Chaurasiya, V.; Singh, J. Numerical estimation of temperature response with step heating of a multi-layer skin under the generalized boundary condition. J. Therm. Biol. 2022, 108, 103278. [Google Scholar] [CrossRef] [PubMed]
  20. Chaudhary, R.K.; Chaurasiya, V.; Awad, M.M.; Singh, J. A numerical study on the thermal response in multi-layer of skin tissue subjected to heating and cooling procedures. Eur. Phys. J. Plus 2022, 137, 1–18. [Google Scholar] [CrossRef]
  21. Cristini, V.; Lowengrub, J. Multiscale Modeling of Cancer: An Integrated Experimental and Mathematical Modeling Approach; Cambridge University Press: Cambridge, UK, 2010. [Google Scholar]
  22. Byrne, H.M.; Chaplain, M. Growth of necrotic tumors in the presence and absence of inhibitors. Math. Biosci. 1996, 135, 187–216. [Google Scholar] [CrossRef]
  23. Amereh, M.; Edwards, R.; Akbari, M.; Nadler, B. In-silico modeling of tumor spheroid formation and growth. Micromachines 2021, 12, 749. [Google Scholar] [CrossRef] [PubMed]
  24. Amereh, M.; Bahri, Y.; Edwards, R.; Akbari, M.; Nadler, B. Asymmetric Growth of Tumor Spheroids in a Symmetric Environment. Mathematics 2022, 10, 1955. [Google Scholar] [CrossRef]
  25. Li, H.; Peng, R.; Wang, Z.A. On a diffusive susceptible-infected-susceptible epidemic model with mass action mechanism and birth-death effect: Analysis, simulations, and comparison with other mechanisms. SIAM J. Appl. Math. 2018, 78, 2129–2153. [Google Scholar] [CrossRef]
  26. Grimes, D.R.; Kannan, P.; McIntyre, A.; Kavanagh, A.; Siddiky, A.; Wigfield, S.; Harris, A.; Partridge, M. The role of oxygen in avascular tumor growth. PloS ONE 2016, 11, e0153692. [Google Scholar] [CrossRef] [Green Version]
  27. Greenspan, H.P. Models for the Growth of a Solid Tumor by Diffusion. Stud. Appl. Math. 1972, 51, 317–340. [Google Scholar] [CrossRef]
  28. Byrne, H. Modelling the role of cell-cell adhesion in the growth and development of carcinomas. Math. Comput. Model. 1996, 24, 1–17. [Google Scholar] [CrossRef]
  29. Casciari, J.J.; Sotirchos, S.V.; Sutherland, R.M. Variations in tumor cell growth rates and metabolism with oxygen concentration, glucose concentration, and extracellular pH. J. Cell. Physiol. 1992, 151, 386–394. [Google Scholar] [CrossRef]
  30. De Araujo, A.L.; Fassoni, A.C.; Salvino, L.F. An analysis of a mathematical model describing acid-mediated tumor invasion. Math. Methods Appl. Sci. 2019, 42, 6686–6705. [Google Scholar] [CrossRef] [Green Version]
  31. Daher, F.B.; Chen, Y.; Bozorg, B.; Clough, J.; Jönsson, H.; Braybrook, S.A. Anisotropic growth is achieved through the additive mechanical effect of material anisotropy and elastic asymmetry. Elife 2018, 7, e38161. [Google Scholar] [CrossRef]
  32. Ramírez-Torres, A.; Rodríguez-Ramos, R.; Merodio, J.; Bravo-Castillero, J.; Guinovart-Díaz, R.; Alfonso, J.C.L. Mathematical modeling of anisotropic avascular tumor growth. Mech. Res. Commun. 2015, 69, 8–14. [Google Scholar] [CrossRef]
  33. Ramírez-Torres, A.; Rodríguez-Ramos, R.; Merodio, J.; Penta, R.; Bravo-Castillero, J.; Guinovart-Díaz, R.; Sabina, F.J.; García-Reimbert, C.; Sevostianov, I.; Conci, A. The influence of anisotropic growth and geometry on the stress of solid tumors. Int. J. Eng. Sci. 2017, 119, 40–49. [Google Scholar] [CrossRef]
  34. Amereh, M.; Akbari, M.; Nadler, B. In-Silico Study of Asymmetric Remodeling of Tumors in Response to External Biochemical Stimuli. Sci. Rep. 2022, 13, 1–14. [Google Scholar] [CrossRef]
  35. Jiang, Y.; Pjesivac-Grbovic, J.; Cantrell, C.; Freyer, J.P. A multiscale model for avascular tumor growth. Biophys. J. 2005, 89, 3884–3894. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  36. Hatzikirou, H.; Breier, G.; Deutsch, A. Cellular automaton modeling of tumor invasion. In Complex Social and Behavioral Systems: Game Theory and Agent-Based Models; Springer: Berlin/Heidelberg, Germany, 2020; pp. 851–863. [Google Scholar]
  37. Kim, Y.; Stolarska, M.A.; Othmer, H.G. A hybrid model for tumor spheroid growth in vitro I: Theoretical development and early results. Math. Model. Methods Appl. Sci. 2007, 17, 1773–1798. [Google Scholar] [CrossRef] [Green Version]
  38. Hughes-Alford, S.K.; Lauffenburger, D.A. Quantitative analysis of gradient sensing: Towards building predictive models of chemotaxis in cancer. Curr. Opin. Cell Biol. 2012, 24, 284–291. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  39. Liu, P.; Shi, J.; Wang, Z.A. Pattern formation of the attraction-repulsion Keller-Segel system. Discret. Contin. Dyn. Syst. B 2013, 18, 2597. [Google Scholar] [CrossRef]
  40. Shojaei, S.; Basso, J.; Amereh, M.; Alizadeh, J.; Dehesh, T.; Rosa, S.D.S.; Clark, C.; Hassan, M.; Tomczyk, M.; Cole, L.; et al. A multi-omics analysis of glioma chemoresistance using a hybrid microphysiological model of glioblastoma. bioRxiv 2022. [Google Scholar] [CrossRef]
  41. Anderson, D.M.; McFadden, G.B.; Wheeler, A.A. Diffuse-interface methods in fluid mechanics. Annu. Rev. Fluid Mech. 1998, 30, 139–165. [Google Scholar] [CrossRef] [Green Version]
  42. Cahn, J.W.; Hilliard, J.E. Free energy of a nonuniform system. I. Interfacial free energy. J. Chem. Phys. 1958, 28, 258–267. [Google Scholar] [CrossRef]
  43. Talkenberger, K.; Cavalcanti-Adam, E.A.; Voss-Böhme, A.; Deutsch, A. Amoeboid-mesenchymal migration plasticity promotes invasion only in complex heterogeneous microenvironments. Sci. Rep. 2017, 7, 1–12. [Google Scholar] [CrossRef] [Green Version]
  44. Hecht, I.; Bar-El, Y.; Balmer, F.; Natan, S.; Tsarfaty, I.; Schweitzer, F.; Ben-Jacob, E. Tumor invasion optimization by mesenchymal-amoeboid heterogeneity. Sci. Rep. 2015, 5, 1–12. [Google Scholar]
  45. Vaziri-Gohar, A.; Cassel, J.; Mohammed, F.S.; Zarei, M.; Hue, J.J.; Hajihassani, O.; Graor, H.J.; Srikanth, Y.V.; Karim, S.A.; Abbas, A.; et al. Limited nutrient availability in the tumor microenvironment renders pancreatic tumors sensitive to allosteric IDH1 inhibitors. Nat. Cancer 2022, 4, 1–14. [Google Scholar] [CrossRef] [PubMed]
  46. Nagar, S.; Argikar, U.A.; Tweedie, D.J. Enzyme kinetics in drug metabolism: Fundamentals and applications. In Enzyme Kinetics in Drug Metabolism; Springer: Berlin/Heidelberg, Germany, 2014; pp. 1–6. [Google Scholar]
Figure 1. Three-dimensional tumor model to study tumor invasion and dissemination. (a) Invasion of patient-derived glioblastoma mutiforme (pdGBM) into patient-derived extracellular matrix (pdECM) and collagen; reproduced with permission [14]. PdECM facilitates the dissemination of cells into the matrix while collagen enhances the collective migration of cells. (b) Invasion of MCF10AT1 spheroids cocultured with patient-derived adipose stromal cells (ASCs) in collagen; reprinted with permission from [15]. 2020, John Wiley and Sons. Both cases show that the collective migration of cells led to the formation of a cluster of cells that preserved a distinct interface with the tumor core and the surrounding matrix, shown by red arrows.
Figure 1. Three-dimensional tumor model to study tumor invasion and dissemination. (a) Invasion of patient-derived glioblastoma mutiforme (pdGBM) into patient-derived extracellular matrix (pdECM) and collagen; reproduced with permission [14]. PdECM facilitates the dissemination of cells into the matrix while collagen enhances the collective migration of cells. (b) Invasion of MCF10AT1 spheroids cocultured with patient-derived adipose stromal cells (ASCs) in collagen; reprinted with permission from [15]. 2020, John Wiley and Sons. Both cases show that the collective migration of cells led to the formation of a cluster of cells that preserved a distinct interface with the tumor core and the surrounding matrix, shown by red arrows.
Symmetry 15 00283 g001
Figure 2. The equilibrium state of separated phases was acquired as f ^ / C ¯ and had two extrema, (a) left axis, and the corresponding flux function J ^ = D f grad ( f ^ / C ¯ ) had two nontrivial roots, (a) right axis. This allowed for the coexistence of high C h and low C l concentrations with a continuous transition (b).
Figure 2. The equilibrium state of separated phases was acquired as f ^ / C ¯ and had two extrema, (a) left axis, and the corresponding flux function J ^ = D f grad ( f ^ / C ¯ ) had two nontrivial roots, (a) right axis. This allowed for the coexistence of high C h and low C l concentrations with a continuous transition (b).
Symmetry 15 00283 g002
Figure 3. Schematic representation of a typical in vitro tumor platform used to study the growth and invasion of tumors. Tumor-embedded hydrogels were injected into the chamber of a microfluidic chip supplied with nutrients via microchannels that resembled the vasculatures. The size of the chip is typically in the order of a few centimeters, and the smallest dimension of the channels is around a few millimeters.
Figure 3. Schematic representation of a typical in vitro tumor platform used to study the growth and invasion of tumors. Tumor-embedded hydrogels were injected into the chamber of a microfluidic chip supplied with nutrients via microchannels that resembled the vasculatures. The size of the chip is typically in the order of a few centimeters, and the smallest dimension of the channels is around a few millimeters.
Symmetry 15 00283 g003
Figure 4. Steady state distribution of cells in absence of nutrients for β ¯ = { 0.1 , 1 , 2 , 4 } .
Figure 4. Steady state distribution of cells in absence of nutrients for β ¯ = { 0.1 , 1 , 2 , 4 } .
Symmetry 15 00283 g004
Figure 5. Cell distribution corresponding to different values of γ , V ¯ , η ¯ , ψ and c. (a,b) Higher values of γ and V ¯ increased the growth of cells moving outward the tumor core, yielding a peak in cell concentration with a distinct interface. (c,d) Higher values of η ¯ and ψ reduced cell concentration and eliminated the peak. Parameter c controlled the formation of a peak and the concentration of cells within the invasion zone. (e) A higher c corresponds to a higher sensitivity of the proliferation rate to the gradient of nutrients giving rise to a larger peak at the leading interface.
Figure 5. Cell distribution corresponding to different values of γ , V ¯ , η ¯ , ψ and c. (a,b) Higher values of γ and V ¯ increased the growth of cells moving outward the tumor core, yielding a peak in cell concentration with a distinct interface. (c,d) Higher values of η ¯ and ψ reduced cell concentration and eliminated the peak. Parameter c controlled the formation of a peak and the concentration of cells within the invasion zone. (e) A higher c corresponds to a higher sensitivity of the proliferation rate to the gradient of nutrients giving rise to a larger peak at the leading interface.
Symmetry 15 00283 g005aSymmetry 15 00283 g005b
Figure 6. The dynamic progression of tumor invasion within a hydrogel matrix in the form of the shell-type motion of invasive cells over time. The shell-type cluster of cells moved up the gradient of nutrients due to chemotaxis while sustaining its structure due to adhesive cell–cell forces. The ability of cells to preserve the shell as they move was modeled with nonconvex free energy that included interfacial stress.
Figure 6. The dynamic progression of tumor invasion within a hydrogel matrix in the form of the shell-type motion of invasive cells over time. The shell-type cluster of cells moved up the gradient of nutrients due to chemotaxis while sustaining its structure due to adhesive cell–cell forces. The ability of cells to preserve the shell as they move was modeled with nonconvex free energy that included interfacial stress.
Symmetry 15 00283 g006
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

Amereh, M.; Struchtrup, H.; Nadler, B. Mathematical Modeling of Spherical Shell-Type Pattern of Tumor Invasion. Symmetry 2023, 15, 283. https://doi.org/10.3390/sym15020283

AMA Style

Amereh M, Struchtrup H, Nadler B. Mathematical Modeling of Spherical Shell-Type Pattern of Tumor Invasion. Symmetry. 2023; 15(2):283. https://doi.org/10.3390/sym15020283

Chicago/Turabian Style

Amereh, Meitham, Henning Struchtrup, and Ben Nadler. 2023. "Mathematical Modeling of Spherical Shell-Type Pattern of Tumor Invasion" Symmetry 15, no. 2: 283. https://doi.org/10.3390/sym15020283

APA Style

Amereh, M., Struchtrup, H., & Nadler, B. (2023). Mathematical Modeling of Spherical Shell-Type Pattern of Tumor Invasion. Symmetry, 15(2), 283. https://doi.org/10.3390/sym15020283

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