Elucidating the Role of Matrix Porosity and Rigidity in Glioblastoma Type IV Progression

: The highly infiltrating nature of glioma cells is the major cause for the poor prognosis of brain malignancies. Motility, proliferation, and gene expression of cells in natural and synthetic gels have been analyzed by several authors, yet quantitative studies elucidating the role of matrix porosity and rigidity in the development of whole malignant masses are missing. Here, an experimental-computational framework is introduced to analyze the behavior of U87-MG cells and spheroids in compact hyaluronic acid gels (HA), replicating the brain parenchyma; and ﬁbrous collagen gels (COL), resembling the organized structures of the brain. Experimentally it was observed that individual U87-MG cells in COL assumed an elongated morphology within a few hours post inclusion (p.i.) and travelled longer distances than in HA. As spheroids, U87-MG cells rapidly dispersed into COL resulting in inﬁltrating regions as large as tumor cores ( ≈ 600 µ m, at 8 days p.i.). Conversely, cells in HA originated smaller and denser inﬁltrating regions ( ≈ 300 µ m, at 8 days p.i.). Notably, COL tumor core size was only 20% larger than in HA, at longer time points. Computationally, by introducing for the ﬁrst time the e ﬀ ects of matrix heterogeneity in our numerical simulations, the results conﬁrmed that matrix porosity and its spatial organization are key factors in priming the inﬁltrating potential of these malignant cells. The experimental-numerical synergy can be used to predict the behavior of neoplastic masses under diverse conditions and the e ﬃ cacy of combination therapies simultaneously aiming at killing cancer cells and modulating the tumor microenvironment.


Introduction
Glioblastoma type IV (GBM) is the malignancy of the central nervous system with the poorest prognosis [1]. Standard treatments for glioblastoma type IV currently consist of surgical resection followed by radiotherapy and chemotherapy [2]. The median survival for GBM patients is nearly 14 months and the 5-year survival rate is only of 5%. A distinctive feature of GBM cells is their ability to infiltrate the host tissue forming multiple micro-metastases which cannot be surgically removed or efficiently killed by molecular anticancer therapies. Spreading of GBM cells is believed to occur through two major pathways: across the extracellular matrix (ECM), which is mainly composed of a dense based on purely physical reasoning we can explain and model the experimentally observed behavior of U-87MG cells and aggregates.
To this aim 3D scaffolds are realized using hyaluronic acid hydrogels (HA), mimicking the dense and gel-like brain parenchyma, and collagen hydrogels (COL), resembling, at first approximation, the ordered structures of the brain. Then, morphology and migration patterns of individual malignant cells of human glioblastoma type IV cells (U87-MG) are monitored via fluorescent and electron microscopies. The complex growth behavior of spheroids encapsulated in HA and COL are also documented over time, discriminating between the dense tumor core and the surrounding infiltrating region. Finally, a multiphase computational model is employed to predict tumor growth and quantitatively elucidate the role of matrix morphology and porosity on the macroscopic behavior of GBM cells.

Cell Culture and Spheroid Formation
Human glioblastoma U87-MG cells (American Type Culture Collection, ATCC) were selected for their high invasiveness and used as a preclinical model of glioblastoma type IV. Cells were cultured in standard Essential Eagle's Minimum Essential Medium (EMEM) supplemented with 10% fetal bovine serum (FBS) and 1% antibiotics (Euroclone). Medium was changed every 2 days and cells in all experiments were used at passage 8, or lower. U87-MG GFP + were kindly donated by Dr. Davide De Petri Tonelli at the Italian Institute of technology in Genova. U87-MG GFP + spheroids were generated using nonadhesive GravityTRAP TM 96 multiwell plates according to the manufacturer's protocol. In brief, cells were washed twice with PBS, detached by trypsin and counted. Then, 500 cells in 60 µL were dropped inside each well of the 96 multiwell plates. Multiwell plates were centrifuged at 250 RCF (relative centrifugal force) for 2 min at room temperature. Spheroids were kept in culture for 24 h prior encapsulation into hydrogels.

Preparation of Hydrogels and Cell Encapsulation
Collagen type I solution (3 mg/mL) was mixed with a ratio of 4:1 (v/v) with a phosphate buffered saline solution with EMEM (5×) with phenol red, 0.1 M NaOH, and DNA/RNA free sterile water adapted to cell culture to yield a collagen concentration of 2 mg/mL at pH of 7.4. Corgel™ hyaluronic acid hydrogel was prepared according to supplier instructions (Lifecore Biomedical LLC, Chaska, MN, USA) with minor adjustments [28]. Briefly, polymer was dissolved in horseradish peroxidase in PBS (10 U/mL) solution provided by the supplier to obtain a solution of 5 mg/mL. Afterwards, 50 µL of polymer solution was mixed with cells resuspended in 60 µL of EMEM completed medium. Crosslinking reaction was initiated by addition of hydrogen peroxide (H 2 O 2 ) with a final ratio of 1:25 (v/v). In both hydrogel preparations, U87-MG GFP + were encapsulated at a final concentration of 2 × 10 6 cells/mL, as a single cell suspension. Volume calculations were performed to achieve the above-mentioned cellular concentration on both gels. In the case of spheroids, one individual spheroid was encapsulated within each single hydrogel. To this aim, spheroids were added individually to the different hydrogels (COL and HA) in a volume of 60 µL each. Hydrogels were kept in culture and medium was changed every 2 days.

Mechanical Characterization of Hydrogels
Collagen type I and hyaluronic acid hydrogels were characterized mechanically via confined compression tests on a dual column universal testing machine (Instron 3365). Both hydrogels were deposited on 12-well culture plates (diameter 22 mm) and tested using a 20 mm diameter stainless steel piston with a loading rate of 1 mm/min. This allowed the liquid to permeate through the gap between the well and piston. The aggregate modulus H a was extracted from the initial slope of the stress-strain curves, corresponding to a deformation smaller than 5%. Then the aggregate modulus was used to calculate the Young's modulus E according to the equation , where the Poisson ratio was assumed to be ν = 0.45.

Scanning Electron Microscopy (SEM)
At different time points, hydrogels were collected for SEM analysis. In brief, samples were carefully washed twice with PBS and afterwards fixed in a solution of 2% Glutaraldehyde and 2% Paraformaldehyde in Na-Cacodylate buffer 0.1 M for 2 h at room temperature. Then, hydrogels were immersed in Na-Cacodylate 0.1 M buffer solution with 1% OsO 4 and 1.5% Hexayanoferrate for 2 h, at room temperature. Samples were then washed in distillated water and dehydrated by serial immersions in ethanol solutions (30%, 50%, 70%, 90%, 96%, and 100%), for 10 min each. Then, hydrogels were immersed in sequential series of ethanol/bis (trimethylsilyl) amine (HMDS) mixtures (3:1; 1:1; 1:3) for 20 min each. Samples were kept overnight in pure HMDS solution until total evaporation of the solvent. Afterwards, samples were sputtered with 10 nm of gold and images were acquired with JEOL JSM 6490-LA (JEOL, Tokyo, Japan) microscope with accelerating voltage of 10 kV.

Optical Microscopy Analyses
Spatial distribution and morphology of encapsulated U87-MG GFP + cells within collagen type I and hyaluronic acid hydrogels was accessed by the use of a fluorescent inverted microscope (Leica DMI6000 B, Milano, Italy) equipped with a high-speed camera (Leica DFC360 FX, Milano, Italy). Images were acquired by using bright field and 488 nm filters. Spheroid growth, inside the hydrogels, was followed by the use of the same optical microscope. Images were acquired over time using the above cited filters. Tumor size was calculated by measuring the spheroid perimeter via ImageJ64 tools (https://imagej.nih.gov/ij/).

Time Lapse Acquisition and Cell Trajectory Analyses
Cell behavior within collagen type I and hyaluronic acid hydrogels was monitored using an automated all-in-one microscope chamber (Eclipse TiE, Nikon, and Okolab), with controlled temperature and humidity. Prior to time-lapse acquisition, U87-MG GFP + cells were encapsulated within both different hydrogels, as described previously. Hydrogels were prepared inside 48 multiwell-plates precoated with agarose to avoid cell adhesion at the bottom of the wells. Cells were encapsulated at a density of 5 × 10 5 cells/mL in both hydrogels (lower than the initial cell concentration) to allow better cell tracing. Hydrogels were cultured in EMEM supplemented with 10% FBS. Phase-contrast and fluorescent images were captured every 10 min for 24 h. Moving trajectories of 50 cells on each hydrogel (five different hydrogels for studied condition) were analyzed.

MTT Assay
The proliferation of encapsulated U87-MG cells within both gels was evaluated at different time points, using thiazolyl blue staining (MTT, Sigma-Aldrich, St. Louis, MO, USA), in which absorbance data correlate with cell metabolic activity. In brief, at each time point, the medium was removed and replaced with 1 mL fresh medium supplemented with MTT stock solution at a final concentration of 5 mg/mL. After 3 h incubation, the medium was removed, and the converted dye was solubilized by the addition of 1 mL ethanol absolute. To obtain a homogenous solution, hydrogels were disrupted mechanically via a potter [29]. Absorbance was measured at 570 nm with background subtraction at 670 nm.

Cellular Metrics: Single Cell Mobility and Spheroid Progression Kinetics
NIS software from Nikon was used to quantify single cell mobility inside both hydrogels. Tumor spheroid metastatic progression inside collagen type I and hyaluronic acid hydrogels was analyzed via ImageJ64 tools (https://imagej.nih.gov/ij/). At each time point, the region occupied by the tumor was separated into two domains-the Tumor Core and the Tumor Infiltrating Region. The size of both domains was measured over time.

Statistics
Values of U87-MG proliferation inside collagen type I and hyaluronic acid hydrogels obtained over time are expressed as a mean ± standard deviation, with n = 3 for each group of culture. Spheroid metastasis kinetic values are expressed as a mean ± standard deviation with n = 3, for each group of culture from three independent assays. Statistical significance was determined using Student's t-test multiple comparison procedure at a confidence level of 95% ** (p < 0.005), * (p < 0.05), via GraphPad (www.Graphpad.com).
The coefficient of determination R 2 is used to quantify the accuracy of the computational model against the experimental data. It is defined as where y i are the experimental data, y is the average value of experimental data,ŷ i are the results predicted with the model. A R 2 of 1 indicates perfect agreement between computational predictions and experimental data.

Computational Model
Our tumor growth model [30,31] is a continuum model of the multiphysics type and includes aspects of Thermodynamically Constrained Averaging Theory [32] where the model derivation proceeds systematically from known microscale relations to mathematically and physically consistent larger scale relations. Other continuum growth models are of the multiparameter type [33][34][35] based on mixture theory where the relevant balance equations are written directly at the level of interest and the thermodynamic consistency is satisfied at the same level. Our model comprises four phases: the solid phase, i.e., the ECM and three fluid phases: tumor cells, TCs, host cells, HCs, interstitial fluid, IF. In this paper, to simulate the growth of U87-MG cells cultured in collagen type I hydrogels and in HA hydrogels, the system is reduced to three phases, ECM, TCs, IF, by taking the volume fraction of HCs, ε h , equal to zero.
ECM is represented by an elasto-visco-plastic porous solid, with porosity ε, filled by TCs and IF. It is recalled that pore size as used in Yang et al. (2010) is not a direct measure of porosity. Both can be measured in biological tissues via scanning electron microscopy (SEM) and microcomputed tomography and from these a relation between the two can be established. In Karageorgiou and Kaplan [36], there are some values of both pore size and porosity for biological tissues such as collagen gels and hyaluronic acid gels. Our chosen values of the porosity for modeling purposes are consistent with those values.
At each point of the continuum, the representative elementary volume (REV) of the multiphase computational model comprises a certain portion of ECM, TCs, and IF, as schematically shown in Figure 1. The TCs are divided in living tumor cells, LTCs, and necrotic tumor cells, NTCs. In the IF phase, transported species (nutrients) are considered. The volume fraction of the solid phase ECM is ε s = 1 − ε. The other phases, tumor cells (ε t ) and interstitial fluid (ε l ), occupy the rest of the volume so that the volume fractions for all phases add up to unity.

Governing Equations
The governing equations of the model used here are the mass balance equation for ECM, IF, TCs, nutrients, and the linear momentum balance equation for the mixture ECM+IF+TCs. The following equations are obtained from the general forms by introducing some simplifications and closure relationships (e.g., a Fickian type equation for diffusion of species, a generalized Darcy's equation for flow of the fluid phases, etc.). The primary variables of the model are differential pressures between phases t, h, l, p hl and p th , IF pressure p l , nutrient mass fraction, and displacement u s of solid phase (ECM). It is recalled that p th + p hl = p tl .
The saturation degree of a fluid phase α is S α = ε α /ε, with α = t, h, l for TCs, HCs, and IF, respectively. In this paper the saturation degree of HCs is S h = 0. From Equation (1) the constraint condition for saturation degrees for the system with just three phases follows The mass balance equation of the ECM is The mass balance equation of TCs reads The mass balance equation of IF reads where µ α is the dynamic viscosity, k α rel is relative permeability, ρ α is the density. Additionally, d s is the Eulerian rate of strain tensor, k the intrinsic permeability tensor of the ECM, and where D nl e f f is the effective diffusivity of the nutrient species in the extracellular space, ω nl the mass fraction of nutrient species n, and nl→t M is the mass of nutrient consumed by tumor cells via metabolism and growth. In our simulations the oxygen is the sole nutrient considered.
A further governing equation of the model is the linear momentum balance of the mixture expressed in rate form as where the interaction between the solid and the three fluid phases is accounted for through the effective stress t s e f f in the sense of porous media mechanics t s e f f = t s + αp s 1 (8) where 1 is the unit tensor, t s is the total stress tensor in the solid phase, α is the Biot's coefficient α = 1 − K/K s , with K the compressibility of the empty ECM. In the modeled problem, K/K s tends to zero hence we can assume a Biot's coefficient equal to 1. Being S h = 0 the solid pressure p s is given as, For the fluid phases the momentum balance equations can be simplified to where R α is the resistance tensor that accounts for the frictional interactions between phases when there is flow and different phases have different velocities. With the generalized Darcy law, the resistance tensor is with µ α the dynamic viscosity that is defined to model the adhesion between tumor cells and ECM.
The constitutive relationships for the fluid phases and the solid phase as well as the constitutive equations for growth and nutrient consumption have been extensively dealt with in [30,31,37]. Diffusion of the chemical species dissolved in the interstitial fluid follows Fick's law [30,31]. The constitutive relationship for the deformable solid phase (ECM) follows an elasto-viscoplastic law of Perzyna type in large deformation regime. Until a yield limit of the effective stress is reached, the solid behavior is contained in the elastic range and then it becomes visco-plastic with a constitutive flow rule obtained from rate-independent hardening plasticity by the introduction of the viscosity. For the numerical implementation of the equations see Appendix A.

Description of Geometry, Boundary, and Initial Conditions of the Simulations
The geometry of the problems at hand is simulated with a sphere segment in axisymmetric conditions with radius of 2500 µm. At the initial time instant, the multicellular tumor spheroid (MTS) is composed of three phases: (i) the TCs phase, in the red area of Figure 2 with radius of 200 µm for the simulation of tumor growing in collagen and 250 µm for the simulation of tumor growing in hyaluronic acid; (ii) the IF phase which fills the whole domain; (iii) the ECM solid phase which also fills the whole domain. The atmospheric pressure is taken as the reference pressure and the initial IF pressure is zero Pa in the entire domain. The initial volume fraction of TCs is 0.02, the HCs are not present. Boundary conditions are imposed as indicated in Figure 2. To allow IF flux at the outer boundary the IF pressure is fixed there to zero Pa. Due to the symmetry of the problem there is no flux normal to the radius of the sphere segment. Oxygen is the sole nutrient species and its mass fraction is fixed to ω nl env = 4.2 · 10 −6 at the outer boundary and throughout the computational domain at initial time.
All model parameters are listed in Tables 1-3 classified by type. Coefficients γ nl growth and γ nl 0 control the TCs uptake of oxygen; they allow modeling consumption due to TCs growth and their metabolism.

Parameter Symbol Value Unit
Diffusion coefficient of oxygen in interstitial fluid D nl As said above the ECM is treated as a deformable solid phase with a visco-elastic constitutive law. A different stiffness of the solid skeleton corresponds to a different volume fraction of the solid phase. In the cases with a nonuniform porosity, the following relations are assumed for the Young's modulus, the intrinsic permeability k and the parameter a, which is a constant parameter related to ECM porosity, used in the saturation-capillary pressure relationship A and B refer to the cases with different porosity dealt with in Section 4.6. The Young's modulus E = 500 Pa, the intrinsic permeability k = 1.8 × 10 −15 m 2 , and the coefficient a = 590 Pa of zone B, with porosity ε = 0.5, are considered as the reference values.

Morphological and Mechanical Characterization of Hydrogels for the Inclusion of U87-MG Tumor Cells
Collagen is the most ubiquitous protein in the ECM and forms random networks of fibers depending on the assembling conditions [38,39]. In the brain, collagen is mostly found along the vascular network together with other proteins, such as laminins and fibronectin, and provides mechanical support as well as important cues for tumor cell dynamics [3,14,17] Another major component of the brain ECM is hyaluronic acid. This is a naturally occurring glycosaminoglycan that has the capability to bind to cell specific surfaces receptors (CD44), which are necessary for cell anchorage, homeostasis, and proliferation [15,[40][41][42]. For this, both COL and HA have been used in different forms to reproduce different portions of the brain tissue.
Cells of glioblastoma type IV transfected with green florescence protein-U87-MG GFP + -were included into hydrogels of COL and HA either individually or as preformed spheroids (Figure 3). For individual cells, these were first expanded on a dish, then resuspended and added to the natural polymeric hydrogel. For spheroids, cell clusters of about 200 µm were preformed and eventually included in the polymeric hydrogel. Finally, the COL and HA mixtures were deposited at the bottom of polystyrene multiwell plates, precoated with agarose, polymerized, and left in cell culture media for up to 8 days.
Knowing that the microstructural and biomechanical features of the extracellular matrix may severely affect the growth and spreading of malignant tissues, the morphology and deformability of the hydrogel scaffolds were characterized. First, scanning electron microscopy (SEM) was used to assess the microstructure of COL and HA. Figure 4A and its lateral insets show the organization of collagen fibers within the hydrogel. COL fibers appear entangled, with a diameter of approximately 10 nm, and randomly arranged to form a matrix with pores varying in size between about 20 and 70 nm, as from SEM images. Differently, HA presents a more compact organization of thinner (about 5 nm) and curly fibers ( Figure 4B). Importantly, the magnified SEM image for HA does not reveal any obvious, regular porosity of the scaffold. This supports the notion by which HA hydrogels can be used for mimicking the soft and dense brain parenchyma, whereas fibers within COL tend to resemble morphological features of ordered structures in the brain [3]. Subsequently, the mechanical deformability of the hydrogels was assessed under wet conditions by standard compression testing. As expected, COL presented a higher Young's modulus E as compared to HA. Specifically, as documented in the chart of Figure 4C, E was 5.16 ± 1.40 kPa for COL and 0.20 ± 0.03 kPa for HA, assuming a Poisson coefficient ν = 0.45. These values are in good agreement with previously reported characterizations for human brain tissue [43][44][45]. COL presents a fibrous structure with a variable porosity. HA presents a complex dense structure. (C) Young's modulus of COL and HA hydrogels. Scale bar = 5 µm and 500 nm for the low and high magnification images, respectively. Statistical analysis was determined using t-test multiple comparison with a confidence level of 95% (p < 0.05). GraphPad was used.

Growth of Individual U87-MG Cells in Different Hydrogels
In studying the proliferation and migration ability of cells within man-made scaffolds, it is fundamental to confirm cell growth over time. Fluorescent and electron scanning microscopy together with MTT assays were performed to document GBM cell proliferation and morphology within COL and HA.
Fluorescent images of U87-MG GFP + cell dispersions in COL ( Figure 5A) and in HA ( Figure 5B) were acquired from day 1 to day 7, post seeding. In both cases, the fluorescent signal significantly increased over time and spreads quite uniformly over the entire hydrogel, thus documenting cell proliferation. Fluorescent images for all time points, namely days 1, 3, 5, and 7, are shown in Figure 6. Cell proliferation via MTT analysis is shown in Figure 7.   In particular, proliferation was statistically significant for the first three time points, namely from day 1 to 3 ** (p = 0.0001 for COL; p = 0.002 for HA); from day 3 to 5 * (p = 0.0126 for COL; p = 0.035 for HA), whereas no significant difference was observed between the last two time points, corresponding to days 5 and 7. Statistical analysis was performed using t-test multiple comparison with a confidence level of 95% (p < 0.05). GraphPad was used.

Cell Imaging
Over the course of the culture, spatial distribution and morphology of encapsulated U87-MG GFP + cells within collagen type I and hyaluronic acid hydrogels was monitored using a fluorescent inverted microscope (Leica DMI6000 B) equipped with a high-speed camera (Leica DFC360 FX). Images were acquired by the use of bright field and 488 nm filters. Images were obtained at 10× magnification.

MTT Assay
The proliferation of encapsulated U87-MG cells within both gels was evaluated at different time points, using thiazolyl blue staining (MTT, Sigma-Aldrich), in which absorbance data correlate to cell metabolic activity. In brief, at each time point, the medium was removed and replaced with 1 mL fresh medium supplemented with MTT stock solution at a final concentration of 5 mg/mL. After 3 h incubation, medium was removed and the converted dye was solubilized by the addition of 1 mL ethanol absolute. To obtain a homogenous solution, hydrogels were disrupted mechanically via a potter [29]. Absorbance was measured at 570 nm with background subtraction at 670 nm.
With time, cells were also observed to change their morphology as they were proliferating and moving about the hydrogels. This is clearly documented by SEM images in Figure 5C,D, respectively, for COL and HA. Initially, U87-MG cells presented a rounded shape on both scaffolds. However, at later time points, cells cultured in COL were mostly elongated presenting a bipolar, spindle-shaped morphology, whereas cells cultured in HA mostly preserved the original rounded shape. Figure 8 shows similar morphologies for U87-MG cells in COL already at day 3. It should also be noted that U87-MG cells in HA have been observed to cluster together forming larger aggregates, differently from cells in COL which are rather isolated. This is already evident at day 3 as evidenced by the right column images representative of higher magnifications of those on the left side ( Figure 8). This behavior is also confirmed by the fluorescent images presented in Figure 5E,F. Again, U87-MG cells cultured in COL tend to assume an elongated morphology, which preludes to higher motility and migration potential. Differently, the same cells in HA hydrogels assume a rounded shape and are more prone to clustering as shown in the high magnification images (right column) of Figure 8. Morphology of U87-MG cells was monitored over time, from day 1 to 7, via scanning electron microscopy. Cells included in COL gels were observed to alter their initial rounded shape into a more spindle-like, elongated morphology. This behavior is already observed at earlier time points. Cells included in HA gels grew into small clusters preserving their original rounded shape. Images on the right column panel of each gel (Collagen type I and Hyaluronic acid) represent higher magnifications of the related left ones. Scale bar = 50 µm.
To further characterize the hydrogel-dependent cell migration behaviors, individual U87-MG GFP + cells were monitored by registering their position and measuring the length of the path traced by cell nuclei over time. Tumor cells in COL hydrogels migrated over longer distances as compared to cells in HA, within the same time period. This is reflected by the polar diagrams of Figure 5G,H and in the Supporting Movies S1 and S2. These different cell behaviors are in agreement with previous reports showing that HA favors the formation of cell clusters [14], whereas COL supports more cell spreading [46].

Growth of U87-MG Spheroids in Different Hydrogels
Cellular spheroids of glioblastoma type IV were preformed into GravityTRAP TM ULA plates, via the progressive assembling of individual U87-MG cells. After 1 day, these preformed spheroids, presenting a radius of about 200 µm, were included in the natural polymeric hydrogels and monitored over time for their growth and spatial organization, from 18 h up to 8 days post inclusion. Figure 9A,B shows bright field and fluorescent images of tumor spheroids at different time points for COL and HA, respectively. Five spheroids were used under each condition. Interestingly, these figures show that spheroids included into COL rapidly tend to spread with individual cells rapidly moving about the scaffold already at day 1. This can be readily observed in the bright field and fluorescent images were the high-cell density circle (dark area), which corresponds to the Tumor Core region, is surrounded by a low-cell density corona (bright area), which corresponds to the Tumor Infiltrating region. A similar trend is also observed for the HA case, but it can only be appreciated at later time points, after day 6. The overall radius of the tumor spheroids, including the radius of the core and the thickness of the surrounding infiltrating region, and the radius of the tumor core are charted in Figure 9C,D for the different times points, as measured via fluorescent microscopy, respectively, for COL and HA. For spheroids growing in COL ( Figure 9C), the overall radius varies from about 200 µm at time of encapsulation to 1221.88 ± 22.93 µm at day 8, documenting a 6-fold increase. The high-cell density tumor core radius reaches 551.22 ± 52.56 µm after 8 days of culturing, with a 2.5-fold increase. This data also shows how significant the spreading of U87-MG cells within COL is. As reported in Figure 9C, the overall tumor radius, comprising the surrounding infiltrating region, is roughly two times that of the tumor core for the whole observation period, starting already at 18 h. On the other hand, spheroids growing inside HA were observed to grow from about 200 µm at time of encapsulation to 783.01 ± 145 µm at day 8, documenting only a 3.5-fold increase ( Figure 9D). Additionally, the high-cell density tumor core radius reaches 448.27 ± 69.53 µm after 8 days of culturing, with a 2.5-fold increase. Only between day 3 and day 4, an infiltrating region was observed. This data mostly indicates a delayed and overall lower growth of the infiltrating region for the spheroid cultured in HA as compared to COL. Overall this different behavior in tumor spheroid growth and infiltration appears to match with the different behavior registered for the individual tumor cells (Figure 5), whereby U87-MG in COL mostly take, already after a few hours of culturing, an elongated shape priming towards migration whereas the same cells in HA acquire a rounded shape prone to clustering.

Predicting the Growth of U87-MG Spheroids in Different Scaffolds
The multiphysics computational model described in Section 3, is used to predict the temporal development of GBM spheroids within different microenvironmental conditions. The role of two major mechanical parameters of the hydrogels was systematically assessed, namely the Young's modulus E and the porosity ε, defined as 1-V s /V tot where V s is the volume occupied by the solid porous ECM within the whole domain V tot . Therefore, the larger ε is and the smaller the volume fraction occupied by the solid ECM is, the higher is the porosity of the ECM. For HA, a Young's modulus of 200 Pa was considered (see Figure 2) while the matrix porosity ε was varied between 0.8 and 0.3, and kept homogeneous across the whole domain surrounding the initial spheroid.
1 Figure 9. Behavior of glioblastoma type IV spheroids in hydrogels. (A,B) Microscopy images-bright field (BF), fluorescent (GFP), and merged-of U87-MG GFP + spheroids documenting tumor growth over time COL and HA, respectively. The tumor mass appears as a dense (high fluorescence intensity) inner portion-core-surrounded by a lighter (low fluorescence intensity) outer portion-infiltrating region (scale bar = 50 µm). (C,D) Tumor spheroid progression over time, expressed in terms of overall radius of the tumor mass (grey bars) and radial dimension of the infiltrating region (dark bars) for COL and HA, respectively. Statistical analysis was done using t-test multiple comparison with a confidence level of 95% (p < 0.05). GraphPad was used. Figure 10A shows the variation of the predicted tumor core size as a function of time for different values of ε. For ε = 0.8, a good agreement (R 2 = 0.78) was found between the computational (solid line) and the experimental data for HA (dots). On the other hand, as the porosity ε reduces, a significant decrease in tumor growth rate was observed. At 8 days, the tumor size for ε = 0.5 and 0.3 was, respectively, 18.52% and 39.48% smaller as compared to the ε = 0.8 case (389.62 and 289.35 µm as compared to 478.18 µm, respectively). Indeed, dense matrices (low ε) favor pressure build-up within the porous solid which eventually limits cell proliferation and induces cell death. Then, starting from the conditions resembling the growth of the tumor core in HA, the effect of the Young's modulus E was investigated. Figure 10B reports on the tumor growth over time for E ranging between 100 and 5000 Pa, the latter being the Young's modulus measured for COL (see Figure 4). As expected, an increase in matrix rigidity leads to a decrease in tumor growth. Indeed, following the same argument presented above for ε, a higher mechanical constriction results in lower cell proliferation. Notably, despite the huge variation in Young's modulus, from 100 to 5000 Pa, the corresponding change in tumor size is quite moderate. At 8 days, the tumor size for E = 100 and 5000 Pa was, respectively, 6.07% higher and 6.74% lower as compared to the E = 200 Pa case (507.25 and 445.95 µm as compared to 478.18 µm, respectively). This evidently implies that tumor growth is more sensitive to variations in matrix porosity rather than rigidity. For the E = 200 Pa case there is a good agreement between computational results (solid line) and the experimental data (dots).
As shown in Figure 9, tumors in COL tend to grow faster than in HA. Therefore, in the numerical simulations, matrix porosity was changed to reproduce this behavior. In Figure 10C, the dashed and dash-dot lines reproduce the cases for ε = 0.5 (low porosity) and ε = 0.9 (high porosity). As expected, lower porosities were associated with smaller tumor sizes over time. However, in both cases, the computational predictions were far from the experimental data. Thus, since COL are characterized by fibers along which tumor cells can more easily migrate, it was assumed that the matrix surrounding the initial spheroid could be characterized by alternating radial sectors of low (ε = 0.5 and 0.3) and high (ε = 0.85) porosities. This configuration should favor tumor growth along the low porosity sectors as compared to the denser matrix regions. The dotted and dash-two dots lines in Figure 10C present data for these nonuniform porosity configurations. Within the first 3 days, both cases presented a similar behavior, with almost overlapping tumor growth curves. However, at longer time points, both cases predicted a faster growth as compared to the experimental data. Interestingly, at longer time points, the ε = 0.5/0.85 case showed a trend similar to the experimental data but shifted towards larger sizes. As detailed in Figure 11, the dual porosity configurations (ε = 0.3/0.85 and ε = 0.3/0.5) predicted quite similar tumor growth behaviors, with the latter giving overall smaller tumor volumes because of the overall lower porosity of the matrix. Both simulations with dual configuration of porosity are not able to reproduce the experimental data (dots). More interestingly, the ε = 0.3/0.85 case resulted in tumor core volumes very similar to those estimated in the case of a uniformly porous matrix with ε = 0.9 (high porosity) ( Figure 10C). These numerical results would imply that matrices with nonuniform porosities experience a more rapid tumor growth as compared to uniform and higher porosity matrices. However, none of the two general porosity configurations so far considered (uniform and nonuniform with radial sectors) accurately predict the tumor growth in COL. This is mostly because of the relatively low growth rate observed within the first 4 days of culturing. Therefore, the more realistic configuration for tumors growing in COL resulted in a combination of the above cited two cases: a rim with uniform porosity surrounding the initial tumor spheroid with ε = 0.5, which serves to modulate the initial growth rate; and alternating radial sectors of low (ε = 0.5) and high (ε = 0.85) porosities, characterizing the most exterior portion of the matrix. This case is identified as ε = 0.5 + 0.5/0.85 in the following. For this configuration, a good agreement (R 2 = 0.91) was achieved between the computational (solid line) and the experimental data (dots) for COL, as documented in Figure 10C. Interesting is comparison of the resulting distribution of the tumor volume fractions for the cases ε = 0.3/0.85 and ε = 0.5 + 0.5/0.85 after 8 days of simulated growth, drawn in Figure 12. Since tumor growth influences the porosity, the porosity distribution after 8 days is shown as the second figure for each case. The initial porosity distribution can be seen in Figure 10C. A more developed heterogeneity of the tumor cells distribution appears clearly in Figure 12B in line with what was expected from the porosity distribution of COL shown in Figure 4A. The blue part is the remaining region of the simulated domain not yet occupied by the core. It contains however the corona.

Discussion
In this work, the proliferation and migration propensity of human glioblastoma type IV cells (U87-MG) was studied within different tridimensional microenvironments: hyaluronic acid hydrogels (HA), whose density and rigidity resembled the gel-like properties of the brain parenchyma; collagen hydrogels (COL), whose spatial organization with fibers and pores resembled the ordered structures in the brain, such as the perivascular space and axonal fibers. Although none of these two simple hydrogels can precisely mimic the tridimensional complexity of a human brain tissue, relevant cell features, such as morphology, migration, clustering, and proliferation, were found to closely match in-vivo observations. In particular, glioma cells tend to spread more rapidly in matrices with organized fibers, such as COL, whereas cell clustering is typical of dense and poorly organized matrices, such as HA. This was documented in Figure 5. It was noticed that individual U87-MG cells rapidly assume an elongated shape, with cell membrane ruffling and long filopodial protrusions, on COL. This spindle shaped cell morphology is a clear sign of migratory phenotype and was observed as early as a few hours post-inclusion. Very differently, individual tumor cells on HA preserve their original rounded shape, which is more prone to form large and dense cell clusters. The different cell morphology could indeed be associated with ECM rigidity, as previously documented by Kumar and collaborators [47], which might induce a cell epithelial to mesenchymal transition and facilitate migration [48]. Furthermore, by tracking over time the motion of individual cell nuclei, U87-MG cells were confirmed to have higher motility and travel longer distances on COL as compared to HA. Overall this behavior agrees with in-vivo observations and the notion that glioma cells would more efficiently migrate along organized brain structures (blood vessel walls and axonal fibers) rather than across the dense extracellular space [3,[49][50][51].
This individual cell behavior is then reflected in the overall development of neoplastic masses, as documented in Figure 9. Tumor spheroids included into COL gels grow at faster rates and, more importantly, infiltrate more deeply the surrounding matrix, as compared to HA gels. At 8 days post spheroid inclusion, neoplastic masses in COL have a more developed infiltrating front as compared to HA. On the other hand, tumors in HA present a more dense and smaller infiltrating front again suggesting that cell clustering is preferred to individual cell migration within such matrices. This can be readily appreciated by comparing the individual pictures in Figure 9A,B: the infiltrating front is more diffuse and spread over the entire region of interest in the case of COL ( Figure 9A), whereas it appears more compact and well defined in the case of HA ( Figure 9B). At 8 days, the radius of the HA infiltrating front is still about two times smaller than for COL: the characteristic sizes are 345.55 ± 139.09 µm for HA and 661.65 ± 39.68 µm for COL. Supporting Movies S3 and S4 readily depict all these differences between COL and HA.
It is worth noting the different evolution in time of the tumor core growing in collagen type 1 and hyaluronic acid: in the first case there appears an accelerated trend while in the second case some tendency to a stable state result. In the medium the growth pattern appears linear. This is more clearly shown in Figure 13, where a direct comparison with spheroids growing in cell culture media is presented too. It is evident that spheroids in suspensions do not develop as rapidly and aggressively as spheroids embedded into more realistic environments: at 8 days, tumor spheroids developing in media have roughly half the size of tumors growing in gels (551.22 ± 52.56 µm in COL, 448.27 ± 69.53 µm in HA, and 278.63 ± 17.08 µm in medium). Obviously, spheroids in suspension cannot develop a truly infiltrating front. This data continues to confirm that cell adhesion and interaction with proper extracellular matrices favor cell duplication and spreading. Figure 13. Tumor spheroid growth over time under different conditions. Variation of the tumor core radius (A) and infiltrating region size (B) over time for spheroids suspended in free medium (black bars), in COL hydrogels (gray bars) and in HA hydrogels (white bar). Note that no infiltrating region can be detected in the case of spheroids growing suspended in free medium. Statistical analysis was undertaken using t-test multiple comparison with a confidence level of 95% (p < 0.05). GraphPad was used.
Finally, a multiphysics computational model was demonstrated to accurately predict the temporal growth of the tumor masses for both hydrogel configurations. A soft matrix with uniform porosity showed to give good results for spheroids grown in HA, whereas a nonuniform porosity distribution was required in order to match more closely the experimental data for spheroids cultivated in COL gels. Interestingly, simulations showed that nonhomogeneous porous matrices, modeled by alternating radial sectors with different values of porosity, could more efficiently support the spreading and growth of tumor masses as compared to uniform porosity matrices. This numerical evidence shows that the computational model is able to capture the higher migratory attitude of glioma cells in matrices with a random network of fibers as compared to matrices with a uniform dense structure. Further, the numerical simulations enable one to evaluate the influence of tumor growth on porosity, as illustrated in Figure 12 at day 8 as an exemplifying case. It is here important to note that the parameters for cell proliferation, adhesion, and ECM remodeling were kept identical ** * * * ** * Figure 13. Tumor spheroid growth over time under different conditions. Variation of the tumor core radius (A) and infiltrating region size (B) over time for spheroids suspended in free medium (black bars), in COL hydrogels (gray bars) and in HA hydrogels (white bar). Note that no infiltrating region can be detected in the case of spheroids growing suspended in free medium. Statistical analysis was undertaken using t-test multiple comparison with a confidence level of 95% (p < 0.05). GraphPad was used.
Finally, a multiphysics computational model was demonstrated to accurately predict the temporal growth of the tumor masses for both hydrogel configurations. A soft matrix with uniform porosity showed to give good results for spheroids grown in HA, whereas a nonuniform porosity distribution was required in order to match more closely the experimental data for spheroids cultivated in COL gels. Interestingly, simulations showed that nonhomogeneous porous matrices, modeled by alternating radial sectors with different values of porosity, could more efficiently support the spreading and growth of tumor masses as compared to uniform porosity matrices. This numerical evidence shows that the computational model is able to capture the higher migratory attitude of glioma cells in matrices with a random network of fibers as compared to matrices with a uniform dense structure. Further, the numerical simulations enable one to evaluate the influence of tumor growth on porosity, as illustrated in Figure 12 at day 8 as an exemplifying case. It is here important to note that the parameters for cell proliferation, adhesion, and ECM remodeling were kept identical for both hydrogel conditions. The in silico tests were performed by varying only the elastic modulus E, the value of porosity ε, and its spatial organization. This confirms that major differences in tumor growth between the two gels have to be ascribed to mechanical cues, such as porosity and rigidity.

Conclusions
One of the key issues in glioblastoma progression is the highly infiltrating nature of the tumor mass, which is pathway dependent and marked by the ability of glioma cells to adopt different migratory strategies. In this work, for the first time, a set of quantitative studies was carried out, by means of in vitro and in silico tests, to elucidate the role of matrix porosity, rigidity, and morphology in the development of glioblastoma malignancy. Proliferation and motility of glioblastoma type IV cells were assessed in different extracellular matrices, namely made out of collagen type 1 and hyaluronic acid. Collagen matrices presented organized fibrous structures with a pore size ranging between 20 and 70 nm and a Young's modulus of 5000 Pa. Differently, hyaluronic acid matrices appeared as a compact and disorganized cluster of short fibers with no obvious porosity and a Young's modulus of 200 Pa. Glioblastoma type IV cells were observed to readily proliferate in both matrices, under similar culturing conditions. However, individual cells in collagen matrices exhibited an elongated morphology, similar to fibroblasts, documenting a propensity to migration whereas cells in hyaluronic acid matrices appeared globular. Culturing of tumor spheroids into the two different hydrogel matrices recapitulated the behavior manifested by individual cells. Glioblastoma type IV spheroids were documented to grow more rapidly in collagen matrices forming already after a few hours cell protrusions, possibly aligned with the collagen fibers. On the other hand, the same spheroids mostly developed as compact masses when cultured in hyaluronic acid matrices. Although infiltrating regions around the tumor core were observed for both matrices, this developed earlier and faster for spheroids included in collagen as compared to hyaluronic acid. A tumor growth multiphysics model, accurately predicting the different experimental behaviors in locally homogeneous environments, was able to match also the evolution in time of the tumor spheroids growing in the two different scaffolds by introducing some matrix heterogeneity: the plateau following the initial growth in hyaluronic acid and the faster trend in collagen. The computational results matching the realistic cases as well as those with different porosity morphology performed for investigative aim, confirmed that matrices with nonuniform porosities accelerate the development and tissue infiltration of neoplastic masses. All together, these data demonstrate that mechanical features of extracellular matrices, such as porosity and its spatial organization, are more relevant than matrix rigidity in supporting rapid tumor growth and infiltration. This experimental-computational framework can be used to quantitatively study the influence of porosity in tumor growth and to clarify the role of porosity morphology in the evolution of malignant masses. Furthermore, in silico tests allow to investigate many combinations of parameters and configurations, even in cases which would be impossible to reproduce in the laboratory. By exploiting the synergy of the tests performed in vitro and in silico, it is possible to predict the behavior of the neoplastic masses and their evolution over time, as well as to elucidate which are the most relevant biomechanical cues for combined therapies aiming simultaneously at killing cancer cells and modulating the tumor microenvironment. Finally, these first results about the importance of heterogeneity point to the need of investigating the impact of the stochastic heterogeneity of physical parameters on tumor growth because this heterogeneity is frequently attributed only to genetics [52,53]. This paper has shown that this is not necessarily the case and paves the way to further investigation on this aspect.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2076-3417/10/24/9076/s1. Video S1-Dynamics of individual U87-MG cells in Collagen hydrogels. The movie shows U87-MG cells moving within the COL gel and assuming an elongated shape. In particular, following cell division (lower left quadrant; right side), U87-MG cells rapidly assume an elongated shape and move away from each other. Video S2-Dynamics of individual U87-MG cells in Hyaluronic Acid hydrogels. The movie shows moderate movement of rounded U87-MG cells with the HA gel. In particular, following cell division (upper right quadrant), U87-MG cells stay close to each other forming small clusters and preserving their original rounded shape. Video S3-Temporal evolution of U87-MG spheroids in Collagen hydrogels. The movie shows rapid development of the tumor spheroid with elongation and migration of U87-MG cells away from the original core towards the infiltrating region. An entire portion of the malignant mass tends to move away from the original core. Video S4-Temporal evolution of U87-MG spheroids in Hyaluronic Acid hydrogels. The movie shows a local development of the tumor spheroid with slow rearrangements in shape. No infiltrating region can be detected.  The governing equations are discretized in space and time using a finite element method in space and a finite difference method in time, adopting a quasi-Crank-Nicolson scheme (θWilson method with θ = 0.52) [30,31]. The weak form of the equations is obtained by means of the standard Galerkin procedure. Within each time step the equations are linearized by means of the Newton-Raphson method. The full set of equations is implemented by using Cast3M software, available at: http: //www-cast3m.cea.fr/.
For the solution of the resulting system of equations, Equation (A1), a staggered scheme is adopted with iterations within each time step to preserve the coupled nature of the system.
. The nonlinear coefficient matrices C ij (x), K ij (x) and f i (x) are given in Sciumè et al. [31]. Three computational units are used in the staggered scheme: the first is for the nutrient mass fraction, the second to compute, p hl , p th , and p l , and the third is used to obtain the displacement vector. Within each coupling iteration, the mass balance equation for the nutrient is solved for the mass fraction of the nutrient. Then the other mass balance equations are solved in a fully coupled way for p hl , p th , and p l , yielding here, de facto, p tl and p l . In this second computational unit, at each iteration i the approximate solution p hl , p th , and p l is used to update the mass fraction of the necrotic tumor cells, the mass exchange terms, and the reaction terms. Once convergence is achieved for the second computational unit, the pressure in the cells phases is used to compute the solid pressure. The solid pressure is needed to solve the momentum balance equations. Once convergence is achieved within a time step the procedure can march forward. It is interesting to note that the general procedure allows also for computations with one degree of saturation equal to zero without additional changes. The modular computational structure allows to take into account more than one chemical species, simply adding a computational unit (equivalent to the first one used for the nutrient) for each of the additional chemical species considered. The procedure has been implemented in the code CAST3M (http://www-cast3m.cea.fr) of the French Atomic Energy Commission.
The input parameters for the model are the domain geometry and values of the physical governing parameters, which are listed in Tables 1-3. As output, the model yields the spatiotemporal evolution of each individual phase and tumor growth over time.