The Impact of Phenotypic Switching on Glioblastoma Growth and Invasion

The brain tumour glioblastoma is characterised by diffuse and infiltrative growth into surrounding brain tissue. At the macroscopic level, the progression speed of a glioblastoma tumour is determined by two key factors: the cell proliferation rate and the cell migration speed. At the microscopic level, however, proliferation and migration appear to be mutually exclusive phenotypes, as indicated by recent in vivo imaging data. Here, we develop a mathematical model to analyse how the phenotypic switching between proliferative and migratory states of individual cells affects the macroscopic growth of the tumour. For this, we propose an individual-based stochastic model in which glioblastoma cells are either in a proliferative state, where they are stationary and divide, or in motile state in which they are subject to random motion. From the model we derive a continuum approximation in the form of two coupled reaction-diffusion equations, which exhibit travelling wave solutions whose speed of invasion depends on the model parameters. We propose a simple analytical method to predict progression rate from the cell-specific parameters and demonstrate that optimal glioblastoma growth depends on a non-trivial trade-off between the phenotypic switching rates. By linking cellular properties to an in vivo outcome, the model should be applicable to designing relevant cell screens for glioblastoma and cytometry-based patient prognostics.


Introduction
Cancer progression is the macroscopic outcome of numerous cellular processes, such as elevated proliferation rates, defects in apoptosis regulation and abnormal angiogenesis [1]. In the development of targeted anticancer therapies, the proliferation, survival and angiogenesis phenotypes are often singled out as the most important. Recently, however, much attention has been given to cancer cell migration as a possible therapeutic target, since it underlies both the local invasive process whereby cancer cells degrade and move through the adjacent tissue, and the formation of distant metastases.
The importance of cancer cell migration is perhaps most evident in the common brain tumour glioblastoma, which is characterised by rapid and infiltrative growth into the surrounding brain tissue. In glioblastomas, neoplastic cells are often found at a long distance (several centimeters) from the main tumour mass. This diffuse growth pattern presents a difficult clinical problem, since residual 'satellite cells' can mediate rapid recurrence of the disease after surgery [2]. Key factors that underlie glioblastoma cell invasiveness include high migration speeds in comparison to other types of cancer (up to 100 mm/h) and the fact that brain parenchyma provides a penetrable substrate for invasion [3]. Thus, inhibition of migration pathways might constitute an interesting complement to standard glioblastoma therapies that seek to inhibit cell proliferation rate or angiogenesis. Several pathways have been suggested to mediate the highly migratory phenotype of glioblastoma cells, including signaling via Focal adhesion kinase (FAK) [4], Phosphoinositide 3-kinase PI3K [5] and Signal transducer and activator of transcription 3 (STAT3) [6]. Other concepts for targeting of migration have also been proposed, including inhibition of integrins [7], perturbing the interactions between ECM components [8], and administering lithium chloride [9]. Further, potential gene targets have also been revealed using molecular profiling efforts [10].
However, the potential of migration as a therapeutic target is complicated by the strong dependency between migration and proliferation phenotypes. Early in vitro experiments by Giese et al. [11] showed that when plated on a substrate, that supports migration, the proliferation rate of glioblastoma cells is markedly reduced. Later, it was shown that cells at the tumour's invasive rim proliferate more slowly than cells in the central parts of the tumour, again suggesting that migration has a 'cost' in terms of reduced proliferation [2]. These and several additional observations led to the so called 'go or grow' hypothesis, stating that migration and proliferation are mutually exclusive phenotypes. Although a molecular explanation for this dichotomy still is missing, it has been suggested that cytoskeleton dynamics could be limiting, as it is involved in both cell division and force generation in migration [12].
In recent experiments, in vivo imaging of fluorescent glioblastoma cells enabled direct observation of phenotypic switching between the 'go' state (migration) and the 'grow' state (proliferation). More precisely it was observed that glioma cells move in a saltatory fashion, where bursts of movement are interspersed by periods of immobility, and it is during these stationary periods that the cells divide [13,14]. Taking these observations into account allows for a more comprehensive understanding of glioblastoma progression, where tissue-level traits, such as progression rate, emerge from cell-level behaviour. Mathematical models at the resolution of individual cells enable a quantitative connection between these scales, and can hence be of great assistance.
Here, we focus on the relationship between cell-level phenotypic switching in glioblastoma, and the properties of the tumour as a whole. In particular we elucidate how the growth rate of the tumour and speed of invasion depends on the specific underlying microscopic parameters, such as phenotypic switching rates, rate of apoptosis et cetera. Please note that the we use the word 'invasion' to denote the process by which glioma cells spread into and displace the surrounding brain tissue, and do not refer to branched finger-like growth patterns. Although several models of glioma growth have previously been proposed (see next section), this model is the first to connect experimentally measurable cell-level traits with gross tumour volume in an analytical way. This yields hope for the future understanding of glioma biology and therapy, since it is the understanding of how drug induced changes on the cell-level scale propagate to the organ scale, that are required in order to accurately predict therapeutic outcome.
In the following we first review previous work in the field of glioblastoma modelling and then proceed by introducing our individual-based (IB) stochastic model of glioma growth. From this model we derive an approximate continuum description of the system, whose properties are compared to the IB-model. We proceed to analyse the continuum model to reveal the influence of the model parameters on the rate of spread of the tumour, and finally discuss our results in the context of other models and experimental results.

Previous work
The growth of glioblastomas was first modelled by means of a continuum approach, which captures the two main features of glioma cells: proliferation and migration ( [15,16], and [17], chapter 11). In that model, the partial differential equation (PDE) that describes the time evolution of the concentration of glioma cells u(x,t) in space and time has the form: where the migration is captured by a diffusion term with diffusion coefficient D (first term) and proliferation of the glioma cells is described by a locally logistic growth function with growth rate r (second term). This equation is known as the Fisher (or Kolomogorov) equation, and was first derived in order to describe the spread of an advantageous gene in a spatially extended population [18,19]. The derivation, originally by Fisher and later refined by Kolmogorov, starts by assuming a contact distribution that describes the probability of migration between two spatial locations, and by then assuming that all moments of that distribution higher than two are negligible (known as the diffusion approximation, see for example [20]) one arrives at the above equation. The Fisher equation has been of particular interest since it gives rise to traveling wave solution u(x,t)~U(x{ct), whose shape is preserved and position in space is shifted at a speed c as time progresses. Significant interest has been devoted to determining the wave speed c, and it has been shown that for reasonable initial conditions (exponentially decaying [20] or of compact support [19]) the wave speed is given by c~2 ffiffiffiffiffiffi ffi Dr p [21]. The speed of propagation thus depends on both the motility, captured in the diffusion constant D, and the rate of proliferation r. Both D and r can be determined from time-course Magnetic Resonance Imaging (MRI) data from actual patients, and it has been shown that their values are of prognostic power [22].
The above modeling approach rests on the assumption that glioblastoma cells follow a random walk (which at the macroscopic scale corresponds to the diffusion of cells). Recently this assumption has been under scrutiny, and this has led to a number of explorations of non-random migration, i.e. where migration is influenced by biological processes such as cell-cell signaling, oxygen pressure, nutrient availability and phenotype switching. In one line of work, Aubert et al. [23] used an individual based (IB) model to show that attraction between glioblastoma cells is likely to influence the dynamics of tumour invasion. Deroulers et al. [24], derived the macroscopic PDE for this case, obtaining a density dependent diffusion equation (D~D(u) in terms of eq. (1)), whose solution deviates significantly from the Fisher-Kolmogorov PDE (see also [25,26]). Khain et al. [27], used IB models to characterise the role of hypoxia in glioblastoma, showing that reduced oxygen levels may down-regulate cell-cell adhesion, leading to increased motility.
The cellular behaviour implied in the 'go-or-grow' hypothesis (see Introduction) is also thought to affect migration and growth dynamics of glioblastomas, in a manner that is not captured by the Fisher-Kolmogorov equation. Hatzikirou et al. [28] proposed a lattice-gas cellular automaton model in which the switching between the proliferative (P) state and migratory (M) state is driven by lack of oxygen, and went on to show that in the corresponding macroscopic (Fisher) equation, there is a tradeoff between diffusion and proliferation reflecting the inability of cells to migrate and proliferate simultaneously. Similar results where obtained by Fedotov and Iomin [29] but with a different type of model known as continuous time random walk model. That model contains two distinct subpopulations (P-cells which are stationary and divide and M-cells that perform random walks), and a cell switches from one compartment to the other after a time t p (and t m respectively) which is exponentially distributed. They analytically show that the spreading rate is smaller than one would expect from the Fisher equation (1). Finally, Lewis and Schmitz have studied the general

Author Summary
In this work, we develop a spatial mathematical model in order to analyse the growth behavior of the brain tumour glioblastoma. Tumours of this type have a diffuse boundary, with considerable local invasion of surrounding brain tissue, making surgery difficult. At the cellular level, the progression of a glioblastoma is known to depend on the balance between cell division (proliferation) and cell movement (migration). Based on recent evidence, our model assumes that each cell in a glioblastoma tumour resides in either of two mutually exclusive states: proliferating or migrating. From a probabilistic model of switching between these two phenotypes, we go on to derive equations that link cellular phenotypes to disease progression. The model has several possible applications. For instance, it could be used to predict the rate of disease progression in an individual patient, and to improve screening methods. relation between organism migration and proliferation and its impact on population spread using reaction-diffusion equations [30]. They show that the system exhibits travelling wave solutions and that the wave speed depends on the rates of switching between the states.
The model that we propose draws from these previous models, but is different in some crucial ways. We consider two distinct subpopulations with a stochastic switching in between (as in Fedotov and Iomin, and Lewis and Schmitz), but instead of starting with a continuum description, we begin with an IB-model in which the cells occupy a lattice and obey size exclusion (as in Deroulers et al.), and from that derive a system of PDEs. This allows for an analytical treatment of the IB-model which establishes a connection between cell characteristics and the macroscopic behaviour of the system previously not demonstrated.

Results
The individual-based stochastic model The cells are assumed to occupy a d-dimensional lattice (we will consider d~1,2,3), containing N d lattice sites, where N is the linear size of the lattice and each lattice site either is empty or holds a single glioma cell. This means that we disregard the effects of the surrounding brain tissue, such as the different properties of grey vs. white matter [31], and the presence of capillaries which might influence the behaviour of the cancer cells. But since the soft tissue in the brain presents little resistance to invading cancer cells and the precise nature of interaction with stromal cells is unclear, focusing on the dynamics of the glioma cells is a reasonable first approximation. Further, the process of angiogenesis, which has been modelled extensively [32], is ignored, and we hence assume the growing tumour to be well vascularised. For the sake of simplicity we do not consider any interactions between the cancer cells (adhesion or repulsion), although this could easily be included, and we also disregard other types of mechanical interactions, such as cell pushing (see Discussion).
The lack of knowledge of the intra-cellular dynamics and extracellular cues that lead to the phenotypic switching behaviour poses a problem, but we will circumvent it by, as a first approximation, considering the switching as a stochastic event. The behaviour of each cell is therefore modelled as a time continuous Markov process where each transition or action occurs with a certain rate, which only depends on the current and not previous states, known as the Markov property. The rates are interpreted in the standard way, i.e. if transition in a variable X (t) from state i to j occurs at rate q ij then the probability of a transition occurring in the time interval ½t,tzDt is given by where o(Dt) means that remaining terms are bounded from above by Dt, and thus that in the limit of small Dt the transition probability is proportional to q ij . Each cell is assumed to be in either of two states: proliferating or migrating, and switching between the states occurs at rates q p (into the P-state) and q m (into the M-state). A proliferating cell is stationary, passes through the cell cycle, and thus divides at a rate a. The daughter cell is placed with uniform probability in one of the empty 2d neighbouring lattice sites (using a von Neumann neighbourhood). If the cell has no empty neighbours cell division fails. A migrating cell performs a size exclusion random walk, where each jump occurs with rate n. Size exclusion means that the cell can only move into lattice sites which were previously empty. If only a motile subpopulation is considered, size exclusion does not affect the macroscopic diffusive nature of a population of random walkers, but if two or more subpopulations are taken into account then, as we shall later see, diffusion becomes non-linear. Lastly, cells are assumed to die, through apoptosis, at a rate m independent of the cell state. Since this type of cell death is associated with cell shrinking and rapid removal of the dead cell, a cell which goes through apoptosis is instantly removed from the lattice and leaves an empty lattice site behind.
The stochastic process is depicted schematically in figure 1. In fact the whole system comprises a continuous time Markov chain with a finite, but very large state space, containing N 3d different states, where N is the linear size of the system and the 3 comes from the three distinct lattice states: empty, P-cell and M-cell.

Parameters
We will consider a lattice of linear size N~200 with a spacing of Dx~20 mm, the typical size of a cancer cell. For the most part we will consider the system in d~2 dimensions, which means that we simulate a lattice, which corresponds to a 4|4mm 2 slice of tissue. This is of course considerably smaller than a clinically relevant glioma, but sufficient to capture the effects of the phenotypic switching on tumour growth rate. The time scale of the model is set to agree with that of the cell cycle (approximately 24 hrs [2]) which means that the proliferation rate a~1 cell cycle {1 , and that we scale all other parameters accordingly. We are mainly interested in the effect of the phenotypic switching rates q p and q m on the growth of the tumour and they will therefore be varied within a biologically reasonable range. It follows from equation (2) that the time spent in one phenotypic state is exponentially distributed with parameter q p,m and thus that the average time spent in each state is given by 1=q p,m (cell cycles). It has been observed that the switching from a stationary to motile state (and back) does not occur faster than on the time scale of one hour [13]. This gives an upper limit on the transition rates, which since time in the model is measured in cell cycles, is given by q p,m v24.
The motility rate is set to n~5 cell cycle {1 . This means that a motile cell on average moves one lattice site, i.e. Dx~20 mm, in a time 1=n~1=5 cell cycles, which gives a linear velocity of 100mm cell cycle {1 , that lies within experimentally determined values of 34mm=24 hours [14] and 500mm=24 hours [13]. The rate of apoptosis is set to the value m~10 {3 cell cycle {1 , which is small compared to the other transition rates in the model.

Simulations
Our concern is the influence of the microscopic cell-level parameters on the growth rate of the tumour as a whole, and we will therefore measure the size the tumour after a fixed time for a given set of initial conditions, as a function of the phenotypic switching rates. More specifically we will measure the tumour mass (the total number of cells), and also later, quantify the rate of spread by measuring the velocity of the tumour interface. The precise initial condition of the model has little impact on the longterm rate of spread (data not shown), but in line with the clonal origin of cancer we initialise the model with a single cell (in the Pstate). All simulations of the IB-model are carried out using the commonly employed Gillespie algorithm [33]. Figure 2 illustrates the results of simulating the model in two dimensions for T~25 cell cycles when (q p ,q m )~(20,10) in three different ways. Panel (a) shows the result of a single simulation, where P-cells are coloured blue and M-cells are red, (b) shows the results of the model averaged across a large number of realisations and gives the occupancy probability Q(i,j) of finding a cell at location (i,j) on the lattice, and (c) shows a slice through this function Q(i,j~50). This figure gives us a general idea of the growth dynamics of the model. The tumour grows with a radial symmetry, and exhibits a solid core, while the tumour margin is diffuse and somewhat rugged. Please note that the time span considered in this simulation is smaller than the time scale of actual glioblastoma growth, which usually occurs on the time scale of months to years, but still sufficient to investigate the dynamics of the model.
In order to quantify the dependence on the phenotypic switching rates we measured the tumour mass at T~50 in the parameter range 0vq p,m v30. The results are displayed in figure 3a and show a strong dependence on the two parameters. For q m~0 all cells are in the proliferative state, and as expected the mass is independent of q p . The other extreme where q p~0 gives rise to tumours with a zero mass, which occurs since the motile cells cannot multiply and eventually die off due to the small but non-zero apoptosis rate m. These results are intuitive, but what is more interesting is that tumour cells with intermediate switching rates are the ones that give rise to the largest tumours. Although migratory behaviour does not directly contribute to an increase in the number of cancer cells it has the secondary effect of freeing up space which accelerates growth compared to the tumours dominated purely by proliferation (q m~0 ). The results suggest that for each q p w0 there is a q m =0 which gives a maximal tumour growth rate. These results also hold for the more biologically plausible 3-dimensional case (see figure 3b). Although the maximal tumour mass seems to occur for a smaller q m , and the region of parameter space giving rise to small tumours is considerably larger (upper left region), the qualitative behaviour is similar. The implications of the observation that q m influences tumour size in a non-monotone way will be discussed later, and we will now proceed to an analytical treatment of the problem.

Derivation of continuum model
In an effort to get a deeper understanding of the somewhat unintuitive relationship between tumour growth rate and phenotypic switching rates we will derive a set of two coupled PDEs which will serve as an approximate way of describing the time evolution of the occupancy probability Q(i,j) (see figure 1b and c). For the sake of clarity we will however constrain the derivation to a one-dimensional system. In fact, radially symmetric travelling wave solutions with constant velocity do not exist for d §2, but instead the velocity of the front depends on the local curvature. However, for large enough times the interface of the circular (spherical) solution has almost zero curvature and its dynamics is  well approximated by the one dimensional solution. A further simplification in our derivation is that we assume the occupancy probabilities of neighbouring sites as independent, which in practice means that we for example allow ourselves to write: Pr(site i empty and site iz1 occupied) = Pr(site i empty)|Pr(site iz1 occupied).
The derivation is carried out in two steps: firstly, a set of coupled master equations, for the two sub-populations, are derived by considering the processes which alter the occupation probabilities at a given site, and secondly these master equations are approximated by a set of PDEs. In brief, the second step is achieved by identifying the on-lattice master equations with a set of coupled PDEs, which when discretised on the length scale of the lattice spacing, equal the master equations. The full derivation is given in Methods and results in the following system of coupled PDEs: Lm Lt~D n ((1{p) Here p(x,t) denotes the density of proliferating cells, and m(x,t) that of the motile cells. In equation (3) we recognise the first term as a diffusion term, modulated by a density-dependent prefactor and the second term as a logistic growth term. The remaining terms are due to the switching between the subpopulations and to apoptosis. In the equation for the motile cells (4) there is also density-dependent diffusion, but of a different type. This is typical of a two species size exclusion process [34], and contains the second-derivative of both species. The values of the diffusion constants are D a~a =2 and D n~n =2, and depend crucially on the choice of spatial scale, which for simplicity is chosen to be that of the cell size Dx. If a coarser spatial scale is considered then the diffusions constants would have to be scaled accordingly (see Methods for details). We will now proceed to investigating the properties of this system of PDEs through both numerical solutions and analysis.

Travelling wave solutions and their velocity
The first question one might ask about a system of equations that presumably describes tumour growth is if it exhibits tumour invasion and hence travelling wave solutions, and further how the model parameters influence the wave speed, i.e. the velocity of the invading tumour front. The results from the IB-model (figure 2 and 3) suggest that the switching rates q p,m strongly influence the tumour mass, and hence we expect them to also have an effect on the speed of invasion.
In order to investigate this, we first solved the continuum model (3)-(4) numerically (which actually corresponds to reverting to the master equations eq. (8)-(9)), for a range of parameter values, in the domain xw0. The initial condition was set to p(x,0)~exp ({bx), b~10 cell width {1 and m(x,0)~0, meant to represent a situation where a tumour is initiated by a small number of proliferating cells (p(x,t)) and no migratory cells (m(x,t)). In fact the balance between p and m in the initial condition is largely irrelevant for the long-term dynamics of the model, the exceptions being the extremes (q p ,q m )~(1,0) and (0,1), when flow between the phenotypes is unidirectional or completely blocked. The boundary conditions of the domain were set to no-flux.
The results can be seen in figure 4 and shows the occupancy probabilities after T~40 and 50 cell cycles.  particular case (q m = 0), when there is no random motion within the cell population. Migration of the cells tends to break up the correlations that build up as the tumour is growing, and as we later shall see, the continuum approximation works better when the cells are more motile.
The observation that the numerical solutions are stationary in a moving frame suggests the existence of travelling wave solutions. In order to close in on these solutions, and get an estimate of their velocity, we will make use of the travelling wave ansatz: p(x,t)~P(z) and m(x,t)~M(z) with z~x{ct, where c is the velocity of the interface. The problem of determining how c depends on the model parameters is solved by applying phasespace analysis (see Methods), and boils down to a four-dimensional eigenvalue problem, namely to find the smallest c such that the eigenvalues of the Jacobian all have imaginary part equal to zero. This problem is analytically intractable, but provides us with a numerically easy way of determining the velocity.

Influence of the parameters on the wave speed
Although phase-space analysis does not yield an analytic closedform expression for the wave speed c, it still provides us with a computationally simple way of determining the velocity of the tumour margin in the model: for a given set of parameter values we start by setting c~0 and calculate the eigenvalues of the Jacobian (18) (or equivalently the roots of the corresponding characteristic polynomial P(l)). If not all eigenvalues are real we increment c slightly and reevaluate the eigenvalues. This procedure is terminated as soon as we find all eigenvalues real, and the value of c for which this occurs corresponds to the wave speed for those parameter values.
In order to test the validity of the wave speed analysis we compared the wave speeds obtained in the continuum and IB models with those from the phase space analysis. For the continuum model an estimate of the wave speed was obtained by, from the initial condition p(x,0)~exp ({bx) (for proliferating cells), b~10 and m(x,0)~0 (for migrating cells), integrating the equations (3)-(4) for 200 time steps (cell cycles). From these solutions we estimated the velocity of the front by measuring the position of a reference point x c , defined as the point where p(x c ,t)zm(x c ,t)~1=2, as a function of time. The comparison between the speed of propagation in the numerical solution and the wave speed obtained from the phase space analysis is shown in figure 5. The agreement is fairly good and the discrepancies are probably due to error in integration and the deviation in the numerical solution from a perfect travelling wave, which from a given set of initial conditions, is only attained in the limit t??. However, since we are interested in biologically relevant scenarios the time frame considered is reasonable.
When it comes to the IB model, we have to take into account the stochastic nature of the model, and therefore need to estimate the average margin velocity from a large number of simulations (100 independent realisations). Each simulation was started with a single P-cell at the center of the lattice and the model was simulated for 100 time steps (cell cycles). In each time step the location of the cells was recorded and from this we calculated the occupation probability Q(i,j,t) of finding any cell at location (i,j) at time t. The wave speed was then approximated by taking the average propagation speed of Q(i,j,t) in the i{ and j{direction (as in the PDE case). In comparing with the two-dimensional simulations we need to rescale the diffusion coefficient D n ?D n =2, since cell movement occurring tangential to the two-dimensional front does not contribute to its propagation. The result can be seen in figure 6, which shows that the analytical result is in good agreement with the discrete individual-based model. The disparity between the IB-model and the analytic answer is largest for small q m , when the dynamics are dominated by proliferation. This is to be expected since for larger q m the movement of the cells decorrelates the sites, and hence our assumption about site independence is closer to truth. The analytical results recapitulates the non-monotone dependence on q m and using this method we found that the largest tumours occur when q max m &q p =2, i.e. when the ratio between the switching rates is 1:2. Naturally the other model parameters also affect the rate of tumour invasion (see figure 7). Increasing the proliferation rate a leads initially (for small a) to an increase in velocity according to c* ffiffi ffi a p , while for a 2 there is a cross-over to a linear dependence with c*ka, with k&0:7. The motility rate n also influences the wave speed in a non-linear way according to the relation c* ffiffi ffi n p , which holds for all nw0. Finally, increasing the rate of apoptosis m, as expected, decreases the wave speed, and does so in a non-linear way. Actually the dependency on m looks very much like that of a second-order phase transition, where the derivative dc=dm diverges at a critical point m c &0:675 cell cycle {1 , and we have for mvm c that c*(m c {m) c (see inset of figure 7c). We observed that the critical apoptosis rate m c , above which no travelling wave solutions exists and hence the tumour disappears, depends on the other parameters of the model, but that the critical exponent c~0:5049+0:0004 is independent of the other parameters.

Phenotypic switching affects growth by altering the tumour interface composition
Our model gives considerable insight into the dependency between five cell-level parameters (switching rates q p and q m , motility rate n, proliferation rate a and rate of apoptosis m) and the macroscopic dynamics of tumour growth and invasion. Focusing on the impact of the phenotypic switching rates we showed that tumour cells with a small q p and large q m (see fig. 3) give rise to small tumours (low c) while those characterised by a large q p and intermediate q m grow into large tumours (high c). To see why this is the case, consider a one-dimensional growth process in which the tumour expands in a narrow channel. If q m~0 , then the tumour expands only through proliferation of the cells at the interface (since interior cells cannot divide), and the interface thus moves with velocity a, equal to the proliferation rate. If q m =0 then cells at the interface spend some time in the motile state, freeing up space and allowing previously blocked cells to proliferate. This process increases the interface velocity, but it is also clear that a large q m has a negative effect on tumour growth, since if q m &q p fewer cells are in the proliferative state and can thus take advantage of the space created via cell migration. From this line of reasoning it is clear that the tumour interface velocity will depend on q m in a non-monotone way. Taken together, our analysis shows that q m and q p affect glioblastoma progression by altering the composition and structure of the tumour interface, and that for each q p =0 the velocity c~c(q m ) attains a maximum, which occurs at q max m &q p =2. The above reasoning, and our model, do however not take into account the effects of mechanical forces between the cells. In particular it is, in real tumours, possible for cells to push one another and hence to divide and move, although there is no free space. This process will most likely lessen the positive effect of cell  migration on tumour growth, but since it has been experimentally established that few cell divisions occur in the core of the tumour due to pressure build-up and hypoxia, we believe that the conclusions of our model still hold to a large extent.
A similar trade-off between proliferation and migration has in fact been observed in the models of Hatzikirou et al. and Fedotov and Iomin [29]. Although a formal comparison with the former model is difficult, the macroscopic equations that Hatzikirou et al. derive show that the number of rest channels (comparable to the likelihood of a cell proliferating), increases the proliferation rate, but at the same time decreases the motility of the cells. In the work of Fedotov and Iomin [29] a similar trade-off is present. Using a continuous time random walk model they showed that if the waiting times in the P-and M-state are exponentially distributed (as in our model) then the margin velocity is non-monotone in the ratio q p =q m and that the maximum velocity is achieved for q p~qm . However it should be noted that their model does not take size exclusion into account, and hence yields an overestimate of the effects migration has on invasion.
A trade-off between proliferation and migration has also been investigated in relation to cancer stem cells and tumour progression by Enderling et al. [35]. They showed that cell migration can lead to the formation of secondary tumour loci, in a process termed self-metastasis, which might accelerate tumour growth, depending on the ratio of migratory and proliferative behaviour. In a related study it was shown that cancer stem cell migration might lead to branched tumour morphology and that it can increase the chance of tumour recurrence [36]. These modelling results together with those presented in this study highlight the importance of cell migration in tumour progression and motivate future experimental studies.

The model recapitulates the wave speed dependency in the Fisher equation
We have also demonstrated that the other parameters in the model affect the speed of invasion. Firstly, the impact of the motility rate and the proliferation rate imply that the wave speed dependence observed in the Fisher equation (1), c* ffiffiffiffiffiffi ffi Dr p , also holds in our system, when equating the diffusion constant D with the motility rate n and the proliferation rate r with a (at least for small and biologically realistic values of a). The Fisher equation has been shown to give an accurate macroscopic description of glioblastoma progression in vivo [22], which also lends support to our model. The correspondence to the Fisher equation is particularly interesting since it allows for a connection between cell level characteristics (n and a) and tissue-scale behaviour, and suggest a means of parametrising our model, not only using singlecell measurements, but also from tissue-level data, such as MRIscans. From consecutive images the position and hence velocity of the invading tumour margin can be determined and compared with the results of the model.

Progression depends on apoptosis rate in a discontinuous fashion
Secondly, we observed a second-order phase transition in the velocity with respect to the rate of apoptosis m. This means that there is a critical apoptosis rate m c above which no tumour can grow and that for mvm c we have c*(m c {m) c , where m c is parameter-dependent, but c&1=2 seems to hold for all parameter values. The discontinuous behaviour of dc=dm is interesting, not only from a theoretical perspective, but also because it implies that if a high enough rate of apoptosis is induced, it may not only retard tumour growth, but in fact lead to regression. However, these results should viewed with caution, since the model would need to be modified and extended in order to properly account for the dynamics of drug delivery and treatment (cf. [37]).

Experimental implications of our model
While data from Farin et al. [13] served as the impetus for our model, we note that a few additional experiments support our modeling assumptions. First, the general observation that glioma cells sampled from invasive fast-growing tumours are characterised by a blend of proliferative and migratory behaviour [2] supports our results, although only in a qualitative way. Second, a recent study on different glioma subclones obtained from the same patient identified a particular cell type as being particularly invasive. Subsequent analysis of proliferation of these clones (determined by Ki67-staining) showed that the most invasive subclone (giving rise to the largest tumours in vivo) had the lowest proliferation rate [38]. Although the subclones were not subject to a motility assay, these results still diminish the importance of cell proliferation in determining tumour growth rate, and future studies that measure both proliferation and migration could be even more useful in this respect.
In order to gain further experimental support for our model, we plan in future work, to measure the five cell specific parameters directly. Such measurements should be possible by applying live imaging microscopy techniques to primary glioblastoma-derived cell cultures. A first application of such measurements could be exploited to develop the model further, to predict progression for an individual patient based on cell-level phenotyping, and to develop chemical compound screens where the impact of a chemical on the model parameters are observed. This might in turn lead to a strategy to define in vivo-relevant compounds more likely to inhibit progression.

Future work
The current model is however far from these highly set goals, and there are a number of extensions that would make the model more realistic. In its current form the model does not account for cell-cell adhesion, which could be incorporated letting the motility rate n be dependent on the neighbourhood of the cell [25,26]. The preferential migration along capillaries and myelin tracts, and the tendency for glioma cells to divide at capillary branch points, is also something that could be included. A further complication is that cancer cells within a real glioblastoma are not identical with respect to their behaviour, but exhibit both genotypic and phenotypic heterogeneity, e.g. cells with a migratory phenotype tend to be located at the tumour boundary whereas dividing cells are commonly found in the main tumour mass, a fact which is not captured by the current model.
Despite this we would still expect the results of our model to hold at least with respect to the large-scale behaviour of the tumour. The real situation is also complicated by the fact that cancer cells are selected for based on their phenotype. One hypothesis which emerges from our model is that selection could drive the behaviour of the cells to the optimal balance between q p and q m , although this hypothesis would require a model that allows for population heterogeneity in order to be tested.
Adding these mechanisms would of course make the model less tractable from an analytical point of view, but this trade-off between simplicity and reality is something that all modellers must deal with.

Derivation of continuum model
Let us consider a one-dimensional lattice indexed by the integers. We let p k (t) denote the probability of finding a P-cell at site k at time t, and equivalently let m k (t) represent the occupation probability of M-cells. The general strategy is to formulate two coupled master equations for the occupation probabilities, which will then be approximated by a set of PDEs, amenable to a wave speed analysis that hopefully will reveal the influence of q p,m on tumour growth.
Let us first consider p k (t). Which are the processes that affect this quantity at a given site?
1. an existing P-cell can die through apoptosis (with rate m) 2. an existing P-cell might switch to an M-cell (with rate q m ) 3. an M-cell residing at the site might switch into becoming a Pcell (with rate q p ) 4. a neighbouring cell might divide placing its offspring in the (empty) considered site (with rate a=2) Summarising all these processes we can write: where the first term is a loss term due to apoptosis, while the second and third term are due to phenotypic switching. The final term is due to cell division from the neighbouring sites, and here we have made use of the independence assumption discussed above. After dividing both sides by Dt and going to the limit Dt?0 we end up with the following expression: In order to simplify the expression and also draw parallels to continuum systems we define a discrete Laplace operator where h corresponds to the spacing of the lattice. Equation (5) can now be written as dp k dt~a If we now turn to the motile cells, the following processes affect m k (t): 1. an existing M-cell can die through apoptosis (with rate m) 2. an existing M-cell might switch to a P-cell (with rate q p ) 3. a P-cell residing at the site might switch into becoming a M-cell (with rate q m ) 4. an existing cell might move away from the considered site (with rate n=2 in each direction) 5. a neighbouring cell might move into the (empty) considered site (with rate n=2) Taking all these processes into account we can write The first three terms can be recognised as apoptosis and switching terms, while the fourth and fifth are due to movement out of and into the site. After a bit of algebra and making use of the discrete Laplacian defined in eq. (6) we get dm k dt~n h 2 2 ((1{p k )D Dm k zm kD Dp k ){(q p zm)m k zq m p k : In summary we have that the time evolution of the occupation probabilities are described by the following coupled equations: dp k dt~a Please note that despite the similarity to PDEs, that describe the changes of a quantity in continuous space and time, these equations are defined on the lattice and describe the probability of finding a cell of a specific type in a certain location. In many instances it is natural to proceed by taking the spatial continuum limit of the discrete master equation(s), but in this case, where we are considering expansion via both cell movement and pure cell division (the case q p~0 ), things are a bit more delicate, because when the size of the cells tend to zero (h?0) so does the contribution of cell division to tumour expansion. In order to achieve a sensible continuum limit, a certain scaling in space and time is required, which implicitly assumes that cell motility occurs on a much faster time scale than cell division [39], something which is generally not the case in the case of glioma biology.
However in order to proceed with the analysis and make use of the toolbox of real analysis we will approximate the above equations with the following PDEs: Lm Lt~D n ((1{p) The motivation behind this choice is that the master equations (8) and (9) are the (space) discretised versions of (10) and (11). The diffusion constants are given by D a~a h 2 =2 and D n~n h 2 =2, where h is the spacing of the lattice, which we for simplicity measure in terms of cell size, and accordingly set h~1. This means that we consider the dynamics on the length scale of a single cell. Please note that the unit of the diffusion constants D a and D n is distance 2 =time, while the unit of the underlying proliferation and migration rates is time {1 . The correspondence between the master equations and PDEs is, however, not rigorous and implies that the analytic results obtained for the PDEs are not in general valid for the master equations, but, as we shall see, still reflect the behaviour of the IB-model to a large extent.

Phase-space analysis
For the sake of clarity let us recapitulate the method applied to the Fisher equation (1) in order to calculate its speed of invasion. The travelling wave ansatz (z~x{ct) turns the Fisher equation into second order ODE in the variable U(z). By introducing the variable V (z)~U'(z) the ODE is turned into a two-dimensional autonomous system. The system has two fixed points (U,V )~(0,0) and (1,0), and by calculating the eigenvalues of the Jacobian (a determinant of partial derivatives) at the two fixed points, one finds that the fixed point at (1,0) (corresponding to the invaded state) is a saddle point (independent of c), while the characteristics of the one at the origin depend on c. For cv2 ffiffiffiffiffiffi ffi Dr p the fixed point is a stable spiral, while for c §2 ffiffiffiffiffiffi ffi Dr p it is a stable node. The heteroclinic orbit connecting the two states goes from (1,0) through the third quadrant (Uw0,U'v0, as in figure 4), and only if the origin is a stable node does it enter the fixed point without spiralling around and attaining negative values of the density U of cancer cells. Negativity would be inconsistent with the non-negative solution of the equation (U(z) §0), and shows that the smallest possible wave speed is given by c~2 ffiffiffiffiffiffi ffi Dr p . In our case the travelling wave ansatz transforms the system of PDEs (10)- (11) to the following set of coupled ordinary differential equations (ODE): where prime indicates derivative with respect to z, and where we have expressed the diffusion coefficients in terms of a and n. Because p and m represent occupation probabilities we seek solutions P(z) §0 and M(z) §0 for all z. In order to perform a phase-space analysis we need to transform the coupled ODEs to an autonomous system by introducing the variables Q~P' and N~M'. This expands equation (12) and (13) into the following four-dimensional system: The boundary conditions reflect the fixed points of the system, which are p 1~( P,M,N,Q)~(0,0,0,0) and p 2~( p ? ,m ? ,0,0), and correspond to the healthy and invaded state respectively. In the limit m?0 the invaded fixed point simplifies to (q p =(q p z q m ),q m =(q p zq m ),0,0), in which case only the relative magnitude of the switching rates q p,m determines the equilibrium occupation probabilities (cf. the values of p and m at x~0 in figure 4b and c). What will help us determine the wave speed c is the characteristics of these fixed points, or more precisely the one at the origin. This method only gives a lower bound on the wave speed, but this minimal c turns out to be the one attained for relevant initial conditions for the Fisher equation [21], although this requires further proof [40].
We will now apply the same kind of reasoning of non-negativity as for the Fisher equation in order to obtain a minimal wave speed c min for our system The roots of this equation l i have, for the biologically relevant parameter values and cw0 non-zero real part, <(l i )=0, implying that the fixed point is hyperbolic and thus that its characteristics are fully determined by linear stability analysis [41]. The aim is now to find the smallest c such that all roots of P(l) have imaginary part equal to zero, since only then are we guaranteed trajectories which do not oscillate around the origin, and remain positive in the variables P and M, which is required since these variables represent non-negative occupation probabilities. Determining the smallest such c turns out to be intractable from an analytic point of view, and we will therefore resort to numerical solution of the eigenvalue problem.