Speed Switch in Glioblastoma Growth Rate due to Enhanced Hypoxia-Induced Migration

We analyze the wave speed of the Proliferation Invasion Hypoxia Necrosis Angiogenesis (PIHNA) model that was previously created and applied to simulate the growth and spread of glioblastoma (GBM), a particularly aggressive primary brain tumor. We extend the PIHNA model by allowing for different hypoxic and normoxic cell migration rates and study the impact of these differences on the wave-speed dynamics. Through this analysis, we find key variables that drive the outward growth of the simulated GBM. We find a minimum tumor wave-speed for the model; this depends on the migration and proliferation rates of the normoxic cells and is achieved under certain conditions on the migration rates of the normoxic and hypoxic cells. If the hypoxic cell migration rate is greater than the normoxic cell migration rate above a threshold, the wave speed increases above the predicted minimum. This increase in wave speed is explored through an eigenvalue and eigenvector analysis of the linearized PIHNA model, which yields an expression for this threshold. The PIHNA model suggests that an inherently faster-diffusing hypoxic cell population can drive the outward growth of a GBM as a whole, and that this effect is more prominent for faster-proliferating tumors that recover relatively slowly from a hypoxic phenotype. The findings presented here act as a first step in enabling patient-specific calibration of the PIHNA model.


Introduction
Glioblastoma (GBM) is the highest grade of glioma from the World Health Organization (Louis et al. 2016). It is uniformly fatal with an average survival time from diagnosis of only 15 months with standard of care treatment (Stupp et al. 2009). The standard therapy regime for this disease is a combination of resection, radiation and chemotherapy (Stupp et al. 2005(Stupp et al. , 2009. Magnetic resonance imaging (MRI) is the standard imaging modality for GBMs and is used routinely to monitor tumor growth and development throughout the progression of the disease. Different MRI sequences such as gadolinium-enhanced T1-weighted (T1Gd) and T2-weighted (T2) are used to identify the gross tumor volume. T1Gd shows gadolinium that has leaked into brain tissue, and T2 shows water that has done the same, which is known as edema. However, these MRI sequences together do not show a complete picture. Infiltrating tumor cells also exist beyond the resolution of these MRI sequences. In fact, malignant glioma cells have been cultured from histologically normal healthy tissue at a distance of 4cm from the gross tumor volume identified by MRI scans (Silbergeld and Chicoine 1997).
Hypoxia has been shown to induce more migration in glioma cells (Keunen et al. 2011;Zagzag et al. 2006). There is also evidence that glioma cells follow a dichotomy of migration and proliferation (Giese et al. 2003) and evidence of a lower proliferation marker for cells that exist in hypoxic regions of GBMs (Brat et al. 2004). Tumors in hypoxic conditions release angiogenesis-promoting factors to encourage vessels to grow toward them and provide nutrients (Gordan and Simon 2007;Korkolopoulou et al. 2004;Zagzag et al. 2000). This process also occurs in normoxic conditions at a lower level (Zagzag et al. 2000). Necrosis occurs in the vast majority of GBMs and presents in the core of the tumor (Louis et al. 2016). Necrotic cells can lead to an unfavorable local microenvironment that injures nearby cells and subsequently spreads cell death (Raza et al. 2002;Yang et al. 2013).
Over the past 20 years, there have been many partial differential equation models that simulate GBM cell density and have provided various insights into this disease (Hawkins-Daarud et al. 2013;Martínez-González et al. 2012;Massey et al. 2012;Swan et al. 2018;Swanson 1999;Swanson et al. 2000Swanson et al. , 2003aSwanson et al. , 2008Swanson et al. , 2011. One such model is the Proliferation Invasion Hypoxia Necrosis Angiogenesis (PIHNA) model, which has been used to analyze the mechanistic properties of GBMs that lead to observed imaging features and has shown similar growth and progression patterns to those seen in patient tumors (Swanson et al. 2011). We carry out a traveling wave analysis on the PIHNA model to determine which parameters drive the outward growth of the tumor as a whole and compare these analytical predictions with computational simulations in the cases of varying relative rates of migration between hypoxic and normoxic tumor cells. We find that the traveling wave dynamics only depend on the equations for the normoxic and hypoxic tumor cell densities. We find that the normoxic cell migration and proliferation rates, D c and ρ, respectively, drive the minimum wave-speed in the PIHNA model, which is given by and also depends on the initial background vasculature in the model, v 0 , relative to the spatial carrying capacity, K . We find that s min holds for published results using the PIHNA model as they have not allowed for different hypoxic cell and normoxic cell migration rates. We allow these migration rates to be different in the model and observe the effect of this variability on simulated tumor growth rates. We find that a fasterthan-minimum wave-speed is achieved when hypoxic cells migrate sufficiently faster than normoxic cells and find a threshold above which these dynamics can occur. This threshold depends on the proliferation rate ρ, the switching rate back from hypoxic cells to normoxic cells γ , and v 0 /K . We denote this threshold k, and it is given by These results are then confirmed and explored computationally through further model simulations. The PIHNA model therefore suggests that hypoxic cell migration, if sufficiently fast, is able to drive the outward growth of the tumor as a whole.
The PIHNA model has been used in various settings to explore possible mechanistic explanations for clinical observations of GBM (Hawkins-Daarud et al. 2013Swanson et al. 2011). It was built out of a much simpler model, the Proliferation Invasion (PI) model, a basic diffusion/logistic growth model whose simplicity has allowed for patient-specific calibration (Swanson 1999;Swanson et al. 2000Swanson et al. , 2003b. Despite its simplicity, the patient-specific calibration of the PI model has proven clinically prognostic for many aspects of clinical care (Adair et al. 2014;Baldock et al. 2014a, b;Harpold et al. 2007;Neal et al. 2013a, b;Singleton et al. 2019;Swanson et al. 2008). The increase in variables and parameters of the PIHNA model enables it to capture a wider range of clinical scenarios and questions; however, it has limited the ability for patient-specific calibration. The results presented here represent a first step in enabling patient-specific calibration of the PIHNA model as it shows the mathematical relationship between the wave speed and a handful of critical parameters, the same key relationship used to calibrate the PI model. Further, this relationship sheds light on expected tumor behavior based on the degree of aggressive hypoxia which is imageable through positron emission tomography (PET) scans, which, with further study, may influence clinical decision making.
We introduce the PIHNA model in the next section before calculating the expression for the minimum wave-speed in Sect. 3. In this section, we also find the threshold, k, on the relative migration between hypoxic and normoxic cells under which the minimum wave-speed is achieved. We then move onto PIHNA simulations in Sect. 4 to computationally validate our findings.

The PIHNA Model
The PIHNA model (Swanson et al. 2011) simulates five different species and their interactions: c-the density of normoxic tumor cells, h-the density of hypoxic tumor cells, n-the density of necrotic cells, v-the density of vascular endothelial cells, a-the concentration of angiogenic factors.
The dimensions of c, h, v and n are cells/mm 3 of tissue. The angiogenic factor, a, is a diffusing concentration with dimensions μmol/mm 3 tissue.
Normoxic cells proliferate with rate ρ and migrate with rate D c , whereas hypoxic cells do not proliferate but migrate with rate D h . Cells convert between normoxic and hypoxic phenotypes depending on the ability of the local vascular density to provide nutrients at their location; hypoxic cells in the model become necrotic if they remain in such a region. When any other cell type meets a necrotic cell, they become necrotic with rate α n . Previous publications on the PIHNA model have set the migration rate of hypoxic cells to be equal to that of normoxic cells, such that D h = D c . However, hypoxia has been shown to promote GBM cell migration, so we have allowed for this to be varied in the PIHNA model (Keunen et al. 2011;Zagzag et al. 2006). and The term V models the relationship between the vasculature and its effect on the tumor. Note that V take values in [0, 1] such that it affects the switching rates between the populations c, h and n. A value of V (c, h, v) ≈ 0 corresponds to a very inefficient vasculature that cannot provide sufficient nutrients to the local tumor population; this would increase the conversion of normoxic cells to hypoxic cells and in turn necrotic cells. A high V (c, h, v) ≈ 1 promotes a normoxic phenotype. It is worth noting that once necrotic cells are present in a simulation, they will always increase in population due to the contact necrosis in the model, which represents the injury of nearby cells and promotion of their necrosis. Further definitions can be found in Table 1.
The expression for T defined in Eq. (5) is a spatiotemporal measure of the relative density of the cells in a region. It is used to limit growth and migration and used as a threshold to determine which densities would appear on different MRI sequences. Substituting the set of Eqs. (3) into Eq. (5) gives from which it is clear that at T = 1 the reaction and diffusion terms vanish, which implies T is restricted by the upper bound of 1 (as long as T (x, 0) ≤ 1). As T is a sum of nonnegative components and K > 0, we have that T ≥ 0. Therefore, we have that T ∈ [0, 1] for sufficient initial conditions, for all x and t ≥ 0. Following the literature, we have assumed that a total relative cell density of at least 80% is visible on a T1Gd MRI, and a total relative density of at least 16% is visible on a T2 MRI (Swanson 1999;Swanson et al. 2011). In the PIHNA model, this translates to T ≥ 0.8 being visible on T1Gd MRI and T ≥ 0.16 being visible on T2 MRI. By construction, the T1Gd radius is always less than or equal to the T2 radius, which agrees with patient data (Harpold et al. 2007).
For the purposes of the wave-speed calculations, we consider the PIHNA model in a one-dimensional spherically symmetric case with zero-flux boundary conditions at the end points, r = 0 and r = r end . This does not take into account the full anatomy of the brain, but it is useful to gain insight into the behavior of the PIHNA model. The initial condition is given by A justification of parameters can be found in the supplementary material of Swanson et al. (2011). * We have altered these rates in this formulation of PIHNA, which have not been changed previously to simulate a small initiating population of normoxic tumor cells decreasing away from the core of the tumor. We also have h(r , 0) = 0, n(r , 0) = 0, v(r , 0) = 0.03K and a(r , 0) = 0. We run the PIHNA simulations with the parameters found in Table 1. In all simulations, the tumor and necrotic cell densities spread outwards. A peak in normoxic cell density leads and is followed by a peak in hypoxic cell density and then a zone of necrosis develops in the core of the tumor, as can be seen in Fig. 1; this figure also shows how we calculate the wave-speed values from simulations. c The T2 radius shown over time for the same simulation. This radial growth is nonlinear for small tumor sizes, but settles to a linear rate, which is the wave speed of the simulation we carried out a wave-speed analysis to find an analytical expression for the tumor wave-speed in the PIHNA model. This wave speed is what has enabled patient-specific calibration of the PI model for GBM patients, and we expect that this similar analysis will eventually allow for the patient-specific calibration of the PIHNA model as well. Note that in spherically symmetric coordinates, the wave speed asymptotically approaches that of a planar wave-speed.
(3)] and discarding nonlinear terms leads to the following set of equations: where we have used T = v 0 /K and V = 1. The equations forĉ andĥ decouple from the system, and it is these two equations that dictate the outward growth rate of the tumor. We will analyze these two equations to look for traveling wave solutions of the form where s is the wave speed. Substituting Eq. (13) into Eqs. (8)- (9), gives rise to the following equations Looking for non-trivial solutions to this system of equations yields eigenvalues as functions of the wave speed, s. We have four eigenvalues, given by We have also found the corresponding eigenvectors for all of our eigenvalues, which we shall denote V i for each λ i . These are given by the following expressions, up to a proportional constant: The terms λ 1,2 are both negative as s > 0 by assumption. Due to positive restrictions on the state space (negative populations do not make any biological sense), a spiral approach around the point (0, 0, 0, v 0 , 0) cannot occur. Therefore, we need the discriminant of the set of quadratic λ 1,2 solutions to be greater than or equal to zero. In other words, Therefore, we have a minimum wave-speed of There is no minimum wave-speed associated with the eigenvalues λ 3,4 .
The PIHNA model will follow this minimum wave-speed if the eigenvalue λ 1 evaluated at this minimum gives the smallest possible negative eigenvalue of λ 1,2,3,4 . If there exists some s > s min such that 0 > λ i (s) > λ 1 (s min ) for some i = 1, 2, 3, 4, we will see the emergence of a solution with a larger wave-speed. In this section, we will compute a threshold below which the minimum wave-speed is achieved but above which other dynamics may emerge.
We will call each eigenvalue evaluated at s min , λ min i , for i = 1, . . . , 4. We start by noting that λ 2 ≤ λ 1 and λ 3 > 0, so neither of those can be negative with a smaller magnitude than λ min 1 to change the PIHNA wave-speed dynamics. As λ 4 becomes less negative for increasing values of D h , there is a threshold value of D h /D c that leads to λ min 1 = λ 1 = λ min 4 for which the minimum wave-speed is still achieved. For values of D h /D c that are smaller than this threshold, the minimum wave-speed will still be achieved. However, larger values of D h /D c may lead to a faster wave-speed, as the eigenvalues become smaller than λ min 1 . We have Setting λ min 1 = λ min 4 and using the expression for s min [Eq. (18)] leads to Solving for D h /D c gives the non-trivial solution We will define this threshold of D h /D c as k. Note that as v 0 ≤ K , k ≥ 2. So for the faster wave-speeds to occur, the hypoxic cell migration rate needs to be at least twice as fast as the normoxic cell migration rate. For D h /D c = 1, as is the case in previous PIHNA publications, we do not expect faster wave-speeds to occur, regardless of other simulation parameters.

Simulation Results
To calculate the simulated wave-speed in numerical simulations, we thresholded the total cell density at T = 0.16, which is a commonly assumed cell density threshold for visible tumor-related abnormalities on T2 MRI (Swanson 1999). Following the establishment of a wave front, the simulated wave-speed levels out to a fixed value, see Fig. 1. We analyze the wave speed of large tumors to ensure we are analyzing established wavefronts while minimizing numerical error. We are particularly interested in the effect on the wave speed of varying hypoxic cell migration rates, more specifically the change in their migration rate compared with normoxic cells (D h /D c ), which has been allowed to vary in the PIHNA model for the first time. We also want to observe the effect of the D h /D c threshold, k, that allows for faster wave-speeds. Although we have already observed that the equations for the normoxic and hypoxic cell densities decouple in the linearized form of the PIHNA model, simulations presented here are for the full PIHNA model. Numerical simulations are run on a spherically symmetric domain, with a step size of 0.01mm. All simulations were run in Matlab 2018a using the inbuilt solver pdepe.m.

Relatively Fast Hypoxic Cell Diffusion Rates Increase Wave-speed
The wave speed for PIHNA simulations with D h /D c ≤ k converges toward s min . However, if we compute the wave speed for simulations where D h > k D c , we see that the wave speed can be faster and continues to increase for larger D h /D c values; an example of this can be seen in Fig. 2. Computing the corresponding eigenvalues shows a change in behavior for values of D h /D c > k. We also plot k on Fig. 2, in which case k = 2.60 (three significant figures). These values of D c and ρ are biologically realistic and based on the mean of previous migration and proliferation rate estimates from the PI model applied to patient-specific MRI data (Wang et al. 2009). From these observations and our analysis in Sect. 3, we can deduce that if the hypoxic cell migration is sufficiently faster than the normoxic cell migration (such that Focusing on the eigenvectors corresponding to the least negative eigenvalues, V 1 and V 4 , we see that they influence the dynamics of the model. By plotting the normoxic cell density across space against the hypoxic cell density across space for a fixed time point where each simulation has converged to a stable wave-speed, together with V 1 and V 4 , we can see how the traveling wave trajectory approaches the state ahead of the wave front. We present two simulations with their corresponding V 1 and V 4 eigenvectors in Fig. 3, one for D h /D c = 10 −1 and another for D h /D c = 10 1 . For D h /D c = 10 −1 , where the wave speed follows the predicted minimum value, we see that the model approaches along (c, h) = (0, 0) along the eigenvector V 1 , whereas for D h /D c = 10 1 , the approach is along V 4 . In the linearized regime, we expect that c =c exp(λ(r − st)), such that To provide further evidence concerning the traveling wave trajectory, we compared the gradient of the log of normoxic cells (c) with the eigenvalues λ 1 and λ 4 . We see Fig. 4 As D h /D c is increased for varying γ values, we see an increase in the converged numerical wavespeed that is more pronounced for smaller values of γ ; the corresponding thresholds k for wave speeds faster than s min are indicated. The values used are γ = 0.005, 0.05 and 0.5/day, with corresponding k values of 2.06, 2.60 and 7.95, respectively. Wave speeds taken as the average speed between 8 and 8.5cm of growth on simulated T2 MRI (16% total cell density threshold) and presented relative to s min that for low values of D h /D c , the gradient more closely follows λ 1 and for large values of D h /D c , the gradient closely follows λ 4 . We present examples of these results in "Appendix" (Fig. 6).

Low Switching Rate from Hypoxia to Normoxia, , Amplifies Wave-speed Increase for Large Values of D h /D c
The switching rate from a normoxic cell to a hypoxic cell (β) is not present in the eigenvalues that dominate the behavior of the wave speed, nor in the expression for k. We do, however, note that the switching rate from a hypoxic cell phenotype back to a normoxic cell phenotype, γ , is present in the expression for λ 4 [Eq. (19)] and subsequently in the expression for the D h /D c threshold, k [Eq. (22)]. We ran a similar set of simulations as in Sect. 4.1 with a higher value of γ = 0.5/ day and a lower value of γ = 0.005/day to verify that the wave speed increase, relative to s min , would be affected for varying γ . As expected, higher values of γ increase k and correspond to a lower wave-speed for equivalent D h /D c values. We present these wave speed results in Fig. 4 where we also mark the corresponding values of k. For γ = 0.005, 0.05, and 0.5/day, we find k = 2.06, 2.60 and 7.95, respectively.

Wave-speed Increase is More Pronounced for Faster-Proliferating Tumors
We also varied ρ to explore its effects on the increase in wave speed for large values of D h /D c . We chose two more values of ρ = 10 1.25 /year (lower ρ), and ρ = 10 2 /year (higher ρ) and refer to the previous simulations with ρ = 10 1.5 /year as a mid-range ρ. Throughout all simulations, we set D c = 10 1.5 mm 2 /year and γ = 0.05/day, leading Eq. (22) to give threshold values of k = 2.19, 2.60 and 3.06 for higher, medium and lower ρ simulations, respectively. We present the wave speeds normalized against their predicted values of s min [Eq. (18)] in order to compare the simulation results across different values of ρ. We see for values of D h /D c below their respective thresholds that the wave-speeds all follow

Discussion
We have found an expression for the minimum wave-speed for the PIHNA model given by and shown that this predicted wave-speed is attained when normoxic cell diffusion is greater than or equal to the diffusion of hypoxic cells. We therefore have shown that the predicted minimum wave-speed is valid for previous publications of the PIHNA model (Swanson et al. 2011). However, due to the in vitro evidence indicating that hypoxia can increase migration (Keunen et al. 2011;Zagzag et al. 2006), we are interested in increasing the migration rate of hypoxic cells compared with normoxic cells in our extension of the PIHNA model. In the case that the hypoxic cells diffuse sufficiently faster than the normoxic cells, we see that the outward growth of the tumor is faster than the predicted minimum wave-speed value. In fact, we have quantified the value at which these faster rates of growth can occur through the threshold and note that the hypoxic cell diffusion has to be at least twice as fast as the normoxic cell diffusion. The threshold of hypoxic to normoxic cell migration rates is increased if the hypoxic cells can easily convert back to normoxia and decreased for faster-proliferating normoxic cell populations. This result suggests that fasterproliferating tumors that can only slowly recover from hypoxia are pushed to grow even faster by a highly migratory hypoxic subpopulation, more so than slower-proliferating tumors that can easily recover from hypoxia. The γ parameter encoding this recovery from hypoxia is the inherent ability of hypoxic tumor cells to adapt to a nutrientrich microenvironment, switch off any hypoxia-related processes and reinitiate those related to a normoxic cell phenotype. This change in microenvironment is represented through the spatial variation in vasculature as the simulated tumors grow. As hypoxic cells migrate (in some cases faster than normoxic cells), they reach regions of more abundant vasculature and convert back to a normoxic phenotype. This recovery from hypoxia would likely also be influential in model dynamics if an initial condition of varying vasculature or treatment effects such as ischemia were introduced into the PIHNA model. The behavior of the PIHNA model suggests that limiting the lasting impact of hypoxia on phenotypic expression may slow the outward growth of GBM as would decreasing the motility of hypoxic tumor cells. Similarly, decreasing the motility and proliferation rates of normoxic cells would decrease the minimum wavespeed, the latter of which is already a widely exploited treatment mechanism through chemotherapy and radiation. Of course, the spherical symmetry assumed in the model does not fully capture the complex heterogeneity present in GBM. Inherent differences in vasculature and nutrient abundance are present in the healthy human brain (Yamaguchi et al. 1986), and spatial heterogeneity in hypoxia is observed in GBM (Bell et al. 2015;Spence et al. 2008). Genetic differences within GBM could drive heterogeneity in all of the tumor-related parameters influencing the PIHNA wave-speed. This environmental and genetic heterogeneity could lead to varying wave-speeds of GBM growth within individual tumors. Even so, we anticipate eventual patient-specific calibration of this model to provide better clinical insights into individual tumor behaviors beyond what the PI model can do. As quantification of aggressive hypoxic volumes becomes more readily available through PET scans, we anticipate this wavespeed estimate to play a critical role in estimating patient-specific parameters for this model.
The analysis presented here shows that the wave-speed dynamics do not depend on the vascular efficiency term, V , as long as V = 1 in its linearized form. We also do not see a dependence on the switching rate from the normoxic cell density to the hypoxic cell density, β. All of the results presented here are dictated by the equations for normoxic cell density and hypoxic cell density due to their independence from the other three equations in the linearized model. Necrosis, angiogenesis and vascular growth dynamics play no role in the outward growth rate in the PIHNA model. We have shown that this is the case both through theory and model simulations. This concept would hold for similar tumor growth models that decouple in their linearized forms and have different motilities between hypoxic and normoxic cell phenotypes.
Mathematically, the increase in simulated wave-speed corresponds to a change in the asymptotic traveling wave trajectory as D h /D c is increased, which causes an eigenvector associated with the hypoxic cell density characteristics to dominate the behavior of the PIHNA model. Biologically, this suggests that the faster migration of hypoxic cells can drive the growth of the whole tumor, as they migrate toward nutrientrich environments and convert back to normoxic cells. If this conversion rate is high, the model suggests that the outward growth rate of the whole tumor is lower. The model does not predict that the wave speed is affected by the proportions of hypoxic and normoxic cells. However, a reduction in vasculature ahead of the wave (v 0 ) does increase the invasion speed of the tumor due to the appearance of v 0 in the minimum predicted wave-speed. It would be interesting in future work to see how including a normal cell density affects these dynamics. Fig. 6 The gradient of the log of the normoxic cells is plotted for a T2 radius of 30 cm. As described in the main text, the leading edge of this simulated gradient (ignoring boundary effects present close to the edge of the domain) should follow the eigenvalue that controls the dynamics of the PIHNA model. Simulations presented here correspond with those presented in Fig. 3. The simulation with D h /D c = 0.1 agrees more closely with λ 1 , whereas the simulation with D h /D c = 10 follows λ 4 . The results of these support the eigenvalue and eigenvector analysis in the main body of this work