Mathematical formulation and parametric analysis of in vitro cell models in microfluidic devices: application to different stages of glioblastoma evolution

In silico models and computer simulation are invaluable tools to better understand complex biological processes such as cancer evolution. However, the complexity of the biological environment, with many cell mechanisms in response to changing physical and chemical external stimuli, makes the associated mathematical models highly non-linear and multiparametric. One of the main problems of these models is the determination of the parameters’ values, which are usually fitted for specific conditions, making the conclusions drawn difficult to generalise. We analyse here an important biological problem: the evolution of hypoxia-driven migratory structures in Glioblastoma Multiforme (GBM), the most aggressive and lethal primary brain tumour. We establish a mathematical model considering the interaction of the tumour cells with oxygen concentration in what is called the go or grow paradigm. We reproduce in this work three different experiments, showing the main GBM structures (pseudopalisade and necrotic core formation), only changing the initial and boundary conditions. We prove that it is possible to obtain versatile mathematical tools which, together with a sound parametric analysis, allow to explain complex biological phenomena. We show the utility of this hybrid “biomimetic in vitro-in silico” platform to help to elucidate the mechanisms involved in cancer processes, to better understand the role of the different phenomena, to test new scientific hypotheses and to design new data-driven experiments.


Methods
Experimental design. We refer to previous works in our group 45,46 for a further explanation of the details of the experimental design. Briefly, human glioblastoma U251-MG cell line, purchased from Sigma Aldrich, was cultured in high glucose Dulbecco's modified Eagle's medium (DMEM) (Lonza, BE12-614F), supplemented with 10% foetal bovine serum (FBS) (Sigma, F7524), 2mM L-glutamine (Lonza, 17-605C) and penicillin/streptomycin (Lonza, 17-602E). U251-MG cells were stably transduced with green fluorescent protein (GFP)-expressing lentiviral vector, kindly provided by Dr. Prats, University Paul Sabatier, Toulouse, France 47 . Shortly, cells were incubated in 1:1 mixture of lentivirus suspension and Opti-MEM medium (Thermo, 31985062) supplemented with 5 μg/mL protamine sulfate (Sigma, P4505). After 24 h, the transduction medium was replaced with growth medium and the cells were routinely cultured for 2 weeks to remove the viral particles. Transduction efficiency was checked by fluorescence microscopy with more than 90% of the cells found to be EGFP-positive.
In order to form a 3D structure, oxygen impermeable microfluidic devices (BEOnChip Ltd.) consisting of a central chamber and two lateral microchannels were used (Fig. 1). They had different dimensions and were made of SU-8, polystyrene or cyclic olefin polymer (COP), using different fabrication processes 45,46 . 3D distribution of cells was achieved within the central chamber, using collagen hydrogel. To prepare 100 µL of collagen hydrogel mixture with a 1.2 mg/mL final collagen concentration, 31.66 µL of 3.79 mg/mL collagen type I from rat tail (Corning, 354236), 0.79 µL of NaOH 1N (Sigma 655104), 10 µL of DMEM 5x (Sigma D5523), 7.55 µL of sterile distilled water and 50 µL of cell suspension were mixed on ice. The mixture was well resuspended and injected into the central chamber of the microfluidic device using a micropipette. The hydrogel droplet was placed on the top of the inlet to prevent evaporation. The devices were placed into an incubator (37 • C, 5% CO 2 ) for 15 min to promote collagen gel polymerization. Afterwards, pre-warmed growth medium was perfused through the lateral channels, mimicking blood vessels, and refreshed every 24h.
Laser confocal and fluorescence images were acquired at different focal planes within each microdevice using a Nikon Eclipse Ti-E C1 confocal microscope. Images were analysed using Fiji software (http://fiji.sc/Fiji). Fluorescence intensity was quantified, in accordance with the software instructions, by selecting a rectangular region across the central microchamber after creating the SUM projection image. In order to transform fluorescence intensity into cell concentration, the cell concentration is assumed proportional to the fluorescence intensity. The constant of proportionality is calculated assuming that the integral of the initial cell concentration along the chamber equals the total amount of cells.
In order to produce the necrotic core formation, a high density of cells ( 40 × 10 6 cells/mL ) was embedded in the collagen hydrogel and injected within the central microchamber. Growth medium was refreshed every day and the culture was maintained for 6 days. Nutrients and oxygen are not able to reach the central part of the device due to cell consumption close to the microchannels, thus causing cell death in the central region appearing an autoinduced necrotic core (Fig. 2), mimicking the parts of the tumour far from functional blood vessels 46 . Visualisation of the necrotic core was performed by calcein/propidium iodide (CAM/PI) staining. Stock solutions of 1 mg/mL CAM (Life Technologies, C1430) and 2 mg/mL PI (Sigma P4170) were diluted to 2 and 4 µg/mL , respectively, in phosphate-buffered saline (PBS) (Lonza BE17-516F). CAM/PI solution was perfused through the lateral microchannels and incubated for 15 min. CAM becomes fluorescent once it reaches the cytoplasm of viable cells and PI stains dead cells, with destroyed membrane.
To promote pseudopalisade formation, cells were seeded at a low density ( 4 × 10 6 cells/mL ) within the central microchamber and one lateral channel was blocked, while constant medium flow was perfused through the other lateral channel. As the region next to the sealed channel was hypoxic, cells migrated towards the perfused channel, rich in nutrients and oxygen (Fig. 3). In the control device, both lateral channels were left open and migration was not observed 45 .
Finally, in the case of double pseudopalisade formation, cells were seeded again at low density ( 4 × 10 6 cells/mL ) within the central microchamber. In this case, the medium was perfused through both www.nature.com/scientificreports/ lateral channels and refreshed every day during 21 days. Hypoxic conditions in the centre of the microchamber induced cell migration towards the perfused channels and invasion of both of them (see Fig. 4).

Mathematical methods. Model equations.
Next, the mathematical model used here for modelling GBM evolution is presented. Even though the mathematical model can be defined for general multidimensional regime, due to the typology of the experiments and for simplicity, the problem may be approximated as one-dimensional, disregarding differences along the direction parallel to the lateral channels. We consider two cell phenotypes (dead cells and alive cells) interacting in the microfluidic device with one chemical species, i.e. oxygen, acting as a regulator of cell processes. These assumptions come from previous experiments in our group 46 that showed that the distribution of other nutrients (glucose) is not responsible of changes in the cells configuration, being oxygen the main (and almost unique) stimulus for cell changes. The variable defining the number of cells for each population at each point and time is their respective concentration u i = C i ( cells mL ), i = 1, 2 , where i = 1 for alive cells and i = 2 for dead cells. Similarly, we call u 0 = O 2 the continuum field of oxygen concentration (in mmHg ). Thus, we shall denote by u = (u 0 , u 1 , u 2 ) T the vector of field variables with 3 rows. The master equation that regulates each variable evolution is the transport equation including source terms: Necrotic core formation. U251 MG cells were seeded at a concentration of 40 × 10 6 cells/mL within the central microchamber. Growth medium was perfused every day through the lateral channels. Viable cells were stained green with calcein AM and dead cells were labelled red with propidium iodide (modified from previous work 46 ). Scale bar is 400 µ m. To mimic thrombotic conditions, medium flow was enabled to flow only through the right microchannel. Under unrestricted conditions, medium was refreshed once a day, through both lateral microchannels (modified from previous work 45  www.nature.com/scientificreports/ with f i the flow term that will include diffusion and chemotaxis for cells and diffusion for oxygen and s i the source term associated to production (proliferation) or loss (death of cells and consumption of oxygen). Note that Eq. (1) is, in general, nonlinear and should embed the coupling between the evolution of the different cell populations regulated by the oxygen concentration that may influence proliferation, migration and death, and oxygen consumption kinetics. In our case, the cell flow term depends on the "random" movement of cells, only driven by differences in their concentration, that is, a diffusion term, and the chemotaxis induced by differences in oxygen concentration (oxygen gradient). For the oxygen pressure, only the diffusion was considered. Then, where D i = D i (u) is the diffusion coefficient and K i = K i (u) the chemotaxis coefficient for population i. Recall that K 0 = 0 . In Eq. (2), we have considered that both coefficients may depend, in general, on the local densities of each cell phenotype and the local concentration of oxygen.
Taking into account the two effects already mentioned for the source term of proliferation and differentiation in cell concentrations and consumption in oxygen, we can write: where α j = α j (u) is the oxygen rate consumed by the cell population j, τ ii = τ ii (u) is the characteristic time of proliferation for population i and τ ij = τ ij (u) the characteristic time of differentiation from population i to population j, that, again, and in general, may depend on the chemical conditions, as well as the local densities of cells. Recall that the apoptotic or necrotic processes are included here as specific differentiation types to the specific phenotype of dead cells.
Equation (1) has to be complemented with the corresponding boundary conditions. We assume here the general case of Robin-like boundary conditions, that is: In the previous equation, x * = 0, L , where L is the width of the chamber. I i = I i (x * , t) and J i = J i (x * , t) are functions characterising the boundary permeability to cell movement or oxygen flow through the boundary, and g i (x * , t) and h i (x * , t) functions defining the controlled value of cell or oxygen concentration and flux at the boundaries. Note that, if I i = 1 and J i = 0 , we have Dirichlet boundary conditions (cell population concentration prescribed at the boundary) and, if I i = 0 and J i = 1 , we have Neumann boundary conditions. Finally, the initial conditions for oxygen and each cell population concentration have to be defined: where u 0 i (x) is a known function. In order to particularise the general equations presented for modelling the population and species evolution in the in vitro experiments made on GBM cells, it is necessary to choose a functional relationship between the coefficients of the model, that is, D i , i = 0, 1, 2 , K i , τ ij , i, j = 1, 2 and α j , j = 1, 2 , and the field variables u.
Even though some papers consider three 48 or four 45 cell phenotypes, here, only two phenotype populations (alive and dead cells) have been considered, thus disregarding possible changes in phenotype along the duration of our experiments. This does not mean that all cells in the chamber equally proceed in terms of proliferation, migration or oxygen consumption, since all these processes depend as well on the particular conditions of the surrounding environment, but that all cells respond equally when they are subjected to the same local environmental conditions. We consider with this assumption that cell adaptation requires longer periods under stressing conditions to modify permanently their internal machinery. Another reason for this assumption is that, in absence of gene expression techniques, it is impossible to distinguish between differentiation into a different phenotype or a change in the cell behaviour as reaction to environmental changes, so, considering one single phenotype for alive cells results in fewer parameters and a better understanding of the role of the different phenomena and parameters, an easier calibration and less uncertainty.
In our microfluidic device, with a controlled production of the hydrogel, we can assume that it is homogeneous and remains with the same properties all along the experimental or, alternatively, that the potential changes in those properties do not affect significatively the cell properties nor the oxygen diffusivity. For alive cells, migration is split in oxygen mediated chemotaxis and pedesis. Dead cells are considered as an inert population ( D 2 = K 2 = 1 τ 21 = 1 τ 22 = α 2 = 0 ). Besides, growth and death rates are also assumed to be dependent on nutrients and oxygen environment. With all these assumptions, it is possible to consider a functional dependency for the following parameters: (2) (3) Scientific Reports | (2020) 10:21193 | https://doi.org/10.1038/s41598-020-78215-3 www.nature.com/scientificreports/ Functions F are nonlinear corrections for cell growth, cell death and oxygen consumption kinetics, while functions model how the different cell mechanisms are activated depending on the oxygen level, D O 2 is the oxygen diffusion coefficient, D n is the diffusion of the normoxic cell population coefficient, χ is the normoxic cell population chemotaxis coefficient, τ g is the characteristic proliferation time, τ d is the death characteristic time and α is the oxygen consumed per unit time and cell. Since cell populations adapt their behaviour to oxygen supply and space availability, two major corrections should be considered in the migration term: • Cellular motility is only possible when the surrounding tissue is not cell saturated 49 .
• Migration following the oxygen gradient happens only when the oxygen supply is below a critical threshold, activating the cell motility mechanism 50,51 .
According to these two major assumptions, a rectified linear unit (ReLU) kind activation function 52 was here used to take into account each of the two phenomena, so the chemotaxis corrections may be written as: with where θ is a threshold parameter.
Here O H 2 is the hypoxia-induced migration activation threshold, representing the oxygen level below which cell migration is activated and C sat is the cell saturation concentration.
The proposed model is in line with the go or grow dichotomy established in GBM literature 53 . Cell energetic resources are spent either in cell migration or in cell proliferation. However, cell proliferation also depends on other needs as nutrient supply or availability of space to grow and split. According to this, we propose a model combining logistic growth and the go or grow paradigm based on oxygen supply. We define the growth corrections as: with and ρ is the logistic correction factor: The function gr is responsible of the go or grow dichotomy and the second is the logistic model for cell population growth.
Cell death is a natural process depending on many factors and agents and has an inherent stochastic nature 54 . Anoxia is one fundamental cause of cell death 55 . Here, a two-parameter sigmoid model is used, able to capture necrosis and apoptosis phenomena: where θ is a threshold parameter and �θ is a sensitivity parameter. They can be seen as a pair of location-spread parameters summarising the stochastic behaviour of the considered phenomenon. With this notation: Finally, oxygen consumption is a complex phenomenon related to the oxidative phosphorylation that occurs in the membrane of cellular mitochondria 56 . Many authors have considered a zero-order consumption function, i.e. a constant consumption rate independent of oxygen concentration O 2 [57][58][59] . A more realistic assumption is to describe the consumption function using the Michaelis-Menten model for enzyme kinetics 48,60 . This model is a correction of the linear consumption and states that: where K is a model parameter. This type of equation was observed for the oxygen consumption rate in the late 1920s and early 1930s 61 . This equation describes more accurately the consumption at low oxygen concentrations and is compatible with previous constant consumption rate models, thus allowing the possibility of comparison with previous studies.
Using this notation and the one introduced in our mathematical formulation, we can write: where 2 is the oxygen concentration at which the reaction rate is half of the rate in a fully oxygenated medium, therefore related to the oxidative phosphorylation kinetics, and the cell structure and morphology (size and number of mitochondria, etc.) and the diffusion process in the cytoplasm.
In the microfluidic device, the culture chamber is connected to the oxygen supplying channels by means of small cavities. The volume and the number of these cavities depend on the microfluidic device design and they are directly related with potential cell losses during the experiment. Actually, when cell populations arrive to the interface between these cavities, some of them may reach the channel and leave the culture. To take into account this phenomenon, we have considered Robin boundary conditions. In principle, since both sides have the same design (number and width of the interface microcavities) there is no reason for considering differences in cell losses (in percentage) between both sides. Therefore, as there is no cell supply through the lateral channels the boundary condition writes: With regard to the dead cell population, homogeneous Neumann boundary conditions are considered since this population does not migrate neither by diffusion nor chemotaxis, so we have Regarding the oxygen supply, we shall consider two possibilities, associated to two different conditions: when oxygen is supplied normally, Dirichlet boundary conditions are considered, that is, we shall assume that the oxygen concentration at the channels remains constant and known throughout the experiment, that is: where O * 2 is a known value. On the other hand, when a channel is sealed, we assume that oxygen provision is negligible, so Neumann boundary conditions are considered: Finally, we assume that, at time t = 0 , all cells are alive and the cell population concentration is known throughout the whole culture chamber. That is, is known (measured experimentally) and C 2 (x, t = 0) = 0 . Moreover, the oxygen profile is assumed to be constant along the chamber and equal to the concentration in the channels, due to the small characteristic time of oxygen diffusion within the hydrogel compared to the characteristic time of cell processes: The differential equation (1) with boundary (5) and initial (6) conditions results in a nonlinear parabolic differential equation in time, with only one space dimension. This equation was solved here numerically by means of a time-space integrator based on a piecewise nonlinear Galerkin approach which is second-order accurate in space 62 , and compatible with this kind of nonlinear equations and boundary conditions. The domain length (associated with the microfluidic device) and mesh size used for the simulation of each experiment are summarised in Table 1.

Results
Parametric analysis. In this section, we discuss the values used in literature for each of the parameters in our model (described in Methods section). We found that many of them are essentially unknown or with high ranges of variation. Our effort goes in the direction of discriminating which works define some of these www.nature.com/scientificreports/ parameters in similar conditions and trying to identify the most likely values within the intervals identified in the literature.
Bibliographic review. In the available literature it is difficult to identify the precise values of such parameters, due to the diversity of models and experimental conditions. Consequently, we include this review clarifying the process for the parameters definition or calculation, often after a reference-crossing process. Cell diffusion coefficient ( D n ). The cell diffusion coefficient is a parameter related to the undriven cellular motility. Cell motility is frequently evaluated in experimental works from a global point of view, that is, including random motility and hypoxia-induced chemotaxis. In this work, however, both phenomena are taken into account separately so diffusion acts as a pure regularisation term while chemotaxis is the main driving force in cell migration. Therefore, only diffusion coefficients associated with healthy tissues in perfect oxygenation conditions will be taken into account.
According to Tija et al. 63 , this parameter depends on the substrate mechanical properties. For a standard collagen ECM, similar to the culture hydrogel here used, a value of 1 × 10 −9 cm 2 /s is proposed. Martínez-González et al. propose in one of their works 48 a value of 6.6 × 10 −12 cm 2 /s , while in another 64 a value of 5 × 10 −10 cm 2 /s is assigned, one order of magnitude lower than the mean of the values reported by Rockne et al. 65 ( 5 × 10 −9 cm 2 /s ). Wang et al. 66 discuss this value for different locations in the brain, observing that glioma cells migrate quicker in white matter than in grey matter, highlighting also the important variation of this coefficient with the tumour stage and after chemotherapy and radiotherapy, ranging all values from 3 × 10 −7 cm 2 /s to 5 × 10 −5 cm 2 /s (median of 3 × 10 −6 cm 2 /s ). Hathout et al. 67 , use values from 5 × 10 −7 cm 2 /s to 2 × 10 −6 cm 2 /s during the tumour initial state.
Chemotaxis coefficient ( χ ). This coefficient is difficult to estimate when considering chemotaxis as an isolated phenomenon 68 . Ford et al. 69 with χ 0 ranging from 1.5 × 10 −5 cm 2 /s to 7.5 × 10 −4 cm 2 /s depending on the complex affinity while several expressions are proposed for f. For example, a hyperbolic tangent dependence is presented, based on a probabilistic mechanobiological model for individual bacteria 70 .
Many other works define the chemotaxis coefficient with respect to the normalised 2 is the vessel oxygen supply pressure. Agosti et al. 68 assume a value of 1.5 × 10 −4 cm 2 /(mM · s) for an oxygen concentration in vessels O v 2 of 0.07 mM 71 . Therefore, χ is of the order 2 × 10 −7 cm 2 /(mmHg · s) assuming an oxygen supply in vessels of 40 − 60 mmHg 72, 73 . With the same conversion between oxygen concentration and pressure, a value between 3 × 10 −10 cm 2 /(mmHg · s) and 1 × 10 −9 cm 2 /(mmHg · s) is adopted by Agosti et al. in their work 74 . Finally, Bearer et al. 75 propose a chemotaxis coefficient of 10 5 µm 2 /d , which gives an equivalent value of 2 × 10 −10 cm 2 /(mmHg · s).
Hypoxia-induced migration activation threshold ( O H 2 ). In our model, hypoxia induced cell migration is relevant only when the oxygen pressure is under a certain threshold O th 2 . According to previous works on GBM simulation 48,64 , cells are considered under hypoxia conditions, when the oxygen pressure is under 7 mmHg (approximately 12 − 18% of the blood vessel oxygen pressure). Agosti et al. consider in one work 68 a threshold of 15 − 50% of blood vessel oxygen pressure and later 74 a threshold of 30% is used. In the review paper 76 , a ratio of 12 − 25% between healthy and tumorous tissue oxygen pressure is considered.
Growth characteristic time ( τ g ). This is also a very context-dependent parameter, since the cell metabolism highly varies between cell types and individuals. In addition, our proposed logistic model implies that the measured growth time in the particular experimental conditions depends on the cell concentration, and therefore could vary with considered values reported in literature. Nevertheless, some growth characteristic times reported in literature for logistic, exponential or Gompertz growth models are here discussed. Gerlee et al. 77  In our work, based on the go or growth assumption, the growth characteristic time is infinite in absence of oxygen and decreases until the oxygen concentration exceeds the hypoxia threshold. Thus, our model captures this variability from hypoxic to normoxic media, where growth is accelerated and therefore characteristic times are smaller.
Cell concentration saturation ( C sat ). An important variability is found in the literature when referring to this parameter, with a range that covers several orders of magnitude. For example, Rockne et al. 65 , propose a value of 10 11 cell/cm 3 whereas Hathout et al. 67 use the value of 10 8 cell/cm 3 according to previous experimental works 85 . www.nature.com/scientificreports/ This parameter depends on the mechanical and structural properties of the medium and on nutrients supply so its variability is natural. In any case, it does not have a major impact in simulations for cell concentrations much lower than the saturation capacity.
Death characteristic time ( τ d ). Even in the case where no cell-concentration dependence is considered, death characteristic time also varies between studies since it is directly measured, without considering, for example, oxygenation conditions, as for the growth characteristic time. In the automaton model from Gerlee et al. 77 , an average apoptosis probability of 0.18 is obtained, resulting in a death characteristic time of 72 h as proposed by Frieboes et al. 80 . In their works, Agosti et al. propose values between 160 h and 400 h 68 , and of 600 h 74 . Finally, Martínez-González et al. 48,64 use two different phenotypes to model the tumour population (normoxic and hypoxic), but they assume that once the cell has arrived to hypoxic conditions, its death characteristic time is 48 h 48 or 7 d 64 .
We model death with a sigmoid function, which integrates both death causes: apoptosis, that is mainly stochastically mediated and necrosis, induced by oxygen lack. This model explains better the variability found in literature, ranging from 72 h in anoxia to 600 h in normoxia, via the two parameters regulating cell switch, discussed below.
Anoxia-induced death location parameter ( O A 2 ). Even though in many mathematical models it is assumed that the hypoxia threshold, inducing migration or proliferation (and therefore the fundamental parameter explaining the go or grow dichotomy), and the anoxia threshold (as an indicator of necrosis) are the same 68 Oxygen diffusion coefficient ( D O 2 ). The oxygen diffusion coefficient is classically known to be around 10 −5 cm 2 /s at 37 • C . Daşu et al. 87 propose a value of 2 × 10 −5 cm 2 /s according to previous studies 57,58 that assign an intermediate value between oxygen diffusion in water and muscle 88 . Other works 77,89 use a value of 1.8 × 10 −5 cm 2 /s . Recent computational patient-specific studies 74 assume a value of 10 −5 cm 2 /s . It is important to note that, in the present work, hydrogels used in microfluidic devices try to reproduce soft human tissue, so similar values can be used.
Oxygen consumption rate ( α ). The maximum value of α is much debated 76,90-92 ranging from 2 µL/(min · g) 76 to 55 µL/(min · g) 90 . There are several possible explanations for this large range of values reported as explained by Daşu et al. 87 , such as the influence of the tissue metabolic characteristics in the consumption rate, the variations of the temperature and pressure conditions when measuring the oxygen volume or experimental reasons associated with the measuring method. The most often quoted value is 15 mmHg/s for the maximum consumption rate in healthy tissue 57,87,93 . This consumption rate gives a maximum diffusion distance of 143 µm 93 , for a blood vessel with 40 mmHg . Assuming that a healthy tissue has a concentration of 0.2C sat 64 , we obtain 7.5 × 10 −7 (mmHg cm 3 )/(cell s) . Assuming the same ambient cell concentration, using the value proposed by Agosti et al. 68 96 . This value has been chosen in recent simulation models 48,64 .
In Table 2 all numerical parameters of the mathematical model are shown with their corresponding variation range extracted from the bibliography.
Model parameter fitting. The value of each parameter was initially chosen to stay within the ranges found in the bibliography ( Table 2). In order to calibrate the specific value for each parameter, the formation of the necrotic core experiment (described in Methods section) was selected as case study. Robin boundary conditions were chosen allowing alive cells to escape from the device. In particular, the following values were chosen for the boundary conditions (5), according to the experimental results and the symmetry of the experiment, we select K 1 (x * = 0, L, t) = 1 , g 1 (x * = 0, L, t) = h 1 (x * = 0, L, t) = 0 and J 1 (x * = 0, L, t) = 1.0 × 10 6 s/cm . Only the parameter J 1 , explaining cell leakage at boundaries, was fitted to capture the corresponding results. For dead cells, homogeneous Neumann conditions were established, assuming no migration of dead cells through the boundaries, that is, K 2 (x * = 0, L, t) = 0 , h 2 (x * = 0, L, t) = 0 and J 2 (x * = 0, L, t) = 1.0 s/cm . With respect to oxygen concentration, Dirichlet boundary conditions were chosen, so oxygen pressure in the channels was assumed to remain constant throughout the whole experiment, since medium flow provision through the channels was sufficiently frequent to keep that pressure without important variations despite oxygen diffusion and cell uptake. With that, we write K 0 (x * = 0, L, t) = 1 , g 0 (x * = 0, L, t) = 7.0 mmHg and J 0 (x * = 0, L, t) = 0 . Finally, as initial conditions, we assume that the initial monitoring of the process starts after getting a uniform oxygen pressure in the whole chamber, equivalent to the one present in the lateral channels, that is O 2 (t = 0) = 7 mmHg.
A heuristic approach was followed to fit the simulated curves with the experimental results, in order to determine the model parameters. This approach tried to get the best fit of the necrotic core (central region) due to its biological relevance, giving less importance to the fitting around the boundaries since the cell distribution here is not representative of the in vivo situation. To quantify the quality of this fitting procedure, we defined two objective cost functions: where u e j is the experimental measurement of the j phenotype ( j = 1 , alive cells, j = 2 , dead cells), u j the simulated one and L the chip length.
After the fitting process for the different results got in the necrotic core experiment, we arrived to the value set for the parameters shown in Table 3, yieliding the results shown in Fig. 5. The fitting process provided values of T = 0.17 and D = 0.73, which reinforce the good agreement between the experimental and simulation results.  www.nature.com/scientificreports/ As it may be observed, the computed results are qualitatively equal and quantitatively very similar to the experimental ones, although there are some significant discrepancies in the alive cell profile, mainly at the centre of the chamber and in the dead cell profile at the boundaries. These differences are unavoidable due to the effects and interactions missed in the model, such as the heterogeneous distribution of hydrogel, cells and oxygen in the initial state and boundary conditions, and to the highly non-linear character of the equations. As it can be seen in Table 2, the range of variation of the parameters in the bibliography is extremely wide in some cases, being therefore tricky to tune the value of the parameters to obtain better numerical results. All this makes essential to perform a sensitivity analysis to understand the quantitative impact of each parameter on the representative results and to assess the mathematical model robustness with respect to the parameter fitting.
Sensitivity analysis. This sensitivity analysis is performed subsequently, with the aim of defining a robust methodology capable of providing a sufficiently accurate reproduction of different in vitro experiments, while avoiding the common overfitting in biological problems. D n , D O 2 , α , χ , τ g and τ d were the parameters studied. The threshold parameters were not considered in this analysis since they mainly affect the location and spread of the necrotic core. Each parameter, say θ , was individually perturbed as follows: if U is a uniform U([−1; 1]) random variable, we generated 100 samples of � = θ × 10 U , that is, θ was perturbed but maintaining its magnitude order. The non-perturbed value of θ corresponds to the value obtained after the fitting procedure. Percentiles 5, 50 (median), and 95 were extracted from the resulting concentration profiles. All results are shown in Figs. 6, 7 and 8. According to the sensitivity analysis performed, a final range of variation for D n , D O 2 , α and χ was chosen (see Table 3), since they are the parameters with a higher impact on the computed results while presenting important variability in the bibliography. The interval bounds were selected to guarantee a value of T < 0.15 for the best simulation in each experiment. τ g and τ d do not have a variation range as their impact on the results is low.
For the simulations presented from here to the end of the paper, each parameter was considered as a random variable with � i ∼ U[θ 1 i ; θ 2 i ] , with θ 1 i and θ 2 i in agreement with the previous discussion. To estimate the correlation coefficients, a set of 100 simulations were performed varying the parameters within their order of magnitude (as done in the sensitivity analysis), and the sets that provide solutions with T < 0.2 were kept. Then, the Kendall correlation coefficient was computed, obtaining a value of τ = 0.9 for D O 2 and α and a value of τ = 0.2 for D n and χ . The rest of the parameters are supposed to be mutually independent, and therefore uncorrelated.
A more in-depth analysis of this multiparametric correlation and its effects on the results is possible, but is out of the scope of this work. Here, the aim is to illustrate that, from the sensitivity analysis, it is possible to get a higher insight into the mathematical model that could be taken into account when performing Montecarlo simulations.
In silico replication of the in vitro experiments. We analyse the performance of the mathematical model presented in Methods section when using one single set of parameters (Table 3), applied to the three experiments described in Methods section: formation of a necrotic core in high concentrated cultures, pseudopalisade formation due to oxygen gradient and double pseudopalisade formation in a symmetric configuration. As there is an inherent uncertainty in the parameters identification, a run of 100 Montecarlo simulations was performed for varying values of the model parameters according to the ranges defined in Table 3 and the parameter correlations discussed in the previous section.
The cells boundary conditions remained the same in all the experiments except for the value of J, which depends on the cell leakage observed for each microfluidic device. In experiment 1, J = 1.0 × 10 6 s/cm ; in experiment 2, J = 1.0 × 10 9 cm/s ; and in experiment 3, J = 1.2 × 10 7 cm/s . Regarding the oxygen boundary conditions, they were adapted to each experimental configuration: in the formation of the double pseudopalisade they were identical to those already explained for the necrotic core formation; whereas in the pseudopalisade formation, impermeability condition (no flux) was imposed at the sealed channel, while, again, the Dirichlet condition was imposed at the right channel. The value of the oxygen pressure, both at the right side and at the initial time for the whole chamber, was fixed to O * 2 = 2 mmHg instead of O * 2 = 7 mmHg as in the other experiments. This is justified by the fact that in this experiment the medium was not renewed as in the previous cases, so the oxygenation was assumed to be lower.
The obtained results for each experiment together with the 90% confident band (between 5th and 95th percentile) of the simulations, are shown in Figs. 9, 10 and 11.
The figures show good agreement between the experimental and simulated results, which will be further explained in Discussion section. Moreover, considering the parameter variability improves the estimation of the dead cell profile (Fig. 9) .

Discussion
Glioblastoma is one of the deadliest tumour types, as it is very heterogeneous and resistant to therapy [97][98][99] . Most research studies have been performed in 2D models, on Petri dishes, but these are not able to represent the real situation. Many 3D culture models are now being developed, such as spheroids, organoids, different scaffold-based models, etc., which better mimic cell-cell and cell-extracellular matrix interactions 100,101 . With the development of microfluidics for biological applications, apart from including different components of the tumour microenvironment, we are also able to control the physico-chemical conditions, so microfluidic models are considered the most biomimetical in vitro models nowadays 102,103 . Special devices have been developed for GBM research to study the behaviour of GBM cells and therapy response within a biomimetic and controlled microenvironment. With them, realistic behavioural patterns have been obtained, similar to the ones Scientific Reports | (2020) 10:21193 | https://doi.org/10.1038/s41598-020-78215-3 www.nature.com/scientificreports/ observed in patients 46,[104][105][106] . In our case, we were able to reproduce autoinduced necrotic core and pseudopalisade formation 45,46 , some of the most important characteristics of GBM, in agreement with the GBM formation model of Brat et al. 107 , which identifies blood vessel occlusion and subsequent hypoxia-induced migration towards the functional blood vessel, following oxygen and nutrients gradient, as one of the main drivers of GBM invasion. From our models of pseudopalisades and necrotic core formation in GBM in vitro models, we set the following phenomena in the model mathematical equations: • GBM cells migrate in an O 2 pressure gradient following a chemotactic cue from lower to higher pressure values, and with a migration speed that depends on the specific local O 2 pressure value. • Very high cell concentrations prevent the arrival of sufficient oxygen to regions far from oxygen provisions sources, due to the oxygen uptake in the transition zones, thus resulting in an auto-induced necrotic core in those far regions. Lee et al. 108 explained the importance of shortage of oxygen and nutrients in necrotic core www.nature.com/scientificreports/ formation. Also, a similar process was observed in spheroid cultures 109,110 , where the gradient depends on the spheroid diameter, observing the appearance of necrotic core in spheroids bigger than 500 µm.  www.nature.com/scientificreports/ -The effect of growth and death times, τ g and τ d respectively, have an obvious impact over the duration and speed of the phenomena involved. However, as the chamber is essentially under hypoxic conditions, the death time plays a major role on the density of the necrotic core (but not in its size).
• Cell adaptation may also have an important role in problems like the one described here but has not been considered so far in our work. This will be part of future developments, although it will require new and specific sets of experiments to capture the corresponding mathematical features and parameters.
Since the model parameters fitting was established under a heuristic basis according to our research experience, the values selected should be interpreted qualitatively as the ones which better describe the most relevant phenomena that take place in glioblastoma evolution in vitro, such as necrotic core formation. It is important www.nature.com/scientificreports/ to remark that a preliminary fitting approach based on a formal mathematical optimization did not provide the best fit in the evolution of the necrotic core (data not shown), demonstrating the difficulty of the problem. Despite the parameter uncertainty, it has been possible to reproduce three different experiments with one single set of parameters. Therefore, it seems possible to extract fundamental biological conclusions such as: • The initial cell density has a crucial influence on the necrotic core formation. The simulation results (and the experimental curves) related to the necrotic core and the double pseudopalisade experiments are essentially identical except for the fact that the former was obtained with an initial concentration of C 0 = 40 × 10 6 cells/mL and the latter with C 0 = 4 × 10 6 cells/mL . The corresponding results are however qualitatively very different: the necrotic core only appears when having very high cell concentrations. This conclusion may have important biological and therapeutic consequences. • Cell migration depends strongly on the oxygen level and gradient. Suitable oxygenation conditions do not induce cell migration even for high oxygen gradients. Conversely, under a certain oxygenation level, the oxygen gradient is the main driving agent of cell migration. It has been shown 111,112 that hypoxic conditions lead to stabilization of hypoxia inducing factor (HIF) which regulates many important pathways important for tumour progression, such as invasion and angiogenesis.  www.nature.com/scientificreports/ The presented mathematical model is, therefore, able to capture some of the main features of some essential phenomena occurring during GBM invasion. Moreover, except for the saturation parameter C sat (that is obviously very dependent on the experimental conditions) and for some values of the oxygen consumption α , the parameter variability is in agreement with that found in the scientific literature. Also, most of the general features observed are similar to the ones obtained in previous experimental 45,46 and computational 41,45,48,64,71 works. However, the model also presents some limitations, some of the most important are the following: • Small discrepancies between the experimental and computational results were found. One possible explanation is that the accumulation of cells at the boundaries may obstruct the oxygen diffusion, provoking a nonhomogeneous O 2 diffusion coefficient. This results in an over-estimation of the oxygen level in the central area  www.nature.com/scientificreports/ of the chamber, which could explain the over-estimation of alive cells, faster cell migration to the boundaries and over-estimation of dead cells. • There is, in general but mainly in the results associated to the necrotic core experiment, a certain lag between the computational and the experimental responses of cells to oxygen variations. It seems that changes in the oxygen concentration are felt earlier in the simulations. This may be associated with a certain cell memory: the cell may need to accumulate some damage before undergoing significant changes in its behaviour. Our GBM model is not able to capture this kind of phenomena. Nonetheless, we have presented a framework where this limitation can be overcome if more phenotypes are considered as in other works 45,48,64 . This strategy requires, however, a sound classification of the cell phenotypes, based on biological evidence, in a sufficient number of classes, which is difficult and would considerably increase the number of parameters.

Conclusions
From the results and discussion presented above, we enumerate the main findings and conclusions of the paper: 1. Mathematical modelling and the corresponding computer simulation of complex cell processes, incorporating cell interactions, chemical and physical cues, require an extensive literature review and the analysis of the fundamental properties of the mathematical equations in simplified conditions, and an in-depth analysis of the model parameters, in order to understand the individual and combined effect of each combination of parameters, both qualitative and quantitatively, in the resulting variables. 2. One single type of experiment is not enough to guarantee a proper quantification and understanding of the effect of each model parameter. Some families of experiments have to be used to fit the parameters, while other families are required to validate and discard spurious parametric solutions. This strategy is fundamental to avoid overfitting and to prevent model-induced effects, result of the fitting procedure. 3. There is a huge variation in the range of many of the parameters found in the literature, sometimes with intervals covering several orders of magnitude, which makes it very difficult to get a reasonable accuracy when modelling experimental tests with the only use of values from bibliography. This can be a result of the high heterogeneity of GBM and of the high adaptive capacity of these cells 113,114 . 4. With a proper parameter identification and analysis, if all the main mechanisms involved are properly considered, it is possible to get an accurate reproduction of experimental tests, provided the experiment is well controlled, the associated variables are accurately measured and the results are correctly interpreted. 5. A proper parameter sensitivity analysis is essential to discover hidden effects that cannot be explained by the model (e.g. presence of dead cells close to the lateral channels), disregard wrong value intervals (e. g. the range of the parameters was strongly reduced from Table 2 to Table 3) and identify the actual conditions in which the experiment is performed. 6. Adopting the presented model as a starting point, there is still room for future development. For instance, the measurement of the oxygen profile would allow us to improve the oxygen diffusion model, taking into account, for example, the oxygen flow obstruction that may be induced by high cell densities.
To summarise, the presented framework is general and allows the analysis of many coupled and highly non-linear physical mechanisms. The effort made in the parametric analysis allows to draw conclusions both qualitative (e.g. pseudopalisade and necrotic core formation) and quantitative (e.g. time scale for necrosis or speed of migration structures). This task is fundamental when working with complex multiparametric models. Nonetheless, this analysis is always conditioned by the choice of the mathematical approach, so the intrinsic model uncertainty is, to some extent, unavoidable. Working with models with so many parameters requires always enough experimental data in sufficiently varied conditions. There, we have been able to work with three different families of experiments resulting in cell profiles along space and time. This extensive amount of information gives value to our work, which could lay the foundations for future works in the topic. Finally, we can conclude that, even with all this, we are still far from getting sufficient knowledge of all the mechanisms involved in complex biological processes, as well as the interactions and quantification of the corresponding phenomenological parameters. Only in very specific and well-controlled conditions, and after an extensive analysis of the tests, model and associated parameters, it is possible to expect for accurate results if the initial conditions are well-measured and the main mechanisms and interactions are mathematically represented. Despite these drawbacks, mathematical models are today invaluable tools to better understand underlying mechanisms and interactions, to establish trends, to test new hypotheses and to check "what if " situations that are many times impossible to test experimentally due to the impossibility to isolate single effects, measure particular variables or simply for ethical reasons. www.nature.com/scientificreports/