Vladimir Kalinin

Cell – extracellular matrix interaction in glioma growth. In silico model

De Gruyter | 2020


The study aims to investigate the role of viscoelastic interactions between cells and extracellular matrix (ECM) in avascular tumor growth. Computer simulations of glioma multicellular tumor spheroid (MTS) growth are being carried out for various conditions. The calculations are based on a continuous model, which simulates oxygen transport into MTS; transitions between three cell phenotypes, cell transport, conditioned by hydrostatic forces in cell–ECM composite system, cell motility and cell adhesion. Visco-elastic cell aggregation and elastic ECM scaffold represent two compressible constituents of the composite. Cell–ECM interactions form a Transition Layer on the spheroid surface, where mechanical characteristics of tumor undergo rapid transition. This layer facilitates tumor progression to a great extent. The study demonstrates strong effects of ECM stiffness, mechanical deformations of the matrix and cell–cell adhesion on tumor progression. The simulations show in particular that at certain, rather high degrees of matrix stiffness a formation of distant multicellular clusters takes place, while at further increase of ECM stiffness subtumors do not form. The model also illustrates to what extent mere mechanical properties of cell–ECM system may contribute into variations of glioma invasion scenarios.

1 Introduction

Glioblastoma multiform (GBM) is one of the most lethal types of human cancers with a median survival time slightly above a year [1]. Rapid proliferation and aggressive invasion into surrounding normal tissue is a hallmark of GBM. This ability of long range infiltration of the primary tumor cells through brain tissue even in the absence of distant metastases, makes GBM a primary target of many research projects including in vitro studies and computer models [1], [2], [3].

It has been found [4], [, 5] that purely mechanical characteristics of glioma microenvironment may play pivotal role in the invasion process. The importance of cell–extracellular matrix (ECM) interaction and the role of matrix density, elasticity and porosity have emerged as a particularly strong regulators of migration and proliferation rate of malignant brain tumor cells. There are number of observable trends reported by different authors [3], [7], [8], [9].

The invasive distance was found to be typically biphasic function of ECM pore size increasing in the range of small pore sizes and decreasing at wider pores [7], [, 10]. The maximum invasiveness is shifting towards wider pore size with increase of ECM stiffness [7]. Glioma cells exhibit increase of cell motility in narrow channels and may demonstrate ability of amoeboid squeezing through narrow pores with sizes below their nuclear diameter [11], [, 12]. There are certain difficulties in generalization of those observations in order to reduce them to simple dependences on ECM porosity or density. Multiple attributes of the matrix microstructure, such as channel configuration (width/length ratio), fiber thickness, fiber alignment and crosslink density may significantly affect cell invasiveness and should be taken into account in data analysis [5], [7], [13], [14].

An increase of matrix density, which affects its stiffness as monotonic function, leads typically to formation of denser spheroids with increased invasiveness [15], [16], [17]. ECM stiffness has been found to influence strongly glioma invasiveness [4], [5], [6], [7], [18], [19], [20]. On the other hand, highly elastic glioma cell lines demonstrate more aggressive behavior when compared with stiffer cells cultured in similar substrate [21]. Based on the reported results, one can conclude that the ratio of ECM- to -cells stiffness has a strong modulatory effect on glioma cell invasion and migration characteristics. Typically, a system of ‘soft cells in stiff matrix’ demonstrates high invasive potential.

Nevertheless, about one quarter of the studied cell lines exhibit mechanosensitivity to ECM stiffness different from that typically observed, demonstrating biphasic function of rigidity-dependence [19] or rigidity-insensitive behavior [5]. Interestingly, these cells insensitive to the matrix stiffness exhibited the highest invasive capacity even in soft matrigel 3D cultures [5]. Therefore, in spite of the general trends every case is unique regarding the behavior of each specific glioma cell line. It was also shown in [4], [, 19] that glioma cell behavior, being migratory dependent on purely mechanical effects of stresses and ECM structure, may at the same time be mediated though more complicated pathways involved in glioma invasion process.

The role of ECM stiffness in glioma growth can be partially explained, as it conditions mechanical stresses in tumor microstructure, which has dramatic influence on the growth dynamics. The mutual influence of tumor cells and the microenvironment through the mechanical stresses can modulate migration and therefore invasiveness of tumor [22], [23], [24] and also regulate proliferation and apoptosis of tumor cells [25], [, 26]. Despite the observed correlations of tumor growth with the stresses and cell–ECM rheology, the precise role of the microenvironment in this process is still not well studied [3], [4], [9], [27].

Not only values, but a special form of the local stresses distribution across tumor border provides its progression. That particular distribution has been observed in glioma MTS growth where pushing and traction forces have been measured along MTS radius across its border [28]. The authors suggested that the traction forces may cause direct migration through tracks formed by ECM alignment that further promotes elongation of tumor cells along the traction line. In so doing, the growing MTS exerts an outward ECM pressure at the distances greater than its own size. Similar “pushing-and-pulling” scenario across the edge of MTS core has been demonstrated on murine colon carcinoma cell line by Kopanska and colleagues [29]. The authors showed ECM radial contraction towards MTS surface and tangential alignment of the gel fibers close to the surface. Generally, the presence of traction forces exerted by cells in a 3D matrix was observed in number of works e.g. [15], [29], [30] and measured in the range of 100 and 250 pN.

One of the key questions associated with these observations would be to what extent the observed scenarios are specific to particular cancer cells or particular type of ECM structure? Regarding glioma cells, it was shown that different cell lines may form substantially different invasive patterns [6], [31], [32], resulting from different invasive strategies. Therefore, generalization of the results obtained for a certain cell line to all glioma cells may be seriously restricted by significant variation of mechanobiological properties between different cell lines. One such difference is elasticity of different glioma cell lines, as shown by Louca and colleagues [21].

The relative contributions of stress or morphology related attributes into tumor progression are difficult to rectify in experiments, as they are often interconnected, do not allow independent variation and may affect other physical characteristics of cell–ECM system, such as focal adhesion and matrix stiffness [6], [, 19]. In silico models, being free from these limitations may seriously assist in further analysis of glioma invasiveness. Computer simulation can be particularly efficient in study of purely mechanical properties of tumor, based on clear and well known principles.

Rheological tumor properties have been modeled initially as a single homogeneous continuum [33] and later within multi-phase approach [34], [35], [36], [37]. A general separation of the tension zone in the center of a tumor spheroid from the zone of compression in the spheroid periphery has been shown within the approach of elastic incompressible continuum [33]. The only source of stress in the model was cell proliferation resulting in tumor expansion. The multi-phase models are generally focused on interaction between the phases. Cells inevitably represent one of the phases, while ECM most of the time represents another one.

Two elastic incompressible phases interact in the model of Ngwa and Agyingi [35]. Spherically symmetric cell agglomeration – phase I is expanding and exerting pressure on an external medium – phase II. In return the phase II offers certain resistance to the growing tumor, which is relatively small comparing with the stress effect related to its growth. Nevertheless the resistance of the surrounding medium clearly influenced the saturation size of the tumor. It was also concluded that the conditions of nutrient diffusion may impose significant limitations to the nutrient supply and therefore influence the tumor saturation size.

Further to demonstration of experimental evidence of cell compressibility [38], [, 39] this property has been included into in silico models. One such approach was used in Ref. [36], where cells modeled as elastic and compressible aggregation with Poisson factor ν = 0.2. ECM was modeled in the work as a rigid structure and this model was intended to study the growth of encapsulated MTS under external stress. The calculated radius of tumor growing in confined environment as function of time was in good correlations with experimental results of Ref. [38], as well as with the outcome of their own in silico model based on assumption of elastic extracellular scaffold. Therefore, the assumptions of rigid [36] and elastic [38] ECM did not show significant difference in simulated results. However, another numerical study [37] demonstrated actual difference in tumor development for these two alternatives. The cells have been modeled as uncompressible fluid, while ECM considered as a porous elasto-visco-plastic or rigid solid material [37]. The simulations based on these alternatives showed faster growth of tumor in elastic matrix in comparison with the case of rigid ECM, as the deformable ECM changes its porosity under the pressure exerted by growing tumor. The exercise demonstrated the importance of ECM deformability as a critical factor of tumor growth. The simulation results are in line with in vitro study of human colon, murine mammary cancer and rhabdomyosarcoma cells [39], where the MTS growth was insensitive to ECM stiffness up to a certain threshold when exhibited a sharp step down in the final spheroid sizes. Regarding glioma growth, however laboratory studies generally show an inverse system of ‘soft cells in rigid matrix’ as a condition for most aggressive tumor growth [4], [5], [6], [7], [17], [18], [19], [20], [21]. This comparison stresses the point that mechanobiological theory of tumor progression is still in its infancy and cell–ECM interaction scenarios need further study.

It was found that glioma cells follow “Go or Grow” (GoG) principle as common rule for growing tumor. The GoG concept proposes that cells are unlikely exhibiting both proliferation and migration concurrently and is based on observations of astrocytoma cells rarely dividing while moving [40]. While it was clearly demonstrated in further observations [41] that proliferation and migration of glioma cells are inversely correlated processes, the factors triggering switch of cell phenotype are still not quite clear. A shortage of nutrients is generally considered as factor of GoG switch, causing metabolic stress in glioma cells [42]. In particular, there are many results showing oxygen concentration as critical parameter for the transition of cell phenotype. The “hypoxic” concept of switching is widely exploited in models and analyses of glioma and general tumor growth [43], [44], [45], [46], [47]. Another nutrient often considered as critical for GoG switch is glucose [48]. Various microenvironmental conditions or a threshold of high cell density were also assumed to trigger the change of cell phenotype [40], [49], [50], [51]. And finally some factors particularly specific to glioma, such as decrease of astrocytes density on tumor periphery below certain level [52] may trigger the transition.

The nutrient based transition criteria are mutually connected, as the nutrient profiles exhibit similar shapes. The profiles are generally determined by strong concentration gradients along the MTS radius resulting from diffusion into the area of high cell density and therefore increasing nutrient consumption [53]. The environmental transition triggers do not have obvious connection to nutrient based thresholds and therefore represent different model of tumor growth. Completely apart from the considered above switch conditions is assumption that the phenotype transition is purely stochastic event [54]. This approach is justified to certain extent by poor current knowledge of the actual conditions leading to cell phenotype transition and the absence of a general concept for phenotype triggering criteria. The switch criteria of different nature may generate substantially different invasive scenarios [44], [, 53].

Another important factor of GoG dichotomy is probably the overlapping of phenotype characteristics. The main question here is whether migration and proliferation states exclude each other completely, strictly “Go-or-Grow”, as considered in Refs. [45], [, 54] or, if not then to what extent these inversely correlated features can overlap to be exhibited at the same time by the same cell, as in Ref. [46]. Two polar cases have been modeled in Ref. [44] where the case of “Go-and-Grow” features overlapped completely was compared with pure “Go-or-Grow” case. The results show “Go-or-Grow” strategy as more efficient in terms of glioma progression. However, a comparative analysis of the simple polar scenarios and the cases of partial ‘Go and Grow’ overlapping is still not carried out.

Generally in glioma growth the GoG dichotomy has significant influence on width of the invasive zone and tumor expansion dynamics, which is expected to be sensitive to transition characteristics including transition times and thresholds [44], [50], [53], [54]. The research results demonstrate how important role cell phenotype switch plays in tumor progression and how little at the same time we know about the criteria and characteristics of this process.

Another important factor influencing the growth of glioma is cell adhesion–repulsion mechanism [34], [50], [55], [56], [57], [58]. The interplay between cell phenotype transition dynamics and cell–cell interaction determines to great extent the shape and speed of invasive glioma front and its growth characteristics [44], [51], [55].

Pham and colleagues [51] considered cell–cell interactions by introducing density-dependent regulation of migration/proliferation dichotomy in a reaction–diffusion model and found significant influence of this interaction on the dynamics of tumor front. Number of simulations based on continuous and discrete models showed that glioma invasiveness is generally reduced when cell–cell attractive interactions grow up [58], [59], [60], [61]. These results were consistent with experimental data.

Cell–cell adhesion–repulsion mechanism also plays important role in formation of subtumors [55], [56], [57]. Cell–cell adhesion is responsible for formation of multicellular groups, strands and clusters forming at the tumor interface and growing away from the main tumor mass. It was shown in silico by Frieboes et al. [55] that subtumors are being formed by all cells of the tumor surface and not as a result of hyper-proliferation of a distant single cell [55]. This result demonstrates that cell motility is probably not the key factor in the formation of subtumors, when a single highly motile cell could invade far into nutrient-rich area and hyper-proliferate there. The actually observed scenario is that specific stress distribution forms hydrostatic pushing forces on the front providing all cells drift in a bias towards the invasion front regardless of their phenotype. This shows again that the distribution of stresses and therefore viscoelastic properties of tumor components including cell aggregation and ECM are strong modulators of glioma progression.

Summarizing the major outcome of the previous works one can conclude that variation of purely mechanical properties of ECM, such as its microstructure, deformability and compressibility, may significantly affect tumor growth. The experimental observations, briefly presented here and mostly dedicated to glioma, show dramatic influence of stress factor on its growth. The mechanism of this influence is presumably that stress can potentially modify the structure and density distribution of the major tumor components – ECM and cell aggregation, which in turn modulates tumor growth characteristics.

Overall, rheological properties of compressible tumor components, cell adhesion and cell phenotype transition are probably the key combination of mechanobiological attributes determining the regime and characteristics of glioma progression. All these attributes are included into the in silico model developed in the current work.

One of the objectives was to build a model demonstrating scenarios of cell–ECM interaction driven by mere mechanical principles of continuum that can be used as a basis in analyses of tumor growth processes observed in laboratory tests or a clear first approximation for further simulations of glioma variability, involving signaling pathways.

Another objective of the work was clarification of basic invasive scenarios in the frame of the model and their sensitivity primarily to matrix rheology and oxygen supply.

2 In silico model

We consider growing MTS as a composite system comprising two basic components: tumor сells and ECM. The interstitial fluids (IFs) are also part of this composite system. As first approximation, the interstitial fluids are assumed to represent an incompressible continuum, traveling through the matrix and cellular content without any resistance being always part of these two basic components.

Cell aggregation, as well as ECM, normally exhibits multiple voids [62]. Therefore, the “dry content” of the basic fractions is porous and compressible. IFs fill the voids of the cell aggregation and porous ECM scaffold, but theoretically can be squeezed out of these system components under external pressure, as they both are fully transparent for IFs. The details of IFs dynamics are not considered in the model.

In mathematical approximations, we formally consider only two system components: fluidized cell aggregation and fluidized ECM, assuming the sum of their volume fractions equal to 1.

ECM is considered to be purely elastic compressible scaffold, while the cell aggregation exhibits viscoelastic properties and also effects related to cell–cell adhesion, see the details below.

The model includes three cell phenotypes:

  • – Proliferative cells that can proliferate and have low or zero motility;

  • – Hypoxic cells that have low or zero proliferation rate and exhibit high motility;

  • – Necrotic cells, which are inert-do not exhibit motility and do not proliferate.

The transition between cell phenotypes is determined by two factors: concentration of oxygen as a switch threshold and characteristic transition times. Proliferative- to -hypoxic and inverse transitions are considered in the model along with hypoxic- to -necrotic switch. Removal of necrotic cells is not taken into account as not important process within current simulation times.

The model is based on continuum approximation. Some of the dynamic parameters of the model are being further obtained under the following assumptions:

  1. (a)

    Cancer cells comprise incompressible viscoelastic spheres. In general case, the cell colony may include multiple voids, having relatively low average cells concentration. A massive cell colony is isotropic visco-elastic compressible agglomeration of spheres.

  2. (b)

    Extracellular Matrix comprises isotropic, compressible, linear elastic, porous material. Cells proliferation and their motion through this scaffold produce elastic stresses in ECM, which may change the transport conditions for these cells.

  3. (c)

    Spherically symmetric simulation domain contains two basic components: cells and ECM, and does not present any phenomenologically determined spatial subdivision of its physical content into different parts (e.g. free moving boundaries with preset shedding rates out of them). Only initial conditions determine a clear boundary of the initial tumor spheroid. No further definitions of the border between the core spheroid and ECM are imposed and the distribution of partial densities in the composite system is formed in self-consistent simulation.

  4. (d)

    The cells exert attractive and repulsive forces on each other depending on intercellular distances. The lower limit of cell concentration in the colony is determined by the critical distance between them. Beyond this threshold cells cannot sense or forcibly affect each other.

  5. (e)

    The total stress distribution in the composite system follows iso-strain approximation. This intercellular force field together with ECM stress gradients induces cells motion tending to minimize the total potential energy of the system.

  6. (f)

    Cell motility is modeled through their diffusion. The diffusion rate is specific for each cell phenotype.

2.1 Basic processes in details

2.1.1 Attraction and repulsion of cells, cell compressibility

There are rather limited data available in literature on compressibility measurements for cell aggregation, as a separate phase independent of ECM. A combination of modeling and experimental measurement of cellular aspiration in a micropipette gives typical Poisson’s ratios of νC = 0.3–0.4 [63] in contrast to typically assumed value of 0.5 [64]. This may result from multiple discontinuities in cell clusters that could provide certain potential for compression. Therefore, apart from single cell compressibility, a cell colony should represent rather compressible aggregation. This presumption is confirmed by in vitro studies of MTS growth [40], which showed that general compressive stress may lead to an increase of cellular packing density. However, even a study of single cell size variations under pressure showed substantial compressibility. Melanoma cells for instance exhibit KC ≈ 2.5 kPa [38]. The authors demonstrated variation of the available data for bulk moduli of different cancer cells in the range 2.5–30 kPa. At the same time, the other sources show extremely low general cell compressibility KC = 2.5 × 109 Pa [65]. Having no specific data on a single glioma cell compressibility we are not actually bound to this particular characteristic. The current model does not consider dense cell aggregations behaving as a homogeneous substance with bulk modulus close to same of a single cell. The porosity of glioma cell aggregation always contributes into its compressibility. The correspondent rheological parameters are being used in the model, see section 2.3 “Choice of Parameters”. This compressible cell colony can be described in first approximation as accumulation of visco-elastic incompressible spheres. Each sphere consists of homogeneous elastic medium with a smooth boundary (similar to that assumed in Hertz or Sneddon models) covered with a layer of membrane protrusions providing cell–cell adhesion similar to “brush” models [66].

Table 1:

Mobility and random motility of Glioma Cells.

Characteristic Value Source
Collagen permeability, k 10−12 m2 [85]
Dynamic viscosity of cancer cells, μ 300 Pa s [83]
Mobility of cells in collagen, M 2.9 × 10−4 mm2/(Pa d) Current. Calc.
Hypoxic cells, Dm 2.0 × 10−2 mm2/day [47]
Proliferative cells, Dp Dp = Dh/10
Table 2:

Rheological parameters.

Parameter Value Source
Poisson ratio for collagen, νECM 0.13 [89]
Young’s modulus for collagen, EECM 0.5–2.5 kPa [89]
Bulk modulus for collagen, KECM 0.3–1.13 kPa Current. Calc.a
Poisson ratio for cells, νC 0.4 [63]
Young’s modulus for cells, EC 300 Pa [91]
Bulk modulus for cells, KC0 500 Pa Current. Calc.a

Table 3:

Oxygen transport and consumption.

Characteristic Value Source
Oxygen diffusion coefficient DO2 86.4 mm2/day [94]
Oxygen diffusion coefficient DO2a 40.0 mm2/day Used for comparison purposes
Oxygen uptake for proliferative cells Cp
αp [(mm3/Cell)mg/(L*Day)] αp = 1.38 × 10−2 [47]
Oxygen uptake for hypoxic cells Cm: αh αh = αp/5 [47]
Critical oxygen concentrations
Threshold to hypoxia ξh = 1.0 mg/L [99]
Threshold to necrosis ξn = 0.8 mg/L [99]
Table 4:

Proliferation and phenotype transition times.

Characteristic times Value Source
Proliferative to hypoxic phenotype τph = 1 h [47]
Hypoxic to proliferative phenotype τhp = 96 h [47]
Hypoxic to necrotic phenotype τhn = 32 h [52]
Proliferative cells doubling time τp = 24 h [47]
Hypoxic cells doubling time τm = 48 h [47]
Hypoxic cells death time τn = 48 h [47]

A theoretical maximum Cmax of cell density is fixed and determined by the volume of a single cell. Taking the figure VCell = 1200 μm3 of the EMT6/Ro rat brain tumor cell volume [67] as typical, we can assess Cmax as 8.3 × 105 mm−3. This dense cell packing does not include any fluids or ECM. The volume fraction of cells in MTS can be defined as ΦC = Ctot/Cmax, where the concentration Ctot includes all types of cells.

Separation distances between sphere centers determine the local cell density. An important characteristic of that cellular system is intercellular equilibrium distance req. In the absence of any other attractant but the cells, they should move in a way that minimizes the potential energy between them. Cells may move towards one another at r > req due to adhesion. In case of compression r < req, the cells should move due to increasing repulsive forces towards the area offering more free space. Another important parameter for the system of interacting spheres is the cut off distance r, which is the minimum separation distance between sphere centers when cells lose completely their mutual connection. In other words, at r > r any two cells in the model do not affect each other. At uniform cell distribution this situation corresponds to cell density Ccell < C.

2.1.2 Parameters of cell–cell interaction

As long as the current numerical model operates with cellular continuum all critical parameters of the discrete cell agglomeration can be interpreted in the form of cell concentrations. An equilibrium cell concentration Ceq corresponds to equilibrium separation distances between sphere centers r = req, so that at cell densities Ccell > Ceq, the cell agglomeration undergoes expansion and otherwise at Ccell < Ceq – contraction.

Experimental results for glioma MTS [70] show the inner quasi-stationary cell density in dynamically growing spheroid about Ceq = 4 × 105 cells/mm3, which agrees well with their assessment of Ceq based on the volume VC = 1200 μm3 of EMT6/Ro tumor cells. It was assumed that cell spheres in the colony at their equilibrium density are in direct contact with maximum number of neighbors, but still not stressed by these neighbors. Slightly lower equilibrium volume fraction of 0.39 for tumor cells was taken in Ref. [69] for a model of GBM growth, which corresponds to Ceq = 3.2 × 105 cells/mm3. We used value Ceq = 3.5 × 105 cells/mm3 in the current work. The Force–Distance (FD) curves for human glioma cells show maximum attraction force at ∼3 μm separation [70]. This force exhibits ten-fold decrease at separation distances about 20 μm and at 50 μm the cells are regarded as detached [70]. Assuming 20 μm separation between cell surfaces is a technical cut off, one can obtain the cut off distance of r  33 μm between the centers of spherical cells once VCell = 1200 μm3. Then the average cut off concentration of spherical cells can be assessed as C = 2.8 × 104 cells/mm3, while maximum attraction force should correspond to ratt ≈ 16.2 μm and Catt ≈ 2.35 × 105 cells/mm3.

The model of spheres is used here only to assess rheological parameters of a cell colony. For the purpose of numerical simulation, the cell colony is described as visco-elastic continuum with density distribution C(r) and ability to build inner positive pressure, as well as inner contractile forces. Having determined the key parameters C, Ceq and Cmax, we can further consider cell colony as visco-elastic continuum with bulk modulus KC defined as:

K C ( C tot , C eq ) = C tot d σ C d C tot
where C tot is total cell concentration including all cell phenotypes.

This further gives linear approximation for intercellular stress σC

(1) σ C = K C 0 ( C tot C eq 1 ) ,
which is feasible within limited vicinity of equilibrium value C eq, when C tot  <C max. Here, the constant K C0 is the bulk modulus at equilibrium cell concentration C tot  = C eq. Obviously, this linear approach is not valid at high local raise of cell concentrations C tot  ∼ C max. It is reasonably assumed that tumor expansion process has quasi-static character, when C tot does not grow too high above local equilibrium density determined by cell interactions along with the stress of extracellular matrix. This assumption has been justified by further MTS growth calculations. The other works e.g. [ 68], [ ,  69] also show smooth gradients of cellular density in GBM tumor, which typically correspond to quasi-static growth character and smooth time isolines.

In the case of volume expansion Ctot ≤ Ceq the expression (1) yields growing negative stress. This corresponds to an increase of intercellular distances at decreasing cell concentration.

As determined above, maximum attraction force fmax between the two cells is achieved at their homogeneous concentration Catt ≈ 2.35 × 105 cells/mm 3. And FD curve for human glioma cells shows fmax ∼ 430 pN at this maximum [70]. Further, assuming the number of neighboring cells roughly equal to 12, as kissing number for close-packing of equal spheres [71] and the surface area of a single sphere Scell ≈ 550 μm2 for VCell = 1200 μm3, we can assess the maximum contractile stress as max σatt ∼ 10 Pa. To satisfy contractile force characteristics at low cell densities another term γ should be introduced in (1) for Ctot ≤ Ceq, so that

(1a) σ C = γ K C 0 ( C tot C eq 1 ) ,
γ = { 1 C tot C eq 27 4 σ att K C 0 C eq ( C eq C ) 3 ( C tot C ) 2 C eq > C tot C 0 C tot < C

The approximation of contractile stresses at Ctot < Ceq is based on FD curve for glioma cells and does not depend on bulk modulus K0 for the cell aggregation. It satisfies characteristic parameters and general behavior of the curve.

2.1.3 Invasive cells motility

The motility of invasive cells is described as diffusion in the model and follows GoG concept. Proliferative cells can transition into hypoxic cells within characteristic time τpm if the local concentration of oxygen drops below hypoxic level. In turn, hypoxic cells can transition to necrotic cells within characteristic time τmn since the local oxygen concentration decreases further below necrotic barrier. At the same time, hypoxic cells can transition back to proliferative phenotype once the local oxygen concentration exceeds the transition barrier. A two switch approximations have been modeled in the present work.

Stiff switch approach strictly follows GoG rule, where all hypoxic cells exhibit invasive migration and are assumed non-proliferative, while all proliferative cells are non-invasive.

Overlapping switch allows certain motility of proliferative cells and attributes some proliferation to hypoxic cells as well. It was further shown in simulations that the “Overlapping switch” provides better agreement with experimental results. The switch parameters used for the calculations are presented in “Choice of Parameters” section.

The third phenotype-necrotic cells are totally inert. Total cell density is a sum of the three partial cell concentrations:

C tot = C p + C m + C n

Cp – concentration of proliferating cells;

Cm – concentration of hypoxic cells;

Cn – concentration of necrotic cells

2.1.4 Stresses in cells–ECM composite

MTS growth builds up stress gradients in the cells–ECM composite. This process conditions a self-consistent dynamics in the cell aggregation that represents motion of individual cells, as well as ECM displacements in the way that minimizes these stress gradients.

According to (1a), the local stress σc in the cell colony can be negative due to cell–cell adhesion, which may cause local contraction of the cell aggregation. Contrariwise, a local increase of total cell concentration above Ceq produces an increase of the local stress followed by expansion of the cell colony. ECM displacements make their own contribution into the distribution of total stress in the system σtot through building up stresses in ECM σECM. Cell invasion into ECM, as well as general growth of cell colony due to their proliferation is taken into account as a volumetric force that produces displacements of ECM, see equation (5).

The hydrostatically conditioned cell motion is described in the model as a drift with velocities vD proportional to the total stress gradient (2), which is typical for continuum models of avascular tumor [72], [, 73].

(2) v D = M σ tot r

σtot – total stress in the system. M is cell mobility that reflects to certain extent the effects of cell adhesion. As long as cell affinity to each other and to the matrix contributes to their viscosity one should reasonably expect a decrease of viscosity for necrotic cells [74]. However, the current model does not quantify the contribution of this effect into cell mobility and as a first approximation the mobility is taken constant for all cell phenotypes. The approximation of the total stress function σtot(σc, σECM) in the two-component media is based here on iso-strain model as most suitable for the description of transition layer (TL):

(3) σ tot = Φ c σ tot + ( 1 Φ c ) σ ECM

The Transition Layer can be described as spherical layer around external edge of MTS core, where the value of cell fraction ΦC drops down significantly, forming the front of spheroid expansion in the radial direction out of MTS center. The total cell concentration on MTS front exhibits a sharp decrease from C > Ceq down to C < Ceq. In other words, TL provides the transition from MTS core to invasive zone (IZ) represented by various forms of invasive patterns with low cell concentration compared to spheroid core. While the inner edge of TL is in the core, its outer edge is in IZ. As shown further in this work, TL plays significant role in MTS growth. This above description of MTS zones is purely phenomenological and based on data distributions obtained further in the current work, while no actual separation of MTS and its growing volume into zones or layers is assumed a priori by the current model.

2.1.5 ECM permeability

One of the objectives of the present study is to verify assumption that the ECM stress built by the cells traveling through the scaffold, may change the transport conditions for these cells. This change may trigger a positive feedback loop in cell invasion when the invading cells may affect ECM structure in a way that stimulates their further invasion.

A substantially non-uniform stress distribution across the MTS edge has been reported in number of works [28], [, 29]. One of the key ECM characteristics influenced by the stress is the permeability of its scaffold. Permeability of collagens and other porous scaffolds reduces at compressive stresses [74] and is reasonably expected to increase when tensile forces are applied to the scaffold volume. The case of compression is well studied for diffusion of fluids through porous scaffold in tissue engineering and study biomechanics of biomaterials in Orthopedics [74], [, 75]. The data obtained by these researchers can be used as an indication of ECM permeability for cells, as function of stress. The mean pore size and consequently scaffold permeability k was shown to decrease with increasing compressive strain ε for Matrigel [76], as well as collagen-glycosaminoglycan (CG) scaffolds [74], so that k ∼ (1 − ε)2.

The permeability of stretched porous matrixes is not so well studied compared to the case of compression. However, a measurement on CG scaffolds demonstrates well expected result where cell migration increases with increasing pore size [77]. Also, a study of stretched porous membranes [78] showed a positive shift in the mean pore size and permeability that accompanies membrane stretching.

We can extrapolate the simple function presented in Ref. [74] for the compressive strain to the case of applied tensile strain assuming the mean pore size of the stressed matrix proportional to l (1 − ε). Here the length l is a parameter of the porous microstructure and the strain ε is taken with negative sign for tensile case. Based on this simple approach, we can further describe the permeability of hydrogel porous matrix, undergoing tensile or compressive stress, as

k = k 0 ( 1 ϵ v 3 ) 2 ,
where ε v is ECM dilatation and k 0 is permeability of unstressed hydrogel.

ECM undergoes symmetric tangential stretch and radial compression in front of the core spheroid, which represents the transition layer. As shown further in simulations, the total ECM stress distribution within TL ranges from tensile strain level up to maximum compressive values (Figures 6 and 7). This specific stress distribution may affect the conditions of tumor progression through ECM permeability.

For the case of stress-dependent ECM permeability the function (2) can be rewritten in the form:

(4) v D = M 0 ( 1 ϵ v 3 ) 2 σ tot r ,

M0 is mobility of fluid in unstressed ECM M0 = k0/μ, where μ is the dynamic viscosity of the fluid.

2.2 Dynamics of cellular system

2.2.1 Governing equations

We consider static linear elasticity problem for spherically symmetric system with the origin at the center of the sphere. The cells, moving due to their motility and also being under hydrostatic pressure are forced through hydrogel scaffold and proliferate in the pores, which causes expansion of the composite system. Similar process of solid tumor dilatation, conditioned by cell production and the reaction of cell–ECM composite system on stresses is described in models [73], [, 79].

ECM is considered here as solid elastic porous component. Its expansion results in spherically symmetric tensile stress and non-uniform stress distribution along the MTS radius.

These stresses, together with the intercellular stresses caused by cell–cell interactions, form the total stress distribution in cells–ECM composite in accordance with iso-strain model (3). Gradients of this stress yield hydrostatic conditions for cell motion according to (2), which leads to stress equilibrium in the system.

Here we reasonably assume that the growth rate of Glioma cells is rather slow compared to the viscous relaxation time in cell–ECM composite-quasi-static expansion. Therefore, the system is in stress equilibrium and its dynamics is provided by cell and oxygen balance equations.

The equation for the stress components:

d σ r r d r + 1 r ( 2 σ r r σ θ θ σ ϕ ϕ ) = 0 ,
where σ rr , σ θθ, σ φφ are diagonal components of stress tensor.

Based further on strain-displacement and stress-strain relations, the governing equation for displacement u in the system [80], [81] can be written as

(5) r ( 1 r 2 r ( r 2 u ) ) = ( 1 + ν ) ( 1 ν ) δ v ,
where ν is ECM Poisson’s ratio; δ V  – the inner volumetric strain produced by expanding cell aggregation inside the matrix
δ v = 1 K ECM d σ C d r = γ K C 0 K ECM d d r ( C tot C eq )

KECM – ECM bulk modulus,

KC0 is the bulk modulus of the cell aggregation at equilibrium concentration Ctot = Ceq

C tot = C p + C m + C n ;

Cp – concentration of proliferating cells;

Cm – concentration of hypoxic cells;

Cn – concentration of necrotic cells.

Radial and tangential components of the strain tensor are as follows:

ϵ r r = d u / d r ; ϵ θ θ = ϵ φ φ = u / r

The components of ECM stress can be found as

(6) σ r r = K ECM ( 1 + ν ) ( ( 1 ν ) ϵ r r + 2 ν ϵ ϕ ϕ ) ,
(7) σ ϕ ϕ = K ECM ( 1 + ν ) ( ν ϵ r r + ϵ ϕ ϕ ) ,
as long as σ θθ  = σ φφ in spherically symmetric system the total stress in ECM
(8) σ ECM = 2 σ ϕ ϕ + σ r r 3 ,
where σ ECM – total stress in ECM.

While the system is in quasi equilibrium of stresses its temporal dynamics is governed by cell and oxygen balance equations. Cell balance obeys reaction–diffusion model supplemented by advection, as cells drift in the hydrostatic force field. Equations (9–12) describe proliferation, invasion and transition of phenotype for each of the three cell types.

(9) C i t + 1 r 2 r ( r 2 v D C i ) = 1 r 2 r [ D i r 2 r ( C i ) ] + S i
where C i is one of the three cell phenotypes C p, C m, C n, D i are diffusion coefficients of ith type of cell: D p, D m and obviously D n = 0; v D is drift velocity of cells, as defined by (2) for constant permeability of the scaffold and (4) for the case of stress-dependent ECM permeability. This velocity is the same for all types of cells in the model. S i represents sources S p, S m and S n for proliferative, hypoxic and necrotic cells correspondingly.
(10) S p = C p τ p ( 1 C tot C q ) n o x n max h p m C p τ p m + h m p C m τ m p ,
(11) S m = h m C m τ m ( 1 C tot C q ) n o x n max + h p m C p τ p m h m p C m τ m p h n e c C m τ n ,
(12) S n = h n e c C m τ n ,
where τ p and τ m are characteristic times of proliferation for proliferative and hypoxic cells, τ pm , τ mp , τ n are phenotype transition times: proliferative to hypoxic, hypoxic to proliferative and hypoxic to necrotic cells correspondingly. The maximum carrying capacity C q is considered here as the upper cut-off limit assessed above as C max = 8.3 × 10 5 mm −3. The source of proliferative cells S p represents logistic proliferation less transition of proliferative cells into hypoxic phenotype C p  -C m. Reverse phenotype transition C m  -C p, being rather slow process, also contributes into S p. The term n ox /n max represents the ratio of oxygen concentration to its maximum level in the domain, which currently corresponds to the initial O 2 concentration n 0. This ratio approximates the decrease of proliferation rate at reduced concentration of oxygen; the approach was also used by other authors (e.g. [ 69]). The last two transition processes form hypoxic cell source S m with the opposite sign along with the third term, which describes cell necrosis C m  -C n. The fraction of necrotic cells is formed through necrosis of hypoxic cells only. The proliferative cells do not undergo direct transition to necrotic state. A deprivation of oxygen is the only condition resulting in cell death in the model, which goes through hypoxic stage.

The switch factors hpm, hmp and hnec were introduced to provide a smooth approximation to the step function:

h p m ( n o x ) = { 1 , ( 1 tanh ( n o x ξ h ) λ ) / 2 0 , ( 1 tanh ( ξ h n o x ) λ ) / 2 , h m p ( n o x ) = 1 h p m ( n o x ) , h n e c ( n o x ) = { 1 , ( 1 tanh ( n o x ξ n ) λ ) / 2 0 , ( 1 tanh ( ξ n n o x ) λ ) / 2 ,
where ξ h and ξ n are thresholds of oxygen concentration for C p  -C m and C m  -C n phenotype transitions correspondingly; λ = 1 provided smooth transition.

Factor hm is defined according to one of the two approaches for cell phenotype transition exploited in the study. The first one assumes stiff “go or grow” approach, which means that proliferative cells exhibit no motility at all and hypoxic cells do not proliferate hm = 0.

The soft “go or grow” approach has also been exploited assuming that all non-necrotic cells exhibit certain motility and mitosis. However, the rate of proliferation drops down for hypoxic cells, as well as rather low motility characterizes proliferative cells.

The oxygen concentration nox ≤ ξn is considered as complete deprivation of oxygen (and therefore threshold to anoxia) and defines the switch factor hnec.

Oxygen is critical metabolite and the only nutrient in this model. Diffusion of oxygen and its local consumption by the cells determine the local concentration nox as per equation (13).

(13) n o x t = 1 r 2 r [ D o x r 2 r ( n o x ) ] ( α p C p + α m C m )
n ox is oxygen concentration, D ox – diffusion coefficient of oxygen, α p and α m – oxygen uptake for proliferative and hypoxic cells correspondingly.

Following (3) we can write for the total stress in the composite system:

(14) σ tot = ( C max C tot C max σ ECM + C tot C max σ C ) .

2.2.2 Initial and boundary conditions for the governing equations

Following Stein’s et al. experimental setup [68] the initial radius of the spheroid in simulations has been set to R0 = 250 μm and cell concentration Cinit = 7.0 × 104 cells/mm3.

We prescribe a uniform initial distribution of oxygen in the domain nox(r, t = 0) = n0 = 4 mg/L.

Spherical symmetry dictates symmetry conditions for all profiles Xi about r = 0: dXi/dr (r = 0) = 0.

We assume zero displacement and strain on the distant border of the domain r = R0:

U ( r = R 0 ) = 0 ϵ r r ( r = R 0 ) = 0

To prevent cells from leaving the domain no-flux boundary conditions were imposed for all cell phenotypes in equation (9): dCi/dr (r = R0) = 0.

Two types of BCs were considered for oxygen diffusion (13).

Type 1:

Closed volume or reflection boundary conditions (RBC)

This assumes the same condition on the domain border as set for cells

d n o x / d r ( r = R 0 ) = 0 .

Type 2:

Permanent supply (PSBC)

This condition assumes permanent supply of oxygen throw the border, so that

n o x = n 0 = const .

2.2.3 Numerical implementation

Finite differences approximation has been used to solve the system (5)(14). Equation (5) is discretized by a second order symmetric scheme. Crank–Nicolson scheme was implemented to approximate solution to the oxygen and cell diffusion problems (9) and (13). Time splitting approach is used to approximate phenotype transition and cell transport due to stress gradients. There was special requirement for conservative and stable algorithm to approximate the problem of hydrostatically forced cell drift in the condition when the force-field changes its sign along the MTS radius. The Phoenical LPE version of Flux Corrected Transport (FCT SHASTA) algorithm has been implemented to solve this problem [82]. This symmetric three point diffusive transport scheme provided second order of approximation, which is important to form a correct cell density profile of the invasive zone and also stable solution in stressed MTS core.

2.3 Choice of parameters

2.3.1 Mobility of glioma cells in collagen

Mobility is calculated as the value of collagen permeability k over dynamic viscosity of cells μ: M = k/μ. The results obtained for cultured anaplastic carcinoma cells [83] were taken as a reference value of cancer cells viscosity 300 Pa s.

There are dramatic variations in the experimental data for collagen permeability ranging from 2 × 10−16 m2 for 3% collagen [84] to 1.00 × 10−12 m2 [85]. Other authors typically show data within the above range, for example 10−15−10−14 m2 [86], [, 87]. Taking the top value of permeability 1.00 × 10−12 m2 we can obtain M = 3.3 × 10−15 m2/(Pa s) = 2.9 × 10−4 mm2/(Pa d). This choice of permeability value brought the results of simulation into correspondence with the experimental data.

2.3.2 Rheological parameters

Similar to permeability values the variations in the stiffness of collagens used for culturing GBM cells are also significant. Typically, GBM cells are being cultured in collagens of E = 1.0 and 35 kPa [6], [17], [88]. Some measurements show much lower values. For example, the examination of gel in the absence of cells used further for Glioma MTS experiments [15] showed Young’s moduli of E  100 Pa at concentrations 1.5–2.5 mg/mL. The study of acellular collagens in Ref. [89] gave typical range for Young’s modulus E = 0.5–12 kPa with Poisson ratio ν ∼ 0.13. In particular, acellular crosslinked gels of 3–5 mg/mL exhibited stiffness of 1–2 kPa. These detailed results were taken as a basis for our calculations. The range E = 0.5–2.5 kPa has been studied, which gives for bulk modulus K = 225–1130 Pa at ν = 0.13.

As discussed above, a cell colony represents compressible media. Compared to hydrogels, cancer cell aggregation exhibit lower compressibility with Poisson ratios about νC = 0.36–0.4 [63], [, 90].

Rheological parameters of porous cell aggregation without any other inclusions should ideally represent cells in the model as a component of the general composite system. Our evaluation of these parameters was based on stiffness of the Glial Cells Ecell ∼ 200–300 Pa [91] and Young’s moduli derived from the surface of fresh glioma samples 350 Pa [88]. We take νC = 0.4; E = 300 Pa and K = 500 Pa.

2.3.3 Cell diffusion in collagen

The diffusion coefficients for the model Dh are determined by the random motility of invasive hypoxic cells. The sources are being mostly referred to in literature show invasive diffusivity for hypoxic cells from Dh = 1.0 × 10−3 mm2/day [92] to Dh = 2.0 × 10−2 mm2/day [93].

In glioma MTS experiments of Stein et al. [68] these data were validated through comparison with the measured characteristics. It was found that the value of Dh = 2.0 × 10−2 mm2/day most closely described the measured MTS growth speed. This particular value was used in the current model along with typical approach for normal proliferative cells Dp = Dh/10 [47]. The comparison of the results of simulation based on this particular couple of diffusion coefficients Dp, Dh provided rather good agreement not only with MTS growth characteristics obtained in Ref. [68], but also described well the experimental density distribution on the MTS front.

2.3.4 Oxygen transport and consumption

Diffusion coefficient DO2 varies in sources within the range DO2 = 86.4–160 mm2/day [94], [, 66]. In the current work the range of its variation DO2 = 40.0–150 mm2/day has been studied and the value of DO2 = 86.4 mm2/day provided the best fit to the experimental data of Ref. [68].

Oxygen uptake αpfor Cp varies insignificantly 1.38–2.8 × 10−2 (mm3/cell) mg/(L*day) [47], [, 66]. The value αp = 1.4 × 10−2 (mm3/cell) mg/(L*day) has been used for the calculations, as it keeps proliferation rate close to experimentally observed values.

Oxygen uptake for Chαh = αp/5 was taken as typical approach [47], [, 95].

2.3.5 Proliferation and phenotype transition

Proliferation Rates vary, but insignificantly and typically for Cp the doubling time is τp = 1.0 day while for hypoxic cells τh = 2.0 day [41], [47], [96].

Phenotype switch time (proliferative to hypoxic) τph = 1 h, (hypoxic to proliferative) τhp = 96 h [47]. Phenotype switch time (hypoxic to necrotic) 3.1 × 10−2 h−1 = 0.74 day−1 [52]

Switch thresholds Cphand Chn

In spite of dramatically large scale of data variation for phenotype transition thresholds [97], [98], [99] their actual values did not have noticeable effect on the process dynamics. The basic threshold values of oxygen concentration used in simulations (Table 3) have been chosen simply as mostly used by other authors. The further calculations show that variations of thresholds influences mostly initial stage of spheroid growth, but has rather moderate effect on total cell profiles of 3–5 days grown MTS.

3 Calculation results

In order to study the process sensitivity to the model variation, some model elements and parameters have been varied around the basic model configuration. The basic approach includes the overlapping cell phenotype switch and assumption of constant ECM permeability, as per (2). The basic assumption for oxygen supply is closed domain, which means that no oxygen is being supplied through the domain border and corresponds to reflection boundary conditions (RBC). Another condition for permanent oxygen supply (PSBC) is applied in calculations for comparison purposes only. The influence of variable permeability, as a function of strain according to (4), has been studied as variation of the basic model. Also, the stiff phenotype switch approach that strictly follows GoG rule has been studied as a variation of GoG model. The major calculations and comparison with experimental data is carried out in the frame of the basic approach. The comparative analysis of the variations is presented with regard to the results of the basic model.

3.1 Formation of cell density profiles

The results obtained in calculation are generally in line with the experimental data of Stein and colleagues [68] for two cell lines U87WT and U87ΔEGFR. Based on the observed cell density distribution the authors identified two main MTS zones. These were the bulk part of the spheroid-MTS Core and the Invasive Zone – the peripheral area of the growing spheroid where cell density drops from typical core values down to zero. There were threshold cell densities introduced in Ref. [68] to define the separation border between the core and IZ. The outer boundary of IZ is defined as a cut off cell density level on the outer edge of the Invasive Zone. These levels were different for the two studied cell lines: CIZ1 = 104 cells/mm3 and CIZ2 = 2.5 × 103 cells/mm3 for U87WT and U87ΔEGFR correspondingly. The border between the core part and IZ was formally defined in Ref. [68] as a rim with cell densities CCR1 ≈ 105 cells/mm3 and CCR2 ≈ 8 × 104 cells/mm3 for the two studied cell lines.

The cores and IZs of the U87WT and U87ΔEGFR spheroids exhibit different invasiveness as well as different characteristic cell densities [68]. The current model based on average mechanobiological properties of glioma cells does not take into account the difference in individual characteristics of the two cell lines. Once the model is unable to quantify this difference we can speak here only about qualitative comparison of the simulations with the experimental results. Not carrying out a fitting to the data of any particular one of the two cell-lines studied in the experiment we present here simulation results obtained for the average parameters taken for glioma cells, according to Tables 1, 2, 3 and 4.

As the initial cell density is well below equilibrium value Cinit = Ceq/5, cell–cell adhesion represents dominating intercellular force in the whole spheroid. The distribution of cell–cell traction forces is presented on Figure 1A for t = 0.3 day. The intercellular stress curve exhibits non-monotonic tensile character, which is determined by local cell concentrations and FD function for glioma cells, according to (1a).

Figure 1: (A) Radial distribution of cell–cell traction forces σCell/P0 and cell concentrations CCell/Ceq for t = 0.3 day; Ceq = 3.5 × 105 cells/mm3, P0 = 100 Pa. (B) Proliferative (solid lines) and hypoxic (short-dashed lines) cell density distributions for the moments of time t = 0, 0.3 and 0.5 day. initial distribution (t = 0) represents proliferative cells only-rectangular area on the graph.

Figure 1:

(A) Radial distribution of cell–cell traction forces σCell/P0 and cell concentrations CCell/Ceq for t = 0.3 day; Ceq = 3.5 × 105 cells/mm3, P0 = 100 Pa. (B) Proliferative (solid lines) and hypoxic (short-dashed lines) cell density distributions for the moments of time t = 0, 0.3 and 0.5 day. initial distribution (t = 0) represents proliferative cells only-rectangular area on the graph.

In the beginning of MTS growth the inner traction forces cause contraction of MTS core part. Figures 1B and 2A illustrate this process as formation of high density spherical layer and its shift towards center, which ultimately leads to densification of MTS center within the first 12 h. The initial contraction was also observed in Stein’s et al. experiment as slight decrease of radius of the core spheroid.

Figure 2: (A) Distributions of total cell densities for t = 0.2–5.0 days; (B) distribution of oxygen concentrations for t = 0.2–5.0 days; The line ncr = ξh/n0 shows hypoxic transition threshold level; n0 is initial oxygen concentration, n0 = 4 mg/L..

Figure 2:

(A) Distributions of total cell densities for t = 0.2–5.0 days; (B) distribution of oxygen concentrations for t = 0.2–5.0 days; The line ncr = ξh/n0 shows hypoxic transition threshold level; n0 is initial oxygen concentration, n0 = 4 mg/L..

Two processes are forming the distribution of cell density immediately after the initial spheroid is placed into the hydrogel: spheroid contraction and separation of initially proliferative cells into different phenotype groups.

The hypoxic condition, being a triggering factor of phenotype transition in the model is achieved in the spheroid center during the first 8 h of its growth, when oxygen concentration drops below critical level ξh = 1.0 mg/L (Figure 2B). Time-limited diffusion of oxygen to the center of spheroid forms specific radial distributions of its concentration at different moments of time, which is one of the major factors determining radial cell density profile (Figure 1B and 2A). The oxygen profile also determines partial concentrations of different cell phenotypes in the total cell population. Figure 1B shows growing accumulation of hypoxic cells in the center and initial expansion of proliferative zone for 7.2 and 12 h. Once the phenotype transition began, the cell diffusion out of the core spheroid, which models cell motility, is growing up quickly and forms low density protrusion (Figure 3A). However, cell motility does not affect MTS expansion at its early stage, as the nearest external layer of the core spheroid is formed dominantly by proliferative cells and expands due to cell proliferation in the area of relatively high oxygen concentration in front of this layer. The phase velocity of proliferative layer expansion (Figure 3B) agrees well with in vitro Stein’s et al. data. The solid line on Figure 3B represents a rim of equilibrium cell density in the simulated MTS core where the intercellular traction is compensated by the repulsion forces. Above this isoline, the average intercellular distances exceed the equilibrium value, therefore traction forces dominate.

Figure 3: Isolines of cell concentrations obtained in simulation of the experimental conditions of Stein’s et al. [68]; (A) Solid black line corresponds to cell concentration of 400 cells/mm3, dashed line – 4 × 103 cells/mm3, dashed-dotted line is for 4 × 104 cells/mm3 and presents on both graphs (A) and (B); open circles show experimental data obtained by Stein et al. for U87WT glioma cell line (CIZ1 = 104 cells/mm3), open triangles – the results for U87ΔEGFR cells (CIZ2 = 2.5 103 cells/mm3); (B) solid black line denotes equilibrium cell density in the model Ceq = 4 × 105 cells/mm3, dotted line corresponds to the concentration of 1.6 × 105 cells/mm3, blue dashed – 8 × 104 cells/mm3 and dash-dotted – 4 × 104 cells/mm3; circles and triangles show the core boundary isolines CCR1 ≈ 105 cells/mm3 and CCR2 ≈ 8 × 104 cells/mm3 obtained in experiment for the two above mentioned glioma cell lines correspondingly.

Figure 3:

Isolines of cell concentrations obtained in simulation of the experimental conditions of Stein’s et al. [68]; (A) Solid black line corresponds to cell concentration of 400 cells/mm3, dashed line – 4 × 103 cells/mm3, dashed-dotted line is for 4 × 104 cells/mm3 and presents on both graphs (A) and (B); open circles show experimental data obtained by Stein et al. for U87WT glioma cell line (CIZ1 = 104 cells/mm3), open triangles – the results for U87ΔEGFR cells (CIZ2 = 2.5 103 cells/mm3); (B) solid black line denotes equilibrium cell density in the model Ceq = 4 × 105 cells/mm3, dotted line corresponds to the concentration of 1.6 × 105 cells/mm3, blue dashed – 8 × 104 cells/mm3 and dash-dotted – 4 × 104 cells/mm3; circles and triangles show the core boundary isolines CCR1 ≈ 105 cells/mm3 and CCR2 ≈ 8 × 104 cells/mm3 obtained in experiment for the two above mentioned glioma cell lines correspondingly.

The initial cell motion towards the spheroid center slows down while the central density grows up to Ceq. At the same time, the invasion of motile cells into the surrounding porous scaffold out of MTS core forms extended IZ within the first two days, as shown on Figure 2A and further on Figures 3, 4 and 5. These two processes: spheroid contraction with densification of its center and the fast development of low density IZ provide substantial redistribution of MTS density during the first day of its growth (Figures 1B and 2A). At further MTS growth, this elongated shape of cell density profile across its border does not change significantly.

Figure 4: Comparison the results of calculation with experimental data for day 3 of MTS growth obtained in [68] for the U87WT cell line (solid spheres). The figures present calculated distributions of total cell density obtained for DO2 = 86.4 mm2/day; Proliferative, hypoxic and necrotic cells in different scales for clear comparison (A) and (B).

Figure 4:

Comparison the results of calculation with experimental data for day 3 of MTS growth obtained in [68] for the U87WT cell line (solid spheres). The figures present calculated distributions of total cell density obtained for DO2 = 86.4 mm2/day; Proliferative, hypoxic and necrotic cells in different scales for clear comparison (A) and (B).

Figure 5: Comparison the results of calculation with experimental data for day 3 of MTS growth obtained in [68] for the U87WT cell line (solid spheres). The figures present calculated distributions of total cell density obtained for DO2 = 40.0 mm2/day; Proliferative, hypoxic and necrotic cells in different scales for clear comparison (A) and (B).

Figure 5:

Comparison the results of calculation with experimental data for day 3 of MTS growth obtained in [68] for the U87WT cell line (solid spheres). The figures present calculated distributions of total cell density obtained for DO2 = 40.0 mm2/day; Proliferative, hypoxic and necrotic cells in different scales for clear comparison (A) and (B).

The results of the simulation satisfy well the measured characteristic times of the core development, which exhibits rapid expansion and further saturation. IZ expansion curves also show rapid initial growth of IZ over 200 μm within the first day (Figure 3). Though, the cell densities do not correspond exactly to the levels measured in experiment.

The expansion of IZ saturates within five days (Figure 3A) and the growth of the MTS core also shows signs of saturation (Figure 3B), which is also in line with the experimental results. The two isolines of cell densities 4 × 102 and 4 × 103 cells/mm3 (Figure 3A) correspond well to time-dynamics of IZ borders observed by Stein and colleagues [68]. The distribution of cell density along the MTS radius is generally non-monotonic, see Figures 1A and 2B.

Non-monotonic growth character can be also seen on cell density isolines (Figure 3B) in the density range 8 × 104–1.6 × 105 cells/mm3. This is the range of strongest cell–cell adhesion and therefore intercellular traction forces. The traction forces build local cell agglomerations in IZ that grow fast at relatively high oxygen content in front of the spheroid. Similar local steps and elevations in the invasive zone of tumor can be seen on cell density profiles obtained in experiments (Figures 4 and 5). At certain conditions these agglomerations may grow fast suppressing cell growth around them through consumption of locally available oxygen and forming long protrusions and subtumors. These effects are considered further in this work.

As calculations demonstrated (Figures 4 and 5), cell density profiles are rather sensitive to oxygen diffusion, which does not fully replenish its consumed volumes inside MTS. The diffusion of oxygen determines to great extent the intensity of phenotype transition along the MTS radius and therefore cell density profile. The simulation demonstrates typical dynamics of MTS growth when the proliferative cells quickly turn into hypoxic phenotype in the central area (Figure 1B) and by day 3 the central part is mostly necrotic (Figure 4A). The proliferative cells are confined within a thin layer in the outer part of MTS (Figure 4B at R ≈ 0.48 mm and Figure 5B at R ≈ 0.43 mm). The figure also shows sharp vertical edge of cell density profile observed in experiment and tenfold rise of cell concentration in the center compared to its initial value. The rise cannot be explained by local cell proliferation only. This is presumably a result of compression from expanded ECM scaffold and cell–cell adhesion effect in the initial spheroid. The effect of spheroid contraction is also seen on simulated profiles, but it is not that strong as in experiment. The simulation data shows about sevenfold increase of cells density in the center. There are set of parameters responsible for this rise, such as ECM stiffness and permeability, diffusion of oxygen and barriers of cell phenotype transition. A detailed analysis of their combined effect on maximum cell density has not been carried out. However, a variation of any single parameter did not provide better agreement with the experimental results. The simulated cell density in the very center was always only 60% of that obtained in experiment.

The sensitivity of MTS growth to oxygen diffusion coefficient is illustrated by Figure 5. The diffusion coefficient has been reduced down to DO2 = 40.0 mm2/day, which is less than a half of that reported in literature, see Table 3. This relatively low diffusion coefficient sharpens the edge of the core spheroid (Figure 5A) and slows down its growth. This moderate effect of oxygen diffusion on the cell density profile is well expected and predictable.

Generally, the rise of total cell density in the middle of IZ obtained in simulation is conditioned by thin proliferative rim (Figures 4B and 5B) located outside of the core spheroid at R ≈ 0.43–0.48 mm. The profile measured in experiment for U87WT cell line also exhibits a clear step around the same point of R ≈ 0.43 mm, which is presumably a sign of proliferative rim.

The measured by Stein et al. density distribution of three days MTS for U87WT cells shows another rise of cell density at the end of IZ at R ≈ 0.85 mm, which was not exactly reproduced in simulation.

The simulated Stein’s conditions generally result in non-monotonic radial profiles of cell concentrations (Figure 2A). This situation may theoretically lead to formation of distant cells clusters in the form of additional proliferating rims. Similar scenario was identified by other researchers as formation of subtumors [55]. Cell–cell adhesion–repulsion has been reported as major mechanism responsible for this effect. We also carried out simulation of subtumors and consider conditions for their formation further in this work.

Though the calculated core densities are generally lower in simulation, the comparison with Stein’s et al. data demonstrates qualitative agreement. The calculated cell density profiles exhibit similar dense core with a steep density fall at its edge and elongated IZ presenting a shelf with steps or small distant rise of cell density (Figures 4 and 5). In so doing, the comparison shows better agreement with the data obtained for U87WT cell line. This cell line exhibits relatively smooth IZ, which should therefore closer correspond to our 1D simulation of radial MTS growth than clearly discrete protrusions in IZ obtained by Stein for the another cell line U87ΔEGFR.

It should be noted that all calculations carried out for comparison with the experimental data have been prepared at RBC, as a basic case. These boundary conditions assume closed domain with limited total oxygen content inside it. MTS growth at RBC inevitably saturates and the saturation time is obviously monotonic function of the domain volume. In so doing, the main proliferation rim exhausts completely. In the absence of proliferative cells MTS expansion slows down considerably, while the growth of IZ is observed within 1–2 extra days width decreasing speed (Figure 3). This leads to further elongation of cell density profile, which produces an effect of its slow expansion.

In case of permanent supply of oxygen through the domain border (PSBC case) the proliferative layer is more stable and oxygen consumption is limited by its diffusion only. This case has also been simulated in the current work.

3.2 The effect of stress distribution and variable permeability

The expansion of MTS resulting from cell proliferation and their intrusion into the surrounding porous scaffold creates significant tangential tension (Figures 6 and 7). The tension value exhibits a sharp raise around the spheroid center. The calculated values of ECM tension are based on the assumption of pure hydrogel elasticity, which works well enough in the whole domain except for the very center of the spheroid. The raise of tension in MTS center is obviously limited by hydrogel yield point and should lead at this point to partial stress relaxation in the center. In the model we do not consider ECM plasticity simply keeping the high values of stress in the spheroid center. In calculation this high ECM stress in the middle does not have any noticeable effect on MTS dynamics, as the partial concentration of ECM fraction is rapidly dropping down towards the spheroid center.

Figure 6: Stress distribution along MTS radius at day 3. (A) The total stress distribution σtot/σ0 is calculated within iso-strain model; the graph also shows radial component of the total stress σRtot/σ0, the stress conditioned by cells interaction σCell/P0 and concentration of cells Ccell/Ceq, σ0 = 100 Ps; (B) distribution of ECM stress σECM/σ0, tangential component of ECM stress στ/σ0 and its radial component σR/σ0.

Figure 6:

Stress distribution along MTS radius at day 3. (A) The total stress distribution σtot/σ0 is calculated within iso-strain model; the graph also shows radial component of the total stress σRtot/σ0, the stress conditioned by cells interaction σCell/P0 and concentration of cells Ccell/Ceq, σ0 = 100 Ps; (B) distribution of ECM stress σECM/σ0, tangential component of ECM stress στ/σ0 and its radial component σR/σ0.

Figure 7: Stress distribution along MTS radius at day 5. (A) The total stress distribution σtot/σ0 is calculated within iso-strain model; the graph also shows radial component of the total stress σRtot/σ0, the stress conditioned by cells interaction σCell/σ0 and concentration of cells Ccell/Ceq, σ0 = 100 Ps; (B) distribution of ECM stress σECM/σ0, tangential component of ECM stress στ/σ0 and its radial component σR/σ0.

Figure 7:

Stress distribution along MTS radius at day 5. (A) The total stress distribution σtot/σ0 is calculated within iso-strain model; the graph also shows radial component of the total stress σRtot/σ0, the stress conditioned by cells interaction σCell/σ0 and concentration of cells Ccell/Ceq, σ0 = 100 Ps; (B) distribution of ECM stress σECM/σ0, tangential component of ECM stress στ/σ0 and its radial component σR/σ0.

In the outward direction along the spheroid radius the radial component of total stress in the cell-matrix composite exhibits a sharp transition from tensile to compressive stress (Figures 6A and 7A). This transition is conditioned to a great extent by ECM reaction on the expansion of MTS, forming similar stress distribution in the matrix within relatively thin near-MTS-surface layer (Figures 6B and 7B). The pressure exerted on the matrix from the expanding MTS creates this transition layer (TL), as described above, with relatively high radial compression of ECM and tangential tension at the same time (Figures 6B and 7B).

The expansion of cell aggregation behind TL causes both radial and tangential stretching of the extracellular scaffold in the bulk spheroid (Figures 6B and 7B). The total stress in the system is tensile and exhibits nearly constant value inside the core, see the negative level shelf on Figures 6A and 7A and a sharp step up within the TL. In general, even simple model of a single elastic incompressible continuum [33] shows a tension zone in the center of expanding tumor spheroid and radial compression on its periphery. But that simple approach does not produce the fine structure of stress distributions.

The constant level of the total stress in the composite system is being formed through mechanical interaction between its two constituents – cells and ECM. The hydraulic pressure built by proliferating cells is decreasing along the MTS radius following the decrease of total cell concentration, which is still above Ceq and therefore causes compressive stress inside the spheroid behind TL. In the transition layer on spheroid surface the cell concentration drops below Ceq creating certain tensile forces in cell aggregation due to dominating intercellular adhesion (Figures 6A and 7A). The positive intercellular pressure produces compensatory effect against tensile stress in the ECM scaffold inside the spheroid. This effect is illustrated by two curves σC (R) and σECM (R) demonstrating opposite sign stresses within the two constituents of the composite system (Figures 6B and 7B). As a result of this effect the total stress calculated within the isostrain composite model is nearly constant throughout major part of the spheroid σtot (Figure 7A). This constant value of stress inside MTS is being maintained by self-consistent process that involves interaction between cells and ECM, resulting in redistribution of stresses and cellular drift. This drift modulates the distribution of cell concentration in accordance with the total stress gradients inside the spheroid.

The sharp raise of the total stress across the transition layer builds hydrostatic pressure gradients within this narrow area directed outwards across the TL. This layer represents a barrier for diffusive cell transport out of MTS core. The resulting pressure force steepens the front of cell density distribution, pushing back some of the motile cells traveling out of the core. Nevertheless, as simulations demonstrate, the highly motile cells typically overcome this pressure barrier in cases of soft and moderately rigid ECMs.

Apart from the particular effect of stress on ECM pore size, the matrix structure and its mechanical reordering are not simulated in the current work. The effect of stress on ECM structure is considered here in the frame of simple model and may depend on the matrix type.

The specific stress distribution within TL obtained in simulation may presumably affect the structure of ECM scaffold. In case of randomly arrayed fibers they should expectedly be rearranging in front of TL, taking directions normal to MTS radius. The fibers should undergo further reorientation within TL along spheroid radius, as radial stress components turn tensile inside the spheroid. This scenario of mechanical remodeling of collagen scaffold has been observed in vitro [29] on murine colon carcinoma. The authors have shown the tangential alignment of the collagen fibers around the spheroid surface and general radial tension inside it. These general traction forces (Figures 6A and 7A) are assumed to facilitate tumor invasion providing tracks through mechanical alignment in matrix structure [28], [, 29]. The ability of collagens to remodel and in particular to align under external strain forces applied has been demonstrated in experiments [100] as a typical property of fibrous biological networks. This effect of mechanical remodeling, in the frame of the current model should technically result in local increase of ECM permeability. However, the observed phenomenology is still not studied well enough to quantify its influence within the model. We can only note that the stress distribution obtained in simulation may practically condition the observed experimentally mechanical restructure in ECM.

It must be stressed that there is a principal difference between the current two-constituents’ composite model and a single-constituent only cellular models. In single constituent models, where only hydrostatic pressure inside the proliferating cell colony determines its expansion [68], [69], [73], the elastic reaction of extracellular constituent on this expansion is often ignored. A regime of stationary MTS growth assumes no intercellular pressure σc ≈ 0, keeping cell concentrations around equilibrium value Ceq within a single constituent self-consistent model [69].

In contrast to these simple cases the present model takes into account the mutual interaction of the two compressible constituents of the composite system: visco-elastic cell aggregation and elastic matrix within the whole simulation domain. This interaction forms a specific cell density profile inside the spheroid and a thin transition layer on its surface, where most mechanical characteristics of the composite system undergo rapid transition.

3.3 Formation of subtumors

The next series of calculations has been carried out at high ECM stiffness EECM = 2.5 kPa and permanent oxygen supply through the domain border (PSBC). These specific parameters conditioned formation of subtumor. There are no laboratory data for comparison and the simulations have been carried out in order to demonstrate this effect as potential scenario discussed earlier in Refs. [55], [56], [57]. The matrix stiffness is selected so that provides substantial contraction of the initial spheroid and sharp front at its edge without elongation of cell density profile into peripheral zone typical for soft matrixes, as on Figures 4 and 5. Instead, the poor accumulation of the invading cells out of the core spheroid shrinks to a narrow densified zone due to cell–cell adhesion and the matrix tension, forming another proliferative rim (Figure 8A). By day two shown on the figure, the major part of the first proliferative rim is located in the area of oxygen shortage and undergoes cell phenotype transition, while the small frontal proliferative rim is formed in the area of high oxygen content. By day three over 90% of the main proliferative rim turns hypoxic, as the oxygen level inside the core part drops below the transition threshold (Figure 8C).

Figure 8: Formation of subtumor. Calculation results presented for days 2, 3, 5 of MTS growth at elevated matrix stiffness EECM = 2.5 kPa and PSBC. The graphs show radial cell density distributions (A, C, D): for proliferative cells Cp, hypoxic Cm, necrotic Cn and total cell density Ctot; also the distribution of oxygen concentration is shown on these graphs; the dashed line ncr = ξh/n0 shows the threshold level of cell phenotype transition proliferative-to-hypoxic, n0 = 4 mg/L; (B) shows stress distribution for day 2 including total stress in cells–ECM composite σtot, intercellular stress distribution σcell and radial component of the total stress σRtot.

Figure 8:

Formation of subtumor. Calculation results presented for days 2, 3, 5 of MTS growth at elevated matrix stiffness EECM = 2.5 kPa and PSBC. The graphs show radial cell density distributions (A, C, D): for proliferative cells Cp, hypoxic Cm, necrotic Cn and total cell density Ctot; also the distribution of oxygen concentration is shown on these graphs; the dashed line ncr = ξh/n0 shows the threshold level of cell phenotype transition proliferative-to-hypoxic, n0 = 4 mg/L; (B) shows stress distribution for day 2 including total stress in cells–ECM composite σtot, intercellular stress distribution σcell and radial component of the total stress σRtot.

A radial component of the total stress in cell–ECM system exhibits strong gradient on the inner side of the second proliferative rim (Figure 8B). The peak magnitude of the radial stress component in front of the core spheroid σRtot ∼ σ0 is nearly 10 times higher than that formed at basic ECM stiffness (compare with Figure 6). The resulting radial force pushes all cells within this spherical layer in outer direction, so that it covers a distance of ∼0.1 mm within a day (Figure 8A and c). Cell–cell adhesion and rigid matrix keep this layer compact and rather thin. After three days of MTS growth the mechanics of displacement of this distant proliferative rim is being gradually changed. The inner force, pushing it out is decreasing. At the same time, the oxygen level behind the rim goes down below the threshold of cell phenotype transition (Figure 8D). As a result, only outer cells of this layer retain proliferative type. The other cells are turning hypoxic. This process creates a growing trace of hypoxic cells behind the rim (compare Figure 8C and D). The cell density in the proliferative layer is decreasing and its expansion is conditioned to great extent by the level of oxygen at its outer side, sustaining actively proliferating regime. The speed of the proliferative layer expansion is determined now by oxygen consumption rate, as the layer follows the shrinking outer oxygen zone with its concentration above the threshold level (Figure 8D). Certain amount of hypoxic cells goes ahead of the proliferative layer due to their high motility, but they still do not contribute into the layer displacement, as the reverse phenotype transition hypoxic- to -proliferative requires rater long residence time (≥96 h) in the area of high oxygen concentration (see Table 4).

It’s important to note that the simulated scenario does not have an immediate practical meaning. The distant subtumor layer has been observed in simulation within rather narrow range of ECM stiffness values at specific oxygen concentration profiles and boundary conditions at the domain borders. However, all these parameters, including oxygen diffusion coefficients, cell motility, are feasible and typical for glioma MTSs (see Tables 1, 2, 3 and 4). The observed subtumor layer is formed at the main tumor interface and grown away from the core spheroid. This scenario is generally consistent with earlier obtained results of in silico models and biological experiments [55], [, 101] demonstrates that the subpopulations of tumor cells are being formed by multiple cells of the main tumor surface and not as a result of hyper-proliferation of a distant single cell. However, the models exploited in Refs. [55], [, 101] do not take into account viscoelastic interaction between the two compressible fractions cells and ECM. The influence of microenvironment on cell dynamics is considered without mutual rheological affects. The formation of cell subpopulations has been simulated in earlier models as a result of circumferential morphological instability coming from the interplay between oxygen transport, cell phenotype transitions and cell adhesion. The current model demonstrates that the same processes together with effects of compressibility and mutual elastic influence in cell–ECM system may lead directly to radial discretization of invasive zone and formation of separate cell populations.

It is important to stress that permanent oxygen supply through the domain border (PSBC) plays a pivotal role in formation of subtumors. The basic boundary conditions for oxygen supply RBC, when the total amount of oxygen in the domain is limited to constant value, do not lead to formation of distant cell agglomerations out of the core spheroid (Figure 9). The calculation results for RBC and the same as above ECM stiffness EECM = 2.5 kPa show a sharp edge of the tumor spheroid and no invasive zone or distant cell clusters. By day 5 the proliferation rim disappears completely, as shown on Figure 9, due to general shortage of oxygen in the domain and its low level inside the spheroid in particular.

Figure 9: Calculations carried out for the same conditions as on Figure 8, except for limited oxygen volume in the domain. ECM stiffness EECM = 2.5 kPa, RBC for oxygen supply. Radial cell density distributions are presented for day 5: proliferative cells Cp, hypoxic Cm, necrotic Cn, total cell density Ctot; and distribution of oxygen concentration.

Figure 9:

Calculations carried out for the same conditions as on Figure 8, except for limited oxygen volume in the domain. ECM stiffness EECM = 2.5 kPa, RBC for oxygen supply. Radial cell density distributions are presented for day 5: proliferative cells Cp, hypoxic Cm, necrotic Cn, total cell density Ctot; and distribution of oxygen concentration.

4 Discussion

A two component in silico model of glioma MTS growth, presented in this work takes into account rheological properties of cell–ECM system, as well as the effect of mechanical ECM deformation. Cell colony, as a separate component in the model represents a porous compressible aggregation. The model considers interaction between the two compressible constituents of the composite system, visco-elastic cell aggregation and elastic matrix, within the whole simulation domain.

The transitions between the three cell phenotypes, proliferative, hypoxic and necrotic, were triggered by the local oxygen content, which was the only metabolite considered in the model. It was demonstrated that diffusion of oxygen from the peripheral zone towards the spheroid center and its distributed consumption built a specific profile of oxygen content, which affects in a self-consistent loop the cell density profile and ultimately the character of tumor invasion.

Special consideration has been given to the case of reduced diffusion of oxygen, which sharpens spheroid edge, bringing the result more in line with the observed shapes, but generally reduces cell concentrations in the core part below the experimental values.

The simulated cell density profiles are in good quantitative agreement with these obtained in experiment within the major part of the spheroid, except for the spheroid center, where the simulated density is only 60% of the measured value. Cell densities measured in the center cannot be explained by local cell proliferation only, being well above the figures estimated in that simple way. Another factor affecting this magnitude is a densification of the center through its compression.

The initial contraction observed experimentally, as slight decrease of spheroid radius has been also obtained in simulation. As long as the starting cell concentration in the spheroid is below equilibrium value, so that cell–cell adhesion overcomes repulsive intercellular forces, the model demonstrates certain shrinkage of the initial tumor volume.

Further compression of cell aggregation in the central part of MTS due to the pressure exerted by expanded ECM has also been simulated. Even though, the simulated compression does not provide sufficient densification of the spheroid center to correspond exactly to the measured values.

An artificial increase of ECM stiffness above the values typical for collagens used in MTS growth tests elevates the calculated cell densities in the center much closer to the measured value, but still not high enough to match it exactly. The model presumably underestimates compressive forces, acting within the spheroid or compressibility of the given cellular continuum. These results demonstrate however, an important role of stresses and compressibility of cell-matrix composite system in MTS growth.

A relatively thin transition layer forms on the spheroid surface, where mechanical characteristics of cell–ECM composite undergo rapid transition. This layer is an important part of MTS structure and results from elastic ECM reaction on spheroid expansion. This structure can be simulated only within the two component model that takes into account the mutual influence of the two constituents of the composite system. Their interaction forms a specific cell density profile inside the spheroid and the Transition Layer on its surface. In particular, there is a nearly constant level of the total stress in the MTS core behind the transition layer. The flat level of tensile stress results from compensatory mechanical interaction between cell aggregation, which forms compressive stress and the stretched ECM with growing tensile stress. Outside the core spheroid, within transition layer ECM exhibits rapid relaxation of tensile stress and cell concentration quickly drops down resulting in sharp change of intercellular stress from compressive to tensile character due to cell–cell adhesion.

The simulations demonstrated that sharp stress relaxation within the transition layer should lead to a reordering of ECM structure. In the case of MTS growth in collagen, the fibers are expected to reorder in front of the Layer aligning to the spheroid surface, undergoing tensile tangential and compressive radial forces. Another reorientation of fibers is expected inside the spheroid, along the radius. This type of fiber reordering has been observed in MTS experiments [29] as mechanical remodeling of ECM. The matrix remodeling should have a modulatory effect on collagen permeability and therefore cell invasiveness. However, there is no clear quantitative model to take it into account.

Another effect of mechanical ECM remodeling through matrix deformation has been studied quantitatively in the current simulations. It was reasonably assumed that the cell population growing into ECM scaffold causes its expansion, stretching pores and therefore affecting its permeability and further-mobility of invasive cells. The average effect of this matrix deformation resulted in increase of tumor invasion speed up to 20%. This result demonstrates the importance of mechanical effects on ECM structure and their possible influence on cell invasiveness.

Rather brief consideration has been given in the work to the influence of particular physical factors, such as cell–cell adhesion, rheological properties of the matrix and cell aggregation, oxygen supply and diffusion rate on glioma mechanobiology. Each of these factors deserves closer study including in silico simulations of microscale scenarios.

The simulation results show sensitivity to the type of GoG model. The stiff GoG model, which does not allow any proliferation of hypoxic cells and vice versa, demonstrates relatively low cell density profile and steep edge of the MTS core. The overlapped GoG model shows overall higher cell density, more elongated profile and better agreement with the experimental results.

The simulated glioma spheroids exhibit higher cell density and faster expansion in rigid matrixes that is in line with in vitro studies of glioma MTSs [6], [, 19]. The model shows that an increase of ECM stiffness sharpens the edge of MTS, which tends to form compact structures in rigid matrixes. Wherein, the cell density profile is getting shorter and steeper due to decrease of invasive zone. The shrinkage of the invasive zone in rigid ECM is conditioned by high potential barrier for motile cells to travel across transition layer, against the background of large stress gradients across the surface of the core spheroid.

Two scenarios have been observed in simulations for rigid matrixes: either a single compact spheroid with sharp edge is being formed or the single spheroid is accompanied by additional fragmented cell structures instead of elongated continuous invasive zone of low cell density typical for soft matrixes. The simulations indicate a presence of specific range of ECM stiffness and nutrient concentrations when additional proliferative rim may grow away from the main one. A specific stress distribution forms hydrostatic radial pushing forces at the edge of MTS core spheroid. So, that the second proliferative rim is being pushed further away. This result demonstrates that cell motility is probably not the key factor in formation of distant subtumor clusters, as it would play an important role only in the scenario when a single highly motile cell could invade far into nutrient-rich area and hyper-proliferate there. This scenario did not take place in the current study. The further motion of the distant proliferative layer is provided by active proliferation in the outer zone reach of nutrients (oxygen in the current model). Cell–cell adhesion and ECM stiffness play crucial role in formation of these distant multicellular proliferative clusters.

Over all, three scenarios of glioma invasion at different ECM stiffness values, driven by mere mechanical principles of continuum and iso-strain dynamics of cell–ECM composite, have been shown in simulations. MTS growing in soft matrix exhibits gentle cell density profile with elongated low density invasive zone, while rigid ECM conditions sharp edge of MTS and higher invasion speed. Moreover, the highest speed of tumor invasion is obtained at intermediate values of ECM stiffness, where a compact cell clusters are being formed instead of invasive zone as additional proliferative layers out of the core spheroid.

It is important to mention that the current 1D model simulates spherically symmetric process of tumor growth. In a real 3D MTS growth, elongated protrusions in the form of fingers should form due to morphological instability of tumor border [55], [, 102], generating azimuthal discretization of cellular continuum.

In subtumor conditions, radial discretization of cell density profile is added to the azimuthal one. This case should exhibit random array of proliferating islands separated from the main tumor spheroid, instead of distant spherical layers obtained in 1D simulations, as the uniform layers must fall off due to azimuthal instability. This characteristic dispersed distribution of distant proliferative islands has been observed in experiments [32], [, 68]. These patterns were identified as ‘starburst’ type distributions, typical for aggressive glioma cell lines.

The scenarios considered in the current study are based completely on relatively simple mechanics of elastic compressible composite system. It demonstrates scenarios of tumor invasion that can be used as basic approximation in further analysis of avascular stage of glioma including its signaling context.

Quantitative calculations of the effects of mechanical ECM restructure was currently limited to alteration of pore sizes only, under stress imposed. Detailed simulations of matrix structure and its remodeling effects, which play substantial role in tumor progression, require further evolution of the model.

Another important factor of MTS growth not included into the model, that may have an appreciable effect on mechanobiological aspects of glioma growth is the removal of necrotic cells from the core spheroid. It may particularly affect cells and stress distributions inside the spheroid within relatively long growth times and is worth including into the model in its further development.


1. De Vleeschouwer, S, editor. Glioblastoma. Brisbane: Codon Publications; 2017. https://doi.org/10.15586/codon.glioblastoma.2017. Search in Google Scholar

2. Martirosyan, NL, Rutter, EM, Ramey, WL, Kostelich, E, Kuang, Y, Preul, MC. Mathematically modeling the biological properties of gliomas: a review. Math Biosci Eng 2015;12:879–905. https://doi.org/10.3934/mbe.2015.12.879. Search in Google Scholar

3. Alfonso, J, Talkenberger, K, Seifert, M, Klink, B, Hawkins-Daarud, A, Swanson, KR, et al.. The biology and mathematical modelling of glioma invasion: a review. J R Soc Interface 2017;14:20170490. https://doi.org/10.1098/rsif.2017.0490. Search in Google Scholar

4. Rape, A, Ananthanarayanan, B, Kumar, S. Engineering strategies to mimic the glioblastoma microenvironment. Adv Drug Deliv Rev 2014;79–80:172–83. https://doi.org/10.1016/j.addr.2014.08.012. Search in Google Scholar

5. Grundy, TJ, De Leon, E, Griffin, KR, Stringer, BW, Day, BW, Fabry, B, et al.. Differential response of patient-derived primary glioblastoma cells to environmental stiffness. Sci Rep 2016;6:23353. https://doi.org/10.1038/srep23353. Search in Google Scholar

6. Ananthanarayanan, B, Kim, Y, Kumar, S. Elucidating the mechanobiology of malignant brain tumors using a brain matrix-mimetic hyaluronic acid hydrogel platform. Biomaterials 2011;32:7913–23. https://doi.org/10.1016/j.biomaterials.2011.07.005. Search in Google Scholar

7. Pathak, A, Kumar, S. Independent regulation of tumor cell migration by matrix stiffness and confinement. Proc Natl Acad Sci USA 2012;109:10334–9. https://doi.org/10.1073/pnas.1118073109. Search in Google Scholar

8. Rubenstein, BM, Kaufman, LK. The role of extracellular matrix in glioma invasion: a cellular Potts model approach. Biophys J 2008;95:5661–80. https://doi.org/10.1529/biophysj.108.140624. Search in Google Scholar

9. Umesh, V, Rape, AD, Ulrich, TA, Kumar, S. Microenvironmental stiffness enhances glioma cell proliferation by stimulating epidermal growth factor receptor signaling. PLoS ONE 2014;9:7. https://doi.org/10.1371/journal.pone.0101771. Search in Google Scholar

10. Yang, Y, Motte, S, Kaufman, LJ. Pore size variable type I collagen gels and their interaction with glioma cells. Biomaterials 2010;31:5678–88. https://doi.org/10.1016/j.biomaterials.2010.03.039. Search in Google Scholar

11. Beadle, C, Assanah, MC, Monzo, P, Vallee, R, Rosenfeld, SS, Canoll, P. The role of Myosin II in glioma invasion of the brain. Mol Biol Cell 2008;19:3357–68. https://doi.org/10.1186/1741-7007-10-29. Search in Google Scholar

12. Thorne, RG, Nicholson, C. In vivo diffusion analysis with quantum dots and dextrans predicts the width of brain extracellular space. Proc Natl Acad Sci USA 2006;103:5567–72. https://doi.org/10.1073/pnas.0509425103. Search in Google Scholar

13. Fraley SI, I., Wu, P, He, L, Feng, Y, Krisnamurthy, R, Longmore, GD, et al.. Three-dimensional matrix fiber alignment modulates cell migration and MT1-MMP utility by spatially and temporally directing protrusions. Sci Rep 2015;5:14580. https://doi.org/10.1038/srep14580. Search in Google Scholar

14. Xie, J, Bao, M, Bruekers, SMC, Huck, W TS. Collagen gels with different fibrillar microarchitectures elicit different cellular responses. ACS Appl Mater Interfaces 2017;9:19630–7. https://doi.org/10.1021/acsami.7b03883 . Search in Google Scholar

15. Kaufman, LJ, Brangwynne, CP, Kasza, KE, Filippidi, E, Gordon, VD, Deisboeck, TS, et al.. Glioma expansion in collagen I matrices: analyzing collagen concentration-dependent growth and motility patterns. Biophys J 2005;89:635–50. https://doi.org/10.1529/biophysj.105.061994. Search in Google Scholar

16. Rubenstein, BM, Kaufman, LJ. The role of extracellular matrix in glioma invasion: a cellular Potts model approach. Biophys J 2008;95:5661–80. https://doi.org/10.1529/biophysj.108.140624. Search in Google Scholar

17. Wang, C, Tong, X, Yang, F. Bioengineered 3D brain tumor model to elucidate the effects of matrix stiffness on glioblastoma cell behavior using PEG-based hydrogels. Mol Pharm 2014;11:2115–25. https://doi.org/10.1021/mp5000828. Search in Google Scholar

18. Ulrich, TA, Jain, A, Tanner, K, MacKay, JL, Kumar, S. Probing cellular mechanobiology in three-dimensional culture with collagen-agarose matrices. Biomaterials 2010;31:1875–84. https://doi.org/10.1016/j.biomaterials.2009.10.047. Search in Google Scholar

19. Ulrich, TA, de Juan Pardo, EM, Kumar, S. The mechanical rigidity of the extracellular matrix regulates the structure, motility, and proliferation of glioma cells. Cancer Res 2009;69:4167–74. https://doi.org/10.1158/0008-5472.CAN-08-4859-74. Search in Google Scholar

20. Payne, LS, Huang, PH. Review: the pathobiology of collagens in glioma. Molecular Cancer Research; 2013. Available from: https://mcr.aacrjournals.org/content/early/2013/07/16/1541-7786.MCR-13-0236 [MCR-13-0236]. https://doi.org/10.1158/1541-7786. Search in Google Scholar

21. Louca, M, Stylianou, A, Minia, A, Pliaka, V, Alexopoulos, LG, Gkretsi, V, et al.. Ras suppressor-1 (RSU-1) promotes cell invasion in aggressive glioma cells and inhibits it in non-aggressive cells through STAT6 phospho-regulation. Sci Rep 2019;9:7782. https://doi.org/10.1038/s41598-019-44200-8. Search in Google Scholar

22. Hofmann, M, Guschel, M, Bernd, A, Bereiter-Hahn, J, Kaufmann, R, Tandi, C, et al.. Lowering of tumor interstitial fluid pressure reduces tumor cell proliferation in a xenograft tumor model. Neoplasia 2006;8:89–95. https://doi.org/10.1593/neo.05469. Search in Google Scholar

23. Tse, JM, Cheng, G, Tyrrell, JA, Wilcox-Adelman, SA, Boucher, Y, Jain, RK, et al.. Mechanical compression drives cancer cells toward invasive phenotype. Proc Natl Acad Sci USA 2012;109:911–16. https://doi.org/10.1073/pnas.1118910109. Search in Google Scholar

24. Paszek, MJ, Zahir, N, Johnson, KR, Lakins, JN, Rozenberg, GI, Gefen, A, et al.. Tensional homeostasis and the malignant phenotype. Canc Cell 2005;8:241–54. doi:https://doi.org/10.1016/j.ccr.2005.08.010. Search in Google Scholar

25. Cheng, G, Tse, J, Jain, RK, Munn, LL. Micro-environmental mechanical stress controls tumor spheroid size and morphology by suppressing proliferation and inducing apoptosis in cancer cells. PLoS One 2009;4:e4632. https://doi.org/10.1371/journal.pone.0004632. Search in Google Scholar

26. Desmaison, A, Frongia, C, Grenier, K, Ducommun, B, Lobjois, V. Mechanical stress impairs mitosis progression in MultiCellular tumor spheroids. PLoS ONE 2013;8:e80447. https://doi.org/10.1371/journal.pone.0080447. Search in Google Scholar

27. Enmon, RM, O’Connor, KC, Lacks, DJ, Schwartz, DK, Dotson, RS. Dynamics of spheroid self-assembly in liquid-overlay culture of DU 145 human prostate cancer cells. Biotechnol Bioeng 2001;72:579–91. https://doi.org/10.1002/1097-0290(20010320)72:6<579::AID-BIT1023>3.0.CO;2-L. Search in Google Scholar

28. Gordon, VD, Valentine, MT, Gardel, ML, Andor-Ardó, D, Dennison, S, Bogdanov, A. Measuring the mechanical stress induced by an expanding multicellular tumor system: a case study. Exp Cell Res 2003;289:58–66. https://doi.org/10.1016/S0014-4827(03)00256-8. Search in Google Scholar

29. Kopanska, KS, Alcheikh, Y, Staneva, R, Vignjevic, D, Betz, T. Tensile forces originating from cancer spheroids facilitate tumor invasion. PLoS ONE 2016;11:e0156442. https://doi.org/10.1371/journal.pone.0156442. Search in Google Scholar

30. Meshel, AS, Nackashi, DP, Franzon, PD, Sheetz, MP. Single force measurements on 3D collagen matrix. Biophys J 2002;82:411a. https://doi.org/10.1529/biophysj.105.061994. Search in Google Scholar

31. Vinci, M, Box, C, Eccles, SA. Three-dimensional (3D) tumor spheroid invasion assay. J Vis Exp 2015;99:e52686. https://doi.org/10.3791/52686. Search in Google Scholar

32. Vinci, M, Gowan, S, Boxall, F, Patterson, L, Zimmermann, M, Court, W, et al.. Advances in establishment and analysis of three-dimensional tumor spheroid-based functional assays for target validation and drug evaluation. BMC Biol 2012;10:29. https://doi.org/10.1186/1741-7007-10-29. Search in Google Scholar

33. Jones, AF, Byrne, HM, Gibson, JS, Dold, JW. A mathematical model of the stress induced during avascular tumour growth. J Math Biol 2000;40:473–99. https://doi.org/10.1007/s002850000033. Search in Google Scholar

34. Breward, CJ, Byrne, HM, Lewis, CE. The role of cell–cell interactions in a two-phase model for avascular tumour growth. J Math Biol 2002;45:125–52. Search in Google Scholar

35. Ngwa, M, Agyingi, E. Effect of an external medium on tumor growth-induced stress matthias. IAENG Int J Appl Math 2012;42:4. Search in Google Scholar

36. Ramírez-Torres, A, Rodríguez-Ramos, R, Merodio, J, Bravo-Castillero, J, Guinovart-Díaz, R. Mathematical modeling of anisotropic avascular tumor growth. Mech Res Commun 2015;69:8–14. Search in Google Scholar

37. Sciumè, G, Santagiuliana, R, Ferrari, M, Decuzzi, P, Schrefler, BA. A tumor growth model with deformable ECM. Phys Biol 2015;11:065004. https://doi.org/10.1088/1478-3975/11/6/065004. Search in Google Scholar

38. Roose, T, Netti, PA, Munn, LL, Boucher, Y, Jain, RK. Solid stress generated by spheroid growth estimated using a linear poroelasticity model. Microvasc Res 2003;66:204–12. https://doi.org/10.1016/S0026-2862(03)00057-8. Search in Google Scholar

39. Helmlinger, G, Netti, PA, Lichtenbeld, HC, Melder, RJ, Jain, RK. Solid stress inhibits the growth of multicellular tumor spheroids. Nat Biotechnol 1997;15:778–83. Search in Google Scholar

40. Giese, A, Loo, MA, Tran, N, Haskett, D, Coons, SW, Berens, ME. Dichotomy of astrocytoma migration and proliferation. Int J Canc 1996;67:275–82. Search in Google Scholar

41. Giese, R, Bjerkvig, R, Berens, ME, Westphal, M. Cost of migration: invasion of malignant gliomas and implications for treatment. J Clin Oncol 2003;21:1624–36. https://doi.org/10.1200/JCO.2003.05.063. Search in Google Scholar

42. Godlewski, J, Nowicki, M, Bronisz, A, Nuovo, G, Palatini, J, De Lay, M, et al.. MicroRNA-451 regulates LKB1/AMPK signaling and allows adaptation to metabolic stress in glioma cells. Mol Cell 2010;37:620–32. https://doi.org/10.1016/j.molcel.2010.02.018. Search in Google Scholar

43. Hatzikirou, H, Basanta, D, Simon, M, Schaller, K, Deutsch, A. ‘Go or grow’: the key to the emergence of invasion in tumour progression? Math Med Biol 2012;29:49–65. https://doi.org/10.1093/imammb/dqq011. Search in Google Scholar

44. Böttger, K, Hatzikirou, H, Chauviere, A, Deutsch, A. Investigation of the migration/proliferation dichotomy and its impact on avascular glioma invasion. Math Model Nat Phenom 2012;7:105–35. https://doi.org/10.1051/mmnp/20127106. Search in Google Scholar

45. Kolobov, A, Gubenov, V, Polezhaev, A. Autowaves in the model of infiltrativetumor growth with migration-proliferation dichotomy. Math Model Nat Phenom 2011;6:27–38. https://doi.org/10.1051/mmnp/20116703. Search in Google Scholar

46. Saut, O, Lagaert, JB, Colin, T, Fathallah-Shaykh, HM. A multilayer grow-or-go model for GBM: effects of invasive cells and anti-angiogenesis on growth. Bull Math Biol 2014;76:2306–33. https://doi.org/10.1007/s11538-014-0007-y. Search in Google Scholar

47. Martínez-González, A, Calvo, GF, Pérez Romasanta, LA, Pérez-García, VM. Hypoxic cell waves around necrotic cores in glioblastoma: a biomathematical model and its therapeutic implications. Bull Math Biol 2012;74:2875–96. https://doi.org/10.1007/s11538-012-9786-1. Search in Google Scholar

48. Kim, Y, Roh, S, Lawler, S, Friedman, A. miR451 and AMPK mutual antagonism in glioma cell migration and proliferation: a mathematical model. PLoS ONE 2011;6:e28293. https://doi.org/10.1371/journal.pone.0028293. Search in Google Scholar

49. Giese, A, Kluwe, L, Laube, B, Meissner, H, Berens, ME, Westphal, M. Migration of human glioma cells on myelin. Neurosurgery 1996;38:755–64. https://doi.org/10.1227/00006123-199604000-00026. Search in Google Scholar

50. Tektonidis, M, Hatzikirou, H, Chauvière, A, Simon, M, Schaller, K, Deutsch, A. 2011 Identification of intrinsic in vitro cellular mechanisms for glioma invasion. J Theor Biol 2011;287:131–47. https://doi.org/10.1016/j.jtbi.2011.07.012. Search in Google Scholar

51. Pham, K, Chauviere, A, Hatzikirou, H, Li, X, Byrne, HM, Cristini, V, et al.. Density-dependent quiescence in glioma invasion: instability in a simple reaction–diffusion model for the migration/proliferation dichotomy. J Biol Dynam. 2012;6:54–71. https://doi.org/10.1080/17513758.2011.590610. Search in Google Scholar

52. Lee, HG, Kim, Y. The role of the microenvironment in regulation of CSPG-driven invasive and non-invasive tumor growth in glioblastoma. Jpn J Ind Appl Math 2015;32:771–805. https://doi.org/10.1007/s13160-015-0188-2. Search in Google Scholar

53. Kim, Y, Roh, S. A hybrid model for cell proliferation and migration in glioblastoma. Discrete Continuous Dyn Syst Ser B 2013;18:969–1015. https://doi.org/10.3934/dcdsb.2013.18.969. Search in Google Scholar

54. Gerlee, P, Nelander, S. The impact of phenotypic switching on glioblastoma growth and invasion. PLoS Comput Biol 2012;8:e1002556. https://doi.org/10.1371/journal.pcbi.1002556. Search in Google Scholar

55. Frieboes, HB, Zheng, X, Sun, CH, Tromberg, B, Gatenby, R, Cristini, V. An integrated computational/experimental model of tumor invasion. Cancer Res 2006;66:1597–604. https://doi.org/10.1158/0008-5472.CAN-05-3166. Search in Google Scholar

56. Ilina, O, Friedl, P. Mechanisms of collective cell migration at a glance. J Cell Sci 2009;122:3203–8. https://doi.org/10.1242/jcs.036525. Search in Google Scholar

57. Khain, E, Katakowski, M, Hopkins, S, Szalad, A, Zheng, X, Jiang, F. Collective behavior of brain tumor cells: the role of hypoxia. Phys Rev E Stat Nonlinear Soft Matter Phys 2011;83:031920. https://doi.org/10.1103/PhysRevE.83.031920. Search in Google Scholar

58. Friedl, P, Noble, PB, Walton, PA, Laird, DW, Chauvin, PJ, Tabah, RJ, et al.. Migration of coordinated cell clusters in mesenchymal and epithelial cancer explants in vitro. Cancer Res 1995;55:4557–60. Search in Google Scholar

59. Aubert, M, Badoual, M, Christov, C, Grammaticos, B. A model for glioma cell migration on collagen and astrocytes. J R Soc Interface 2008;5:75–83. https://doi.org/10.1098/rsif.2007.1070. Search in Google Scholar

60. Kim, Y, Lawler, S, Nowicki, MO, Chiocca, EA, Friedman, A. A mathematical model for pattern formation of glioma cells outside the tumor spheroid core. J Theor Biol 2009;260:359–71. https://doi.org/10.1016/j.jtbi.2009.06.025. Search in Google Scholar

61. Turner, S, Sherratt, JA. Intercellular adhesion and cancer invasion: a discrete simulation using the extended Potts model. J Theor Biol 2002;216:85–100. https://doi.org/10.1006/jtbi.2001.2522. Search in Google Scholar

62. Axpe, E, Lopez-Euba, T, Castellanos-Rubio, A, Merida, D, Garcia, JA, Plaza-Izurieta, L. Detection of atomic scale changes in the free volume void size of three-dimensional colorectal cancer cell culture using positron annihilation lifetime spectroscopy. PloS One 2014;9:e83838. https://doi.org/10.1371/journal.pone.0083838, https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3879280/. Search in Google Scholar

63. Rodriguez, ML, McGarry, PJ, Nathan, J. Sniadecki. Review on cell mechanics: experimental and modeling approaches. Appl Mech Rev 2013;65:1–41. https://doi.org/10.1115/1.4025355. Search in Google Scholar

64. Schaller, G, Meyer-Hermann, M. Multicellular tumor spheroid in an off-lattice Voronoi–Delaunay cell model. Phys Rev E Stat Nonlin Soft Matter Phys 2005;71:051910. https://doi.org/10.1103/PhysRevE.71.051910. Search in Google Scholar

65. Hartono, D, Liu, Y, Tan, PL, Then, XY, Yung, LY, Lim, KM. On-chip measurements of cell compressibility via acoustic radiation. Lab Chip 2011;11:4072–80. https://doi.org/10.1039/c1lc20687g. Search in Google Scholar

66. Guz, N, Dokukin, M, Kalaparthi, V, Sokolov, I. If cell mechanics can be described by elastic modulus: study of different models and probes used in indentation experiments. Biophys J 2014;107:564–75. https://doi.org/10.1016/j.bpj.2014.06.033. Search in Google Scholar

67. Li, CK. The glucose distribution in 9L rat brain multicell tumour spheroids and its effect on cell necrosis. Cancer 1982;50:2066–73. Search in Google Scholar

68. Stein, A, Demuth, T, Mobley, D, Berens, M, Sander, L. A mathematical model of glioblastoma tumor spheroid invasion in a three-dimensional in vitro experiment. Biophys J. 2007;92:356–65; https://doi.org/10.1529/biophysj.106.093468. Search in Google Scholar

69. Colombo, MC, Giverso, C, Faggiano, E, Boffano, C, Acerbi, F, Ciarletta, P. Towards the Personalized treatment of glioblastoma: integrating patient-specific clinical data in a continuous mechanical model. PloS One. 2015;10:e0143032. Available from: https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0132887; https://doi.org/10.1371/journal.pone.0132887. Search in Google Scholar

70. Andolfi, L, Bourkoula, E, Migliorini, E, Palma, A, Pucer, A, Skrap, M, et al.. Investigation of adhesion and mechanical properties of human glioma cells by single cell force spectroscopy and atomic force microscopy. PLoS ONE 2014;9:e112582. https://doi.org/10.1371/journal.pone.0112582. Search in Google Scholar

71. Conway, JH, Sloane, NJA. Sphere packings, lattices and groups, 3rd ed., vol. 21. New York: Springer-Verlag; 1999. Search in Google Scholar

72. Byrne, HM. Mathematical biomedicine and modeling avascular tumor growth. In: Mathematics and life sciences. The Mathematical Institute, University of Oxford. Available from: https://core.ac.uk/download/pdf/9695564.pdf. Search in Google Scholar

73. Roose, T, Chapman, SJ, Maini, PK. Mathematical models of avascular tumor growth. SIAM Rev 2007;49:179–208. https://doi.org/10.1137/S0036144504446291. Search in Google Scholar

74. O’Brien, FJ, Harley, BA, Waller, MA, Yannas, IV, Gibson, LJ, Prendergast, PJ. The effect of pore size on permeability and cell attachment in collagen scaffolds for tissue engineering. Technol Health Care 2007;15:3–17. Available from: http://iospress.metapress.com/content/103189/. Search in Google Scholar

75. German, Carrie L, Madihally, Sundararajan V. Applications of computational modelling and simulation of porous medium in tissue engineering. Computation 2016;4:7. Available from: https://doi.org/10.3390/computation4010007: https://www.mdpi.com/2079-3197/4/1/7. Search in Google Scholar

76. McCarty, WJ, Johnson, M. The hydraulic conductivity of Matrigel. Biorheology 2007;44:303–17. Available from: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2562030/. Search in Google Scholar

77. Murphy, CM, Haugh, MG, O’Brien, FJ. The effect of mean pore size on cell attachment, proliferation and migration in collagen-glycosaminoglycan scaffolds for bone tissue engineering. Biomaterials 2010;31:461–6. https://doi.org/10.1016/j.biomaterials.2009.09.063. Search in Google Scholar

78. Siddiqui, MU, Arif, AF, Bashmal, S. Permeability-selectivity analysis of microfiltration and ultrafiltration membranes: effect of pore size and shape distribution and membrane stretching. Membranes 2016;6:40. https://doi.org/10.3390/membranes6030040. Search in Google Scholar

79. Kohles, SS, Wilson, CG, Bonassar, LJ. A mechanical composite spheres analysis of engineered cartilage dynamics. J Biomech Eng 2007;129:473–80. https://doi.org/10.1115/1.2746366. Search in Google Scholar

80. Landau, L, Lifshitz, E. Theory of elasticity, course of theoretical physics. Bristol: Pergamon Press; 1970. Search in Google Scholar

81. Bower, AF. Applied mechanics of solids, vol. 10. Boca Raton: CRC Press; 2009: p. 1439802475. Search in Google Scholar

82. Boris, JP, Book, DL. Solution of continuity equation by the method of flux-corrected transport. In: Killeen, J, editor. Methods in computational physics. Advices in research and applications, v.16 Controlled Fusion. New York: Academic Press; 1976: pp. 85–129. Search in Google Scholar

83. Rianna, C, Radmacher, M. Comparison of viscoelastic properties of cancer and normal thyroid cells on different stiffness substrates. Eur Biophys J 2017;46:309–24. https://doi.org/10.1007/s00249-016-1168-4. Search in Google Scholar

84. Ramanujan, S, Pluen, A, McKee, TD, Brown, EB, Boucher, Y, Jain, RK. Diffusion and convection in collagen gels: implications for transport in the tumor interstitium. Biophys J 2002;83:1650–60. https://doi.org/10.1016/S0006-3495(02)73933-7. Search in Google Scholar

85. Moreno-Arotzena, O, Meier, JG, Del Amo, C, García-Aznar, JM. Characterization of fibrin and collagen gels for engineering wound healing models. Materials 2015;8:1636–51. Search in Google Scholar

86. Serpooshan, V, Quinn, TM, Muja, N, Nazhat, SN. Hydraulic permeability of multilayered collagen gel scaffolds under plastic compression-induced unidirectional fluid flow. Acta Biomater 2013;9:4673–80. https://doi.org/10.1016/j.actbio.2012.08.031. Search in Google Scholar

87. Qazi, H, Shi, Z-D, Tarbell, JM. Fluid shear stress regulates the invasive potential of glioma cells via modulation of migratory activity and matrix metalloproteinase expression. PLoS ONE 2011;6:e20348. https://doi.org/10.1371/journal.pone.0020348. Search in Google Scholar

88. Pogoda, K, Chin, L, Georges, PC, Byfield, FJ, Bucki, R, Kim, R. Compression stiffening of brain and its effect on mechanosensing by glioma cells. New J Phys 2014;16:075002. https://doi.org/10.1088/1367-2630/16/7/075002. Search in Google Scholar

89. Raub, CB, Putnam, AJ, Tromberg, BJ, George, SC. Predicting bulk mechanical properties of cellularized collagen gels using multiphoton microscopy. Acta Biomater 2010;6:4657–65. Available from: https://doi.org/10.1016/j.actbio.2010.07.004: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3373188/. Search in Google Scholar

90. Trickey, WR, Baaijens, FP, Laursen, TA, Alexopoulos, LG, Guilak, F. Determination of the Poisson’s ratio of the cell: recovery properties of chondrocytes after release from complete micropipette aspiration. J Biomech 2006;39:78–87. Search in Google Scholar

91. Lu, Y, Franze, K, Seifert, G, Steinhäuser, C, Kirchhoff, F, Wolburg, H, et al.. Viscoelastic properties of individual glial cells and neurons in the CNS. Proc Natl Acad Sci USA 2006;103:17759–64. https://doi.org/10.1073/pnas.0606150103. Search in Google Scholar

92. Hegedus, B, Zach, J, Czirok, A, Lovey, J, Vicsek, T. Irradiation and taxol treatment result in non-monotonous, dose-dependent changes in the motility of glioblastoma cells. J Neuro Oncol 2004;67:147–57. Search in Google Scholar

93. Demuth, T, Hopf, N, Kempski, O, Sauner, D, Herr, M, Giese, A. Migratory activity of human glioma cell lines in vitro assessed by continuous single cell observation. Clin Exp Metastasis 2000;18:589–97. Search in Google Scholar

94. Mueller-Klieser, WF, Sutherland, RM. Oxygen consumption and oxygen diffusion properties of multicellular spheroids from two different cell lines. Adv Exp Med Biol 1984;180:311–21. https://doi.org/10.1007/978-1-4684-4895-5_30. Search in Google Scholar

95. Griguer, CE, Oliva, CR, Gobin, E, Marcorelles, P, Benos, DJ, Lancaster, JR, et al.. CD133 is a marker of bioenergetic stress in human glioma. PLoS ONE 2008;3:e3655. https://doi.org/10.1371/journal.pone.0003655. Search in Google Scholar

96. Ke, LD, Shi, YX, Im, SA, Chen, X, Yung, WK. The relevance of cell proliferation, vascular endothelial growth factor, and basic fibroblast growth factor production to angiogenesis and tumorigenicity in human glioma cell lines. Clin Canc Res 2000;6:2562–72. Search in Google Scholar

97. Vaupel, P. The role of hypoxia-induced factors in tumor progression. Oncol 2004;9:10–17. https://doi.org/10.1634/theoncologist.9-90005-10. Search in Google Scholar

98. Drasdo, D, Höhme, S. A single-cell-based model of tumor growth in vitro: monolayers and spheroids. Phys Biol 2005;2:133–47. https://doi.org/10.1088/1478-3975/2/3/001. Search in Google Scholar

99. Kim, Y. Regulation of cell proliferation and migration in glioblastoma: new therapeutic approach. Front Oncol 2013. https://doi.org/10.3389/fonc.2013.00053 [Accessed 12 Aug 2020]. Search in Google Scholar

100. Vader, D, Kabla, A, Weitz, D, Mahadevan, L. Strain-Induced alignment in collagen gels. PLoS ONE 2009;4:e5902. https://doi.org/10.1371/journal.pone.0005902. Search in Google Scholar

101. Rejniak, KA. A single-cell approach in modeling the dynamics of tumor microregions. Math Biosci Eng 2005;2:643–55. https://doi.org/10.3934/mbe.2005.2.643. Search in Google Scholar

102. Guiot, C, Delsanto, PP, Deisboeck, TS. Morphological instability and cancer invasion: a ‘splashing water drop’ analogy. Theor Biol Med Model 2007;4. https://doi.org/10.1186/1742-4682-4-4. Search in Google Scholar

Supplementary Material

The online version of this article offers supplementary material (https://doi.org/10.1515/jib-2020-0027).