Simulation of head and neck cancer oxygenation and doubling time in a 4D cellular model with angiogenesis

Tumor oxygenation has been correlated with treatment outcome for radiotherapy. In this work, the dependence of tumor oxygenation on tumor vascularity and blood oxygenation was determined quantitatively in a 4D stochastic computational model of head and neck squamous cell carcinoma (HNSCC) tumor growth and angiogenesis. Additionally, the impacts of the tumor oxygenation and the cancer stem cell (CSC) symmetric division probability on the tumor volume doubling time and the proportion of CSCs in the tumor were also quantified. Clinically relevant vascularities and blood oxygenations for HNSCC yielded tumor oxygenations in agreement with clinical data for HNSCC. The doubling time varied by a factor of 3 from well oxygenated tumors to the most severely hypoxic tumors of HNSCC. To obtain the doubling times and CSC proportions clinically observed in HNSCC, the model predicts a CSC symmetric division probability of approximately 2% before treatment. To obtain the doubling times clinically observed during treatment when accelerated repopulation is occurring, the model predicts a CSC symmetric division probability of approximately 50%, which also results in CSC proportions of 30–35% during this time.

clinical data for HNSCC. While tumor irradiation was not simulated, tumor growth kinetics during accelerated repopulation were obtained by increasing the CSC symmetric division probability.

Methods
The tumor growth model. Simulations of HNSCC tumor growth were performed using a computational model that was developed in-house using Matlab (version R2017a, The MathWorks, Inc.) and has been previously described 10 . The flow chart in Fig. 1 outlines the spatial and temporal features of the model and how they are related. Briefly, each tumor cell is modeled as an ellipsoid and packed into randomized positions in 3D space without overlap (Fig. 2a). The tumor grows over time by cell division, wherein a cell upon reaching the end of its cell cycle time (CCT) divides into two daughter cells, consequently pushing neighbouring cells outward towards the tumor periphery. A hierarchy of cell types is simulated, including CSCs, three generations of transit cells (T1-3) and differentiated cells (Fig. 2b). The probability for CSCs to undergo symmetric division (i.e., divide into two CSCs as opposed to one CSC and one transit cell) is set by the user. The sloughing of differentiated cells, which is characteristic of epithelial tissue, is also simulated.
Angiogenesis is modeled reflecting a connected and chaotic tumor vasculature that grows out with the cells (Fig. 2c,d), with blood vessels represented by consecutive discrete vessel units. Tumors can be grown with different vascularities which are quantified by the relative vascular volume, RVV. Cellular pO 2 is modeled dynamically as a function of distance from the nearest vessel using a diffusion equation (Table 1), with key parameters being the blood oxygenation, p 0 , and the distance from vessels to the onset of necrosis (the necrosis distance, ND). The tumor vascularity (RVV) and the blood oxygenation (p 0 and ND) affect the amount of hypoxia in the tumor (Fig. 2e,f). Hypoxic cells have longer CCTs and cells that become necrotic are gradually resorbed by the tumor. Table 1 summarises the main parameters of the model and their values for HNSCC.
Simulations entail the following. First, a unique 3D mesh of non-overlapping cell/vessel unit positions is generated using Monte Carlo methods. The cell density reached is 2 × 10 8 cells/cm 3 . A connected network of blood vessels is then generated. As a result, a selection of the mesh positions are designated as vessel unit positions. A unique vasculature is generated each time using Monte Carlo methods. The vasculature is chaotic and tortuous, representative of tumor vasculature 1, 2 . As the tumor grows larger, the vasculature grows out by activating more of the vessel unit positions. Mesh positions that are not designated as vessel unit positions are cell positions, meaning tumor cells may occupy them during tumor growth simulation.
Once a blood vessel network has been generated and prior to tumor growth simulation, the oxygen tension at each cell position is determined and used to calculate the CCT. Cells push one another around when a cell divides, a differentiated cell is lost or a necrotic cell is resorbed. When a cell changes position, it retains its age (the time since it last divided), but its CCT changes to the CCT at its new position. When the age of the cell equals its position dependent CCT, it divides. When a cell divides, it pushes a neighbouring cell outward towards the tumor periphery, causing a chain of cell movement outward, making room for the additional daughter cell. Thus, one daughter cell occupies the position where the parent cell used to be, and the other daughter cell occupies an adjacent position.
The daughter cells are always one generation more differentiated than the parent cell (CSC → T1 → T2 → T3 → differentiated), except in the case of CSC symmetric division. A differentiated cell loss frequency of 80% is simulated to model the natural cell death of these cells, i.e., apoptosis. When a cell becomes differentiated, after a time equal to the CCT, there is an 80% likelihood that the differentiated cell is removed from the tumor. If it is not removed, it remains for another period of time equal to the CCT, then there is again an 80% likelihood that it is removed, and so on. When a differentiated cell is removed, there is a chain of cell movement inward to fill the vacant position. The same occurs when a necrotic cell is resorbed from the tumor, which occurs when its age reaches the necrotic cell resorption time (Table 1).
For a more in-depth description of the computational model methods, please refer to ref. 10. Study of tumor oxygenation. The HNSCC tumor model was used in this work to quantify how the vascularity (RVV) and the blood oxygenation (p 0 and ND) affect the tumor oxygenation. Noting that HNSCC exhibit RVV from 2-10% 11,12 , p 0 from 20-100 mmHg 11,13 and ND from 80-300 µm 14,15 , three combinations of (p 0 , ND) were considered, namely (20 mmHg, 80 µm), (40 mmHg, 180 µm) and (100 mmHg, 300 µm), depicting scenarios of poor, moderate and high blood oxygenation, respectively. In each case, the model input parameter RVV 0 (which would be equal to the tumor RVV if the tumor grew with spherical symmetry, but due to preferential growth along the vessels, RVV ends up larger than RVV 0 ) was varied from 0-10% in 1% increments (for the (20 mmHg, 80 µm) combination, RVV 0 values of 0.25% and 0.5% were also used), yielding values of tumor RVV from 0-16%. The tumor RVV was determined at the end of the growth simulation as the ratio of the number of vessel units to the number of living cells + necrotic cells + vessel units (×100%). Tumor growth simulations began with approximately 70 CSCs and ended with a final tumor diameter of 1 mm (10 4 -10 5 cells). A CSC symmetric division probability of 50% was used in the tumor oxygenation study for fast computations, since this affects the doubling time and CSC proportion but does not greatly affect the tumor oxygenation.
The tumor oxygenation at the end of the growth simulation was evaluated using several different descriptors. The hypoxic proportions HP 10 , HP 5 , HP 2.5 and HP 1 were determined, which were the proportions of living cells with pO 2 <10, 5, 2.5 and 1 mmHg, respectively. The mean and median cellular pO 2 in living cells were also calculated. The volume proportion of necrosis (necrotic volume) was evaluated as the ratio of the number of necrotic cells to the number of living cells + necrotic cells + vessel units (×100%).

Study of tumor growth rate and CSC proportion.
The HNSCC tumor model was then used to explore the effects of tumor oxygenation and CSC symmetric division on the doubling time and the CSC proportion. The doubling time, T d , in the final days (in "tumor time") of the tumor growth simulation was evaluated as follows. Let N(t) denote the number of living and necrotic cells in the tumor at time t. Then the average slope, k, of the curve lnN vs t in the final few days of the simulation was used to calculate the final doubling time according to: The effect of tumor hypoxia and necrosis on doubling time was observed in the simulations from the tumor oxygenation study. Since these simulations all used a CSC symmetric division probability of 50%, the relative variation in the doubling time was reported.
To investigate the effects of the CSC symmetric division probability on the doubling time and the CSC proportion, CSC symmetric division probabilities of 2%, 5%, 10%, 25%, 50%, 75% and 100% were simulated for the two extremes of HNSCC tumor oxygenation. The most oxygenated case was RVV = 10%, p 0 = 100 mmHg and ND = 300 µm and the most hypoxic case was RVV = 2%, p 0 = 20 mmHg and ND = 80 µm. In order to achieve approximately these RVVs, the model parameter RVV 0 was set to 8.2% and 0.75% respectively. Again, the simulations began with approximately 70 CSCs and ended with a final tumor diameter of 1 mm. Three simulations (n = 3) were conducted for each value of CSC symmetric division probability for both well oxygenated and severely hypoxic cases (with the exception of the severely hypoxic case with CSC symmetric division probability 10%, for which n = 4 was used). The CSC proportion was calculated as the ratio of the number of CSCs to the number of living cells (×100%). The doubling times and CSC proportions were plotted using the mean value of the 3 (or 4) simulation runs and with error bars corresponding to the standard error of the mean (SEM). Prism (version 7, GraphPad Software, Inc.) was used to determine whether statistical significance had been reached.

Equipment.
Simulations were performed on the Phoenix cluster at the University of Adelaide 16 using as many as 12 cores and 10 GB of RAM.
Data availability. The data that support the finding of this study are available from the corresponding author upon reasonable request. Code availability. The code used to analyze the data is available from the corresponding author upon reasonable request. The code used to perform tumor growth simulations has not been made publicly available at this time.
The effects of hypoxia and CSC symmetric division on the tumor growth rate and the CSC proportion. The doubling time decreased with increasing tumor vascularity (RVV) and increasing blood oxygenation (p 0 and ND) (Fig. 5a). For HNSCC (RVV = 2-10%, p 0 = 20-100 mmHg, ND = 80-300 µm), the doubling time increased by a factor of 3 from well oxygenated tumors to the most hypoxic tumors. The doubling time was considerably affected, even without the presence of necrosis, by low cellular pO 2 effects such as increased CCTs and cell quiescence (Fig. 5b,c).
The doubling time decreased with increasing probability of CSC symmetric division probability (Fig. 6a). The mean doubling time was 2.6-3.3 times larger for the most hypoxic tumors than for well oxygenated tumors of HNSCC across all values of CSC symmetric division probability. The difference in doubling time between severely hypoxic and well oxygenated conditions was significant (p-value < 0.05 using unpaired t-test with Welch's correction) at every value of CSC symmetric division probability.
The CSC proportion increased with CSC symmetric division probability (Fig. 6b). The mean CSC proportion was 1-1.14 times larger for the most hypoxic tumors than for well oxygenated tumors across all values of CSC symmetric division probability. The difference in CSC proportion between severely hypoxic and well oxygenated conditions was significant (p-value <0.05 using unpaired t-test with Welch's correction) at every value of CSC symmetric division probability except 2%.

Discussion
Several studies have reported pO 2 measurements in human HNSCC using invasive polarographic needle electrodes. The results from some of these studies are collated in Table 2. Since accessing a tissue block from the HNSCC primary tumor can be difficult, these studies often took measurements from sufficiently large lymph node metastases originating from a primary HNSCC. The table also lists measurements of necrotic volume in HNSCCs, as assessed by CT scan or MRI. The results from the current study are included at the bottom for comparison.
The values of median pO 2 , mean pO 2 , HP 2.5 , HP 5 and HP 10 produced by the tumor growth model using RVV from 2-10%, p 0 from 20-100 mmHg and ND from 80-300 µm are in line with these clinical measurements. This assists in model validation since these values for RVV, p 0 and ND are based on clinical data for HNSCC. The clinical tumor oxygenation data overall indicate well oxygenated tumors are rare and lower values of RVV, p 0 and ND are typical. The clinical studies sometimes reported large values of necrotic volume outside the range produced by the tumor growth model. This is likely because the clinical studies observed macroscopic regions of necrosis. Macroscopic necrosis occurs when tumors become large and whole macroscopic regions of the tumor lose blood supply. The tumor growth model only produced necrosis at the microscopic scale (between distant blood vessels) in the sub-clinical sized tumors used in this study.
Pre-treatment doubling times of HNSCC have been obtained in studies that measured tumor growth while patients waited for treatment. Jensen et al. 17 found that in the time between diagnostic scan (MR or CT) and treatment planning CT scan (median 28 days, range 5-95 days), the median doubling time was 99 days (range 15 to > 234 days) for 61 patients with HNSCC. Waaijer et al. 18 found that in the time between diagnostic and treatment planning CT scans (mean 34 days), the mean doubling time was 96 days (range 21-256 days) in 13 patients with oropharyngeal SCC. Murphy et al. 19 found that in the time between diagnostic (MRI or CT) and planning or interval CT scan (median 35 days, range 8-314 days), the median doubling time was 94 days (range 16-6931 days) in 85 oropharyngeal SCC. These average clinically measured pre-treatment doubling times are similar to those produced by the presented tumor growth model under moderately hypoxic conditions and with a CSC symmetric division probability of approximately 2% (recall the doubling times obtained with the tumor model using 2% CSC symmetric division averaged 45 days for well oxygenated tumors and 130 days for severely hypoxic tumors).
In the current work, a CSC symmetric division probability of 2% yielded a proportion of CSCs in the tumor of approximately 6% for all HNSCC tumor oxygenation levels. Methods have been established for identifying CSCs in HNSCC. For example, cells that express markers such as ALDH1, CD133 and CD44 exhibit CSC-like properties, while others do not 9, 20-24 . Cells that efficiently efflux Hoechst 33342 dye, termed side-population (SP) cells, are also CSC-like 25,26 . Chinn et al. 27 reported a mean CD44 high content of 10.8% (range 0-84.5%) for 40 patient-derived primary HNSCCs. In 10 human oral SCC tissue samples, Zhang et al. 28 found CD133 + content of 1-3%. Lu et al. 29 identified approximately 2.1% SP cells in 7 human primary HNSCC samples, and all SP cells were also CD133 + . The CSC proportion of 6% obtained by the tumor growth model using 2% CSC symmetric division, is close to these clinical estimates for HNSCC pre-treatment and results from other models (e.g. 5.9% from Marcu & Marcu 30 ). Tumors respond to treatment by undergoing accelerated repopulation. In an analysis of 5 clinical trials containing a total of 2653 patients, Pedicini et al. 31 using an analytical/graphical method arrived at a best estimate of 3.5 days (95% CI 3.1-3.9 days) for the doubling time of HNSCC during radiotherapy. The loss of asymmetric division by CSCs is believed to be a key mechanism behind accelerated repopulation 6-8, 30, 32, 33 . In the tumor growth model, a CSC symmetric division probability of 50% yielded doubling times from 2.3 to 6.1 days, depending on the tumor oxygenation, which are in line with the estimate by Pedicini et al. A 50% CSC symmetric division probability yielded CSC proportions from 30-35% in the current work. To the authors' knowledge, there are no clinical studies in the literature that measured the CSC proportion in HNSCC in patients during accelerated repopulation. In the model by Marcu & Marcu 30 , the CSC proportions obtained were higher than in the current work for the same CSC symmetric division probability. For example, in their work, 10%, 20% and 30% symmetric division yielded 25%, 35% and 45% CSCs, respectively (recall in the current work, 10% and 25% symmetric division yielded approximately 9% and 15% CSCs, respectively). Conversely, the CSC proportions were slightly lower in the HYP-RT model by Harriss-Phillips et al. 33 than in the current work. In that model, 30% symmetric division yielded just 10% CSCs. Most in vitro studies of various cancer types show a 3-5 times increase in CSCs post single irradiation 24 .

Conclusion
The current work established how the tumor oxygenation varies with vascularity and blood oxygenation, how the doubling time varies with tumor oxygenation and CSC symmetric division probability, and how the CSC proportion varies with CSC symmetric division probability in a 4D cellular model of HNSCC tumor growth. The doubling time varied by a factor of ~3 from well oxygenated tumors to the most severely hypoxic tumors of HNSCC. A CSC symmetric division probability of 2% yielded clinically relevant doubling times and CSC proportions for HNSCC before treatment, while a value of 50% produced the doubling times observed in the clinic for HNSCC undergoing accelerated repopulation. This 50% probability yielded CSC proportions from 30-35%.
In future work, the tumor growth model will be extended to a radiotherapy simulation tool for both low and high LET beams. The cellular geometry will be imported into Geant4 34 and irradiated in Monte Carlo track structure simulations. Radiolysis will be simulated along the particle tracks. Ionisation events and generated • OH species will be clustered in the cell nuclei to predict the complexity and extent of DNA damage to each cell. The cellular pO 2 will affect how efficiently • OH attack to the base of DNA is translated to strand breakage 35 . Irradiation will be simulated in fractions separated by time intervals, during which the tumor growth model will translate DNA damage to cell death while also regrowing the tumor.  Table 2. Clinical data for tissue pO 2 and necrotic volume in human HNSCC. *Mean ± SD of two or more subgroups were combined with appropriate error propagation.