Interactions between ploidy and resource availability shape clonal interference at initiation and recurrence of glioblastoma

Glioblastoma (GBM) is the most aggressive form of primary brain tumor. Complete surgical resection of GBM is almost impossible due to the infiltrative nature of the cancer. While no evidence for recent selection events have been found after diagnosis, the selective forces that govern gliomagenesis are strong, shaping the tumor's cell composition during the initial progression to malignancy with late consequences for invasiveness and therapy response. We present a mathematical model that simulates the growth and invasion of a glioma, given its ploidy level and the nature of its brain tissue micro-environment (TME), and use it to make inferences about GBM initiation and response to standard-of-care treatment. We approximate the spatial distribution of resource access in the TME through integration of in-silico modelling, multi-omics data and image analysis of primary and recurrent GBM. In the pre-malignant setting, our in-silico results suggest that low ploidy cancer cells are more resistant to starvation-induced cell death. In the malignant setting, between first and second surgery, simulated tumors with different ploidy compositions progressed at different rates. Whether higher ploidy predicted fast recurrence, however, depended on the TME. Historical data supports this dependence on TME resources, as shown by a significant correlation between the median glucose uptake rates in human tissues and the median ploidy of cancer types that arise in the respective tissues (Spearman r = −0.70; P = 0.026). Taken together our findings suggest that availability of metabolic substrates in the TME drives different cell fate decisions for cancer cells with different ploidy and shapes GBM disease initiation and relapse characteristics.

8 !-# # , # $ , %/ is a flux vector to the neighboring voxels, defined as: , where, k 30:,! is the maximum rate of cell migration for cells with ploidy i, and < 3!$,!<"4 C 3!$,!describe the biphasic dependence of cell migration on glucose availability as: Y-# $ / = h G0 +!( •( ( -1 − h G; +!( •( ( / … `7 `(%) = c @ ⋅ % @ + c g ⋅ %The `(%) = c @ ⋅ % @ + c g ⋅ % We assume that resource perfusion within the brain tissue under normal conditions is maintained at homeostasis via existing vasculature.However, the presence of rapidly growing tumor cells consumes resources at a rate faster than it can be replenished.At the same time, tumors can recruit additional vasculature to bring in additional resources.Under the simplification that each voxel is small enough that vasculature within that volume evenly distributes resources throughout the voxel volume, but large enough that diffusion of nutrients across voxel boundaries can be neglected (characteristic diffusion length of nutrient reaction-diffusion -~100 µm -is much smaller than voxel dimensions: ~4 mm 2 ).Thus, changes in resources (oxygen and glucose individually) within a voxel are modeled as: In equations 2<,2C, F $ , F # are the additional resource generation rates by recruited vasculature and are the excess resource consumption rates by a voxel of tumor cells.Hereby D $@ ≔ 16 ⋅ D $ , i.e. cells consume more glucose when oxygen is low.Vasculature recruitment probability is given by equation 2c, whereby ¢ tracks adjacent voxels.The recruitment of new vasculature is promoted by a voxel being at certain capacity (1/4 shown here) and surrounded by voxels that are stressed due to comparable ratios of live and dead cells.Before triggering cell death, nutrient starvation elicits an autophagic response, which in turn has been shown to mediate pro-angiogenic effects via AMPK/Akt/mTOR signaling 84 [64].The ROS-ER stress-autophagy pathway has also been shown to induce angiogenesis in and ischemic myocardium model 85 .Equation 2c therefor assume that comparable ratios of live and dead cells in a voxel maximize the probability of angiogenesis.As more dead cells accumulate the efficacy of recruiting vasculature is assumed to decrease.If on the other hand the voxel is dominated by alive cells we assume vascular recruitment is unnecessary as the area is sufficiently vascularized.The presence of recruited vasculature ($) is modeled as a binary value of 0 or 1.We assume radiation therapy will immediately kill all angiogenic cells within the treated voxels 86,87 .
The model does have several floating parameters, particularly in terms describing resource and ploidy dependent cell growth (' !, Γ !(e) ) and death (4 !(e, 5)) rates, as well as tissue stiffness and resource dependent cell migration (8 !(e, %)) rates for cells with different ploidy.An overview of all 25 model parameters is given in Table 1.For 13 of these parameters we found evidence from literature supporting a value or range typical for GBM, or we had in vitro data specific for GBM to inform them (Level 1).Another 12 of these parameters had approximations available for multiple other cancer types or other cell types, covering a wider range.In particular parameters representing resource availability and resource-driven phenotypic transitions were not well described and the motivation behind the selected values is described in the following section.

Motivation for literature-based parameters related to resource-driven phenotypic change
Oth -threshold of oxygen concentration below which a metabolic switch from oxidative phosphorylation (OXPHOS) to aerobic glycolysis occurs.The parameter was established based on threshold O2 concentration (0.5%) for which the induction of expression level of hypoxia-inducible factor-1α and -1β (HIF-1α/β), key cellular regulator of metabolic switch from OXPHOS to glycolysis, is maximal, as measured in cultured HeLa cells 37 ;.
Cd -clearance rate of dead cells via autophagy or related mechanisms.The range for this parameter (0.48-0.78) is motivated by the quantification of macrophage engulfment of target cells triggered to undergo apoptosis, necroptosis, or ferroptosis over 24 hours 28 and by the abundance of tumor-associated macrophages and microglia in the GBM tumor mass 29 .
X# -oxygen-dependent death constant per cell type (O2 concentration at which cell survival is half-maximal) was selected based on the results of study evaluating the proliferation of different cell types under oxygen deprivation in vitro 20 .Oxygen concentration as low as 0.5% did not alter the fibroblast or tumor cells survival, irrespective of HIF-1 status, while extreme hypoxia (<0.01%) led to a block in proliferation and a progressive loss of viable cells to 20%-30% at 72 hours, also irrespective of HIF-1 status.Values in range 0.0007-0.034mg/L (0.01-0.5% at 36°C) are supported by other studies of tumor proliferation under extreme oxygen concentrations 37 .ì# -glucose-dependent death constant per cell type (glucose concentration at which cell survival is half-maximal) was selected based on the results of a study examining impact of low-glucose conditions on the survival of lung cancer cells 21 .We digitized Figure 1 using WebPlotDigitizer 88 and arrived at parameter range between 0-2.9 mMol depending on the cell line and pH conditions.The impact of pH on parameter values was not evaluated in S3MB.
X$ -oxygen-dependent growth constant per cell type (O2 concentration at which cell proliferation is half-maximal) was selected based on the results of study evaluating the proliferation of different cell types under oxygen deprivation in vitro 20 .Similar to X#, it was found that cells continue to proliferate with undiminished kinetics under 0.5% oxygen and slow down proliferation in lower O2 concentrations, as reviewed by McKeown et al 37 .;X$ was assumed to be slightly higher compared to X#.
<Tchemo -oxygen at which chemo-shrinkage rate is half-maximum was selected based on temozolomide IC50 fold difference under hypoxia vs normoxia 34 .c0 -radiation-induced change in stiffnes -was the parameter value that led to a reduction of initial brain tissue stiffness by 55%, consistent with radiation-induced oedema 25 .Specifically, reduction in brain tissue stiffness is assumed to happen incrementally after each radiation cycle, without significant temporal delay 77 .A 55% reduction is achieved upon completion of a standard of care radiation therapy regimen, consisting of 60 Gy administered over 30 cycles.c1change in stiffness from remodelling by tumor cells -Given StiffMax, the maximum GBM tissue stiffness of 22.5 kPa 18 and StiffIC, the average stiffness of healthy brain as quantified in our MRE derived stiffness atlas, we assume: StiffIC *(1+c1) = StiffMax and solve for c1.

Deriving oxygen, glucose and stiffness atlases of the human brain
We derive glucose uptake, oxygen and brain tissue stiffness atlases as follows.

Glucose level atlas
An atlas of glucose-corrected standardized uptake values (SUVs), which inversely correlates with regional glucose availability, is available from a public database of 37 healthy adult human brains.It is comprised of [ 18 F]FDG PET, T1 MRI, FLAIR MRI, and CT images for each brain 47 .The PET scans were already coregistered to their respective T1 anatomical scans.To allow us to create an average brain for the purpose of our model simulations, all brains needed to exist in a standard space in which they could be averaged voxel-wise; to complete this, the scans were aligned with the MNI ICBM 152 2009c Nonlinear Asymmetric template 89 space using a linear followed by a nonlinear transformation with the Analysis of Functional Neuroimaging (AFNI) function 3dQwarp.One aged stiffness scan had to be excluded because of poor registration.
Following co-registration and registration to MNI space, PET images were normalized to obtain SUV and Standard Uptake Value ratio (SUVr) images at a ≤ 5% risk of Type 1 error 47 .These are used here and stratified by the patient's age.To approximate glucose levels, we assumed SUVr values are inversely correlated to glucose uptake 48,49 and that glucose concentration in the brain ranges between 1.0-2.7 mMol 50 .This was achieved through reflection of SUVr values, followed by scaling to the known range of brain glucose concentrations: # $ = `B<ih(−%v$Y, [1.0, 2.7 btli]).Since SUVr values are dimensionless, the scaling also assigns a unit to the resulting estimate # $ .

Oxygen level atlas
The quantitative BOLD (qBOLD) data from seven volunteers was downloaded from the Oxford Research Archive site 51 .All image processing was performed with Quantiphyse 52 .To improve signal-to-noise ratio, four slices of 1.25mm thickness had been averaged together and there was slice gap of 2.5mm between resulting 5mm slices.We modified the NIFTI header to correctly report the actual slice interval of 7.5 mm, which allowed the qBOLD scans to align well with the anatomical scans.A linear registration and motion correction available in Quantiphyse 90 , was used to correct any motion defects between timepoints of the qBOLD data.To suppress isolated noisy voxels we performed sub-voxel smoothing with a smoothing kernel of 1.5 mm.The mean reversible transverse relaxation rate (T2), the oxygen extraction fraction (oef) and the deoxygenated blood volume (dbv) were inferred with Bayesian modelling using the qBOLD widget in Quantiphyse 53 .The method deconvolutes three major confounding effects of qBOLD: (i) Cerebral Spinal Fluid (CSF) signal contamination; ii) Macroscopic magnetic field inhomogeneity and iii) Underlying transverse T2 signal decay 53 .To register all data to the Montreal Neurological Institute coordinate system (MNI) a DEEDS deformable registration 91 was performed to the anatomical images and the registration matrix was then applied to all the qBOLD data.Finally, an atlas of dissolved oxygen concentration (dO2) was derived from dbv and oef, in units of mg/L as: , where O2cap = 0.4655 mg/L is the typical average oxygen partial pressure expected in capillaries 54 .One scan had to be excluded because of poor registration.

Brain tissue stiffness atlas
Magnetic Resonance Elastography (MRE) scans were obtained from a database of 12+12 aged and young brains generated at the University of Edinburgh, UK 12 .The T1 scans were first skull stripped with AFNI's 3dSkullStrip function.
Next, because the MRE scans only covered a portion of the brain, a pad of 15 zeros was applied to the edge of the scan using the 3dZeropad function.The MRE scan was then co-registered to the T1 scan via rigid body transformation with the align_epi_anat.pyfunctions, using the "-giant_move" option and NMI cost function; for those subjects who failed alignment, the "-ginormous_move" option (3 brains) or eliminating the move option (one brain) was used instead.The MRE scans were then transformed to MNI space by the same protocol as the PET scans.

Correlation with age
MRI scan demographics for age and sex are tabulated below in Table S1.
We first considered mean voxel value within the white matter (WM), grey matter (GM), and the whole brain parenchyma (WM+GM).We utilized the FSL FAST automated segmentation tool to create masks of the WM, GM, and cerebral spinal fluid for all of the t1 scans associated with patients who received PET and MRE scans.We segmented with FAST options: 3 classes, main MRF parameter of 0.1, 4 iterations for bias field removal, and bias field smoothing of 20.0 mm FWHM.We masked the various functional scans with the segmentations using afni's 3dcalc function and and computed the voxelwise mean (excluding zero values) with afni's 3dBrickStats.
Because the coregistered MRE scans were somewhat distorted on the edges, incorporating values between 0 and the minimum MRE value, we thresholded the t1-aligned MRE scans with the minimum original MRE value.This prevented the incorporation of low value but non-zero voxels in the grey matter average calculation.The rest of the steps were the same as those for the PET scans.
We then proceeded to compute Pearson correlation coefficients between patient age and brain tissue stiffness, glucose, and oxygen in a spatially dependent manner to determine which brain regions were most significantly altered during aging.Supplementary Figure 6: Map of Pearson correlation coefficients between patient age and brain tissue stiffness (A), Glucose (B) and Oxygen (C) across voxels.

Modeling the effect of surgery on oxygen, glucose and stiffness atlas
Because stiffness, glucose, and oxygen maps all share the same MNI coordinate system as the tumor cell density maps (see Supplementary section 4.3), we can model the effect of surgery on these environmental factors.Let RÉ$ denote the resection cavity map after a given surgery and %, # $ , # # are the stiffness, glucose, and oxygen maps.We assume tumor resection decreases stiffness as well as perfusion in the surgery-affected brain region (due to lack of a fractal pattern in blood supply provided by vasculature), as follows:

Cells
Experiments were performed using the near-diploid (2N) SUM-259 cell line.Near-diploid cells have previously been fused to derive a near-tetraploid (4N) lineage of the cell line as described in Miroshnychenko & Baratchart et al. 92 .All fused (4N) cells contained both, a fluorescent Mcherry staining to mark the nucleus and the GFP fluourescent probe to mark the cytoplasm.

Incucyte transwell assay
DMEM media was prepared from powder (Sigma-Aldrich, 0453-0020) to vary glucose and dFBS concentrations (0.25-25 mM and 0.1-10% respectively).Media was filtered through a sterile 0.2mm filter bottle (Nalgene, 450-0020).Cells were starved for 12-24hrs before testing by replacing media with DMEM containing 1mM glucose and 1% dFBS.The Clearview Cell Migration plate (Sartorius, 4582) was used, if necessary, with a coating protocol with 25mg/mL collagen 1 (Corning, 354236) for an hour at 37ºC.Cells were seeded into the upper chamber at 1,000 cells per well in 60mL DMEM with the appropriate glucose and dFBS concentrations.The lower portion of the wells were all filled with 200mL DMEM with 25mM glucose and 10% dFBS.The plate was imaged using the Chemotaxis assay in a Sartorius IncuCyte S3 at 37ºC & 5% CO2 using 10x magnification, and images were taken every 2h for 48h.

Image analysis
We processed 6000 raw images from the Incucyte transwell assay (80 wells; 25 times points per well; and three channels Phase, Top, and Bottom); for this, we colored the top and bottom channels with red and green respectively, then we merged them and stacked the time points images associated with the same well, thus obtaining 80 files equal to the number of wells (80 wells).
After manipulating the raw images, the newly processed images were used for the creation of the training and test datasets.First, we selected only one frame (time point) per well sequentially, so all different times, from 0 to 24, were sampled throughout the files (wells) at least once.Second, we cropped a random region of 200 × 200px from each frame selected.All this process was done twice to create both training and test datasets.Each dataset was then composed of 80 frames.
Succeeding the data preparation, we used the image classification and segmentation software ilastik 44 to obtain cell counts and cells' spatial information.We trained the model (using the training dataset) by first assessing a pixel classification to separate cells from the background.Then, we performed the object classification task.In this last part, cells were classified manually into three labels: top alive, bottom alive, and top dead.We identified top and bottom alive cells with the help of the fluorescent staining -red and green colors we used previously-; meanwhile, we recognized top dead cells by eye inspection and by comparison with the changes observed in the time evolution of the original files for each well.Finally, the classification was validated using the test dataset; for this, we compared the classifier's cell counts output for the test dataset by manually counting the number of cells in each frame of the dataset.

Supplementary Figure 8: Validation of cell counts derived from the chemotaxis Sartorius IncuCyte assay.
We processed 6000 raw images from the Incucyte transwell assay (80 wells; 25 times points per well; and three channels Phase, Top, and Bottom).To create a training and test dataset, we cropped two random regions of 200×200px from each frame selected.We used the image classification and segmentation software ilastik28 to classify cells as "top alive", "bottom alive", and "top dead".Shown here is the comparison between pipeline's output on the test dataset (y-axis) and manual counts (x-axis) per each of the three classes.

ODE model of cell migration in response to glucose deprivation
In order to fit the data (cell counts) obtained from the image analysis described in section 3.3, we developed an ODE model which describes the temporal evolution of three cell compartments: Top alive (5 0 ), Bottom alive (\ 0 ), and Top dead, (5 " ), as follows: In this model, 5 0 cells can proliferate at a rate ] F , migrate to the bottom chamber at a rate ^ and become \ 0 cells, or die at a rate _ and become 5 " cells.Furthermore, since there are enough resources at the bottom chamber, cells at the bottom (\ 0 ) can only proliferate at a rate ] h and do not die.Additionally, all parameters were considered to depend on the glucose concentration in the top chamber o F , but only the proliferation rate of bottom cells ] h remained constant.The reason for keeping ] h constant is that the conditions at the bottom chamber are always the same among wells (10% FBS and 25 mM of glucose), the only difference is the various glucose concentrations at the top chamber.We fitted the ODE model to the data from the five top chamber glucose concentrations (0.1, 0.5, 1, 5, and 25 mM) simultaneously.To do so, we used the lmfit library in Python and the method of Levenberg-Mawquardt.In total, 16 parameters and 15 initial conditions were allowed to vary during the fitting.

Sequencing
We performed whole genome sequencing (WGS) of four surgical specimens collected during the 1 st and 2 nd surgeries of a primary GBM patient.Sequencing adapters were trimmed from whole genome sequence data using CutAdapt 93 .Trimmed sequences were aligned to the GRCh38 human reference using the Burrows-Wheeler Aligner (bwa mem) 94 .Duplicate identification, insertion/deletion realignment, quality score recalibration, and BAM sorting and indexing were performed with PICARD, the Genome Analysis ToolKit (GATK) 95 , and samtools 96 .Germline variants were identified using GATK UnifiedGenotyper on all tumor and normal samples simultaneously.Allele frequencies for all tumor samples were extracted at normal sample heterozygous variant positions from the GATK output VCF file using custom Perl scripts.
Copy-number based subclonal architecture and how it changes between the two surgeries was derived using HATCHet 59 .Briefly, read depth was counted in 50kb bins.Read depth and B-allele frequency were integrated in the 50kb bins, which were then prepared for HATCHet analysis.The clustering step (cluBB) used tolerancebaf=0.08and tolerancerdr=0.2,and clusters were evaluated visually for appropriate cluster assignment.HATCHet was run with a range of 2 to 8 clones, minimum clone proportion of 0.03, ghostproportion of 0.35, and upper bound on relative increase of objective function of 0.6.

H&E image analysis
d. o A ≔ # $ pppp is the average glucose recorded in the glucose atlas # $ (see also Supplementary section 2.1).
e. Θ A ≔ # # pppp is the average oxygen recorded in the oxygen atlas # # (see also Supplementary section 2.2) .
Next, the T1 scan (and its mask) was carried to MNI space using a linear followed by nonlinear transformation.The T2 flair scan was next coregistered to the T1 scan via an affine rigid body transformation and this transformation matrix was applied to the T2 mask.Both the T2 flair scan and the mask were then registered to MNI space using the same transformation as the T1 scan and mask.The tumors distorted normal cortical anatomy so the skull stripping was somewhat imperfect; a binary mask was generated from the standard MNI brains and applied to the MNI-registered brains with tumors to exclude residual skull for our simulations.To approximate relative cell density (RCD) within the patients' brain from T1 post and T2 flair signals we used a previously described multiparametric regression model for mapping glioblastoma cellularity 97 , giving us four relative cellularity matrices before and after 1 st and 2 nd surgery: RCD kH , RCD k@ , RCD lH , RCD l@ .To go from relative to absolute cell density (ACD), we assumed that the maximum RCD value before the 1 st surgery corresponds to the tumor's carrying capacity (see Table 1), giving us a normalization factor: This factor was used to approximate ACD before and after each surgery as: , where x ∈ {B1, B2, A1, A2}.A drawback of this approach is that T2 flair signal after surgery is not only influced by tumor cell density but also by inflammation caused by the surgery itself.To subtract the effect of inflammation, we enforced tumor cell density after a given surgery not to exceed tumor cell density prior to the surgery: The underlying assumption is that surgery can only decrease, and not increase tumor cell density per voxel.Lastly we set cell density inside the resection cavity to zero: where RÉ$ denotes the resection cavity map.Because stiffness, glucose, and oxygen maps all share the same MNI coordinate system as the ACD maps derived above, we have access to proxies of these environmental factors throughout the tumor and the healthy brain tissue.

Unit Tests
We performed a series of unit tests to verify the implementation of eqs.`1-`11.These tests verify the behavior of birth/death (Supplementary Figure 12A), diffusion (Supplementary Figure 12B), angiogenesis (Supplementary Figure 12C) and metabolic switches (Supplementary Figure 12D) in isolation from all other functionality.The exact conditions underlying each test and the respective expected outcome are summarized in Supplementary Figure 12.
For the diffusion test (Supplementary Figure 12B), we seed the initial cell population in a single voxel, precisely the middle one in position (50,50).We also keep constant population size by considering null birth and death rates.We don't consider angiogenesis and we set a positive migration rate for cells, together with constant resources of glucose and oxygen throughout the brain.We then plot the distribution of cells at different times of the cancer progression, and we compare these results with the following theoretical expectation (solution to diffusion equation in 2D from point source in infinite domani): where: • ∆c and ∆u are the distance between each voxel in the equispaced grid (in our case then they are equal to 1); • ≥ A and c A are respectively the size of population and the position of the voxel cells are initially seeded; • I is the number of days for which we are approximating the distribution; • ∂(c, I) is the approximation of diffusion from equation `6 for timepoint I, calculated as: The simulations match the theoretical expectation for different time frames (Supplementary Figure 12B).
For the angiogenesis tests (Supplementary Figure 12C-E), we set cell migration to zero, and we seed the cells at 1/8 of the carrying capacity uniformly throughout the brain, together with low constant resources of glucose and optimal stiffness.We perform two different analyses for understanding the role of angiogenesis in tumor spreading and glucose concentration over time and space.In Supplementary Figure 12C, we show how the average glucose concentration over time and space is correlated to the parameters D = .In the test in Supplementary Figure 12D, we perform an analysis of the correlation between glucose concentration and angiogenesis activity, showing that the first is much higher where angiogenesis is active.The glucose concentration on the y-axis is calculated by taking the mean over time among the voxels with active angiogenesis, more precisely: where m 0 (I) and m &0 (I) are respectively glucose concentration where angiogenesis is active and the one where angiogenesis is not active at time I.

Data fitting
The S3MB domain was initiated from 4mm 3 voxels of the MNI space, with segmented T1 post and T2 flair regions encoding for different tumor cell densities per voxel.A representative slice was taken to include the lateral ventricles and used as the initial plane for the 2D simulations (Figure 3C), along with its WGS-derived ploidy composition (Figure 3A).We performed a parameter search to recapitulate the GBM's cell density and ploidy composition before the 2 nd surgery.Fitting was done in two steps: (i) modeling just a single subpopulation and using tumor cell density alone for parameter optimization; (ii) modeling both subpopulations but fixing a subset of the parameters informed by step (i) to be identical for both populations, while the remaining parameters were further optimized.For (ii) tumor cell density as well as subpopulation composition were used for parameter optimization.

Modeling both subpopulations
Let N, P ∈ ℝ ) H,@ be the predicted and measured proportions of the two tumor subpopulations respectively (see Figure 3A), also on the day the second surgery.We define: The cost is then calculated as: In a first step, the following four parameters were fitted according to the approach described in 6.1: (i) maximum migration speed (k 30: ); (ii) glucose dependent death constant (ì " in eq.`5 ), (iii) glucose dependent growth constant (ì $ in eq.`5) and (iv) starvation shrinkage rate (_ R ).Then, two of these parameters (ì " and ì $ ) were centered according to the results from the first step and fitted at the resolution of the two aneuploid subpopulations according to 6.2.