A Brownian dynamics tumor progression simulator with application to glioblastoma.

Tumor progression modeling offers the potential to predict tumor-spreading behavior to improve prognostic accuracy and guide therapy development. Common simulation methods include continuous reaction-diffusion (RD) approaches that capture mean spatio-temporal tumor spreading behavior and discrete agent-based (AB) approaches which capture individual cell events such as proliferation or migration. The brain cancer glioblastoma (GBM) is especially appropriate for such proliferation-migration modeling approaches because tumor cells seldom metastasize outside of the central nervous system and cells are both highly proliferative and migratory. In glioblastoma research, current RD estimates of proliferation and migration parameters are derived from computed tomography or magnetic resonance images. However, these estimates of glioblastoma cell migration rates, modeled as a diffusion coefficient, are approximately 1-2 orders of magnitude larger than single-cell measurements in animal models of this disease. To identify possible sources for this discrepancy, we evaluated the fundamental RD simulation assumptions that cells are point-like structures that can overlap. To give cells physical size (~10 μm), we used a Brownian dynamics approach that simulates individual single-cell diffusive migration, growth, and proliferation activity via a gridless, off-lattice, AB method where cells can be prohibited from overlapping each other. We found that for realistic single-cell parameter growth and migration rates, a non-overlapping model gives rise to a jammed configuration in the center of the tumor and a biased outward diffusion of cells in the tumor periphery, creating a quasi-ballistic advancing tumor front. The simulations demonstrate that a fast-progressing tumor can result from minimally diffusive cells, but at a rate that is still dependent on single-cell diffusive migration rates. Thus, modeling with the assumption of physically-grounded volume conservation can account for the apparent discrepancy between estimated and measured diffusion of GBM cells and provide a new theoretical framework that naturally links single-cell growth and migration dynamics to tumor-level progression.


Introduction
Tumor modeling is rapidly developing as an approach to understand and predict cancer progression [1][2][3][4][5]. Ideally, a tumor progression model would be consistent across a wide range of spatial-temporal scales [6] to predict how molecular-level perturbations, resulting from either mutation or therapeutic intervention, would affect tumor-level progression and, ultimately, patient outcome.
Many computational cancer models of progression have been established, each with variable predictability and applications. The simplest model describing tumor progression is an exponential growth model. While this type of model is helpful for predicting proliferative growth of benign tumors [7], it lacks spatial information and therefore does not capture the invasive behavior of malignant tumors. Similarly, other non-spatial models have been explored to describe tumor growth, such as the Gompertz model or a power law function [8], but have the same limitations of only describing changes in cell number and not tumor spread.
A common simulation method that captures spatial information is the reaction-diffusion (RD) partial differential equation approach [2]. In some cases, this approach also includes convection to account for chemotactic gradients, i.e. simulations are conducted using RD-convection equations [9]. RD models are based on the continuity equation that describes mass transfer in a differential control volume and can be used to predict changes in particle concentration over space and time. Assuming the simplest case of exponential growth, this RD equation is written as ∂C ∂t = D∇ 2 C + pC (1) where C is the concentration of cells, D is the singlecell random motility coefficient (units: m 2 s −1 ) modeled as random Fickian diffusion, and p is the net proliferation rate of cells (units: 1/s). A fundamental limitation of this approach is that it allows cancer cells to achieve unrealistically high concentrations, higher than the physical dimensions of a cell would allow (~10 µm radius). Following equation (1), in all cases when D ⩾ 0 m 2 s −1 , the tumor will grow exponentially without limit. Therefore, even in the extreme case of non-migrating cells (D = 0 m 2 s −1 ), cells are able to continue proliferating. To address this problem, many have used the Fisher-Kolmogorov partial differential equation to describe tumor growth [9,10], which assumes logistic (i.e. sigmoidal) growth with a maximum carrying capacity, k, to limit proliferation when the local tumor cell concentration is too high to permit the creation of a new cell, so that equation (1) becomes Models using this approach have been shown to recapitulate clinically-realistic tissue-level spreading behavior which is often estimated as a traveling wave with linear radial velocity [2,7,11]. These RD models typically do not account for volume exclusion; however an exception is the study by Gerlee and Nelander [12]. More sophisticated versions of the RD model can also account for the normal, non-cancerous cells in the tissues, as well as cells that could potentially be recruited by growth factor secretion [9]. Discrete particle or AB models are more complex approaches that simulate the spatio-temporal dynamics of individual cells rather than mean tumor behavior. These approaches can include complex cellular-level dynamics such as acquisition of mutations or phenotypic changes [3]. While there are a variety of approaches to discrete particle modeling, including cellular automata models [11,13,14], Potts models [15][16][17], or AB models [18,19] (reviewed in (20)), they all have the common attribute of simulating individual cells or cell clusters. Since simulating individual cells at realistic tumor quantities (N ≈ 10 9 ) [13] is computationally demanding, lattice models are often used to reduce computational demands (though this is not a requirement [21][22][23]).
AB models conducted on a lattice often simplify cell migration dynamics by limiting cell movement to only occur as a result of proliferation events based on the availability of a neighboring lattice site. When migration is included, cells seemingly 'jump' from one lattice site to another. One limitation of this approach is that it introduces an arbitrary length scale, namely the choice for the size of the lattice. Often lattice size is approximated as the dimensions of the cell [12,24]; however, lattices can also be sub-cellular to capture more realistic cell morphologies or movements [25,26]. Because AB models of cancer take increasingly more computational power as the cell numbers and number of nodes rise, AB lattice models with sub-cellular lattice sizes are not frequently applied to cancer (such as those mentioned above). Examples of AB lattice models are the model by Hatzikirou et al which includes cell motility modeled as a random walk on a lattice [11], and the model by Khain et al which used a discrete model to simulate in vitro scratch wound assays of glioma [27].
So-called 'off-lattice' AB models allow researchers to investigate the effects of stochastic dynamics of individual cells without the limitations of an arbitrary lattice scale. This type of approach is crucial for understanding the complexity of cell-cell and cell-substrate interactions such as cell-cell repulsion/attraction, cell pushing, or haptotaxis [28][29][30][31].
Another aspect of current AB tumor progression models is that they frequently employ a probabilistic switching between proliferative and migratory phenotypes, representing the so-called 'go-or-grow' dichotomy [12-14, 19,]. An example of an AB lattice model that incorporates phenotypic switching is the recent study by Waclaw et al [13]. While this switch may be substantiated by in vivo data [12], it also adds additional parameters which can be difficult to measure experimentally.
More complex yet are hybrid continuum-discrete models that incorporate environmental variables such as tissue oxygenation or extracellular matrix concentration as a continuum (often modeled with RD equations) into the discrete cell framework [3,20,28]. Some of these models are quite sophisticated and simultaneously account for several variables that contribute to tumor progression. For instance, Alarcon et al used a hybrid model to incorporate the effects of vascularization, blood flow, growth factors, and cellular interaction of normal and cancerous cells into a discrete model of tumor cell growth [6]. Similarly, in a model of breast cancer growth, Kim et al simulated four distinct cell typescancerous epithelial cells, normal epithelial cells, fibroblasts, and myofibroblasts-at single-cell resolution and incorporated growth factor secretion and mechanical interactions between the tumor and surrounding tissues [22]. These hybrid models are highly complex and able to link cellular, or even subcellular, events to gross macrolevel behavior. A summary of the various model types mentioned above is provided in table 1.
In simulating the specific case of the brain cancer glioblastoma (GBM), which is characterized by the invasive migration of primary tumor cells, it is especially important to consider the effects of cell migration on progression. This highly invasive tumor is characterized by exceptionally poor survival of the patient, which is largely attributed to the invasive cells which invariably lead to tumor recurrence and death [1,2,32,33]. Even with current standard treatment of combined surgery, radiotherapy, and chemotherapy, invasive cells, which are not eliminated by these interventions, typically cause tumor recurrence in 3-6 months [34,35] and patients with GBM only survive approximately 15 months post-diagnosis [36]. Previous studies that focused on simulating GBM progression recognized the importance of modeling its invasiveness. As such, nearly all models of GBM have used an RD framework to capture invasive cell behavior [20,32,[37][38][39]. Exceptions are the AB model by Zhang et al that simulates 'go-or-grow' dynamics by including EGFR signaling in their lattice model [19] and the cellular automata model by Tektonidis et al that demonstrates that in vitro glioma spreading behavior can be recapitulated when including cell-cell repulsion and density dependent phenotypic switching [14].
The RD approach has been used in GBM research to validate the therapeutic benefit of surgical resection [33] or chemotherapy [40] by comparing expected tumor growth in the absence of treatment based on repeated computed tomography (CT) images to actual patient outcomes. Similarly, efficacy of various therapies has been proven by assuming continued linear growth based on pretreatment magnetic resonance (MR) images of patients' GBMs in the absence of treatment and comparing to patients' progression with treatment [20]. While many researchers have found such simulations to have predictive clinical power, the predictions are of mean cellular behavior and the physical dimensions of single cells, and their presumed inability to overlap, is not strictly enforced, potentially leading to a mechanistic disconnect between actual cell movements at the single-cell level and the tumor-level progression dynamics.
Input parameters to the RD modeling approach such as proliferation rate and diffusion coefficient for GBM are commonly estimated from sequential CT or MRI scans using Fisher's approximation [1,39]. Briefly, Fisher's approximation states that the growth of the tumor margin will be proportional to the square root of the product of proliferation rate and diffusion coefficient (v = 2 √ pD). While there are limitations with MRI, mainly that it is an indirect measure of glioma cell location as discussed previously [41], it is currently the best clinical tool available for tracking in vivo tumor progression in humans, and it is not unreasonable to attempt to use these parameter estimates to simulate human GBM (table 2). However, more recently, ex vivo measurements of single-cell migration in animal models have emerged. In these studies, the trajectories of single fluorescently tagged cells within ex vivo brain slice cultures are measured by light microscopy, and the cell diffusion coefficient (also known as random motility coefficient) is calculated based on the mean squared displacement (MSD) via where D is the diffusion coefficient, t is the time interval, and d is dimensionality of migration. A summary of experimentally measured values for D in the literature is given in table 3. Interestingly, the values of glioma migration estimates based on tissue-level observations (table  2) are generally 1-2 orders of magnitude larger than those measured at the single-cell level in animal models (table 3). Such high diffusion values used in the RD models have been noted in the literature before. In their simulation of in vitro glioma growth, Stein et al found that more realistic input values for diffusion (40-800 µm 2 h −1 ) resulted in poor agreement between in vitro growth and simulated growth, and, in fact mention that their model outputs would be in much better agreement with in vitro data if they increased diffusion coefficient by a factor of 10 [46]. Moreover, it has been noted that estimating D and p parameters based on Fisher's approximation does not make independent predictions, but accuracy of the approximation improves as the number of cells is increased [47] (i.e. estimates made from bulk behavior). Therefore, we sought to understand the cause of this quantitative difference between parameter estimation and direct measurements by examining the assumptions of the RD approach. Particularly, we considered the unrealistic assumptions of the basic mass continuity equation, which assumes particles are present at low concentrations and occupy no volume (i.e. are point-like).
To evaluate these assumptions, we used a Brownian dynamics (BD) approach. Briefly, BD approaches describe the motion of diffusive particles by gradually incrementing time, assigning randomly selected displacement vectors based on mean diffusive behavior to each of N particles, and tracking the resulting particle positions over time [48]. Such methods have been used to simulate nanometer scale events such as protein folding [49] and filament assembly [50]. Altough cells are much larger by comparison, they have also demonstrated random walk migration behavior [51], hence we used the BD simulation method to simulate migration, growth, and division of individual cells in three different models:  Anderson et al [3] Generic tumor growth 3600 Jiang et al [42] Generic tumor growth 3.6 Tracqui et al [40] Anaplastic astrocytoma (estimated by CT) 43 000-46 000 Woodward et al [33] Glioma 1800-5400 Burgess et al [43] High grade glioma (Fisher's approximation) 5400 Low grade glioma (HG Glioma/10) 540 Swanson et al [32] GBM gray matter (Fisher's approximation) 5400 GBM white matter (Fisher's approximation) 27 000 Rockne et al [39] Glioma (estimated by MRI) 600-37 000 Wang et al [20] GBM  Using these models, we explore how the inclusion of non-zero cell volume affects overall tumor progression rates, such that relatively low diffusion coefficients enable realistic tumor margin velocities, and are sufficient to explain the discrepancies between estimated and measured cell diffusivity.

Results
Comparison of glioma migration parameter estimates obtained from estimates of bulk tumor growth (table 2) to single cell migration parameters obtained by fluorescence microscopy (table 3) shows a disconnect between single-cell behavior and RD model-estimated diffusivity. Therefore, we first assessed whether the estimates for cellular diffusion used in RD models were physically realistic. Among the fastest cells are human embryonic mesenchymal stem cells, which have been shown to move along a 1D adhesive track at a velocity of 5.2 µm min −1 and are highly directional in their movements [53]. From this velocity and the relevant time scale of typical directional switching (~10 min), the theoretical diffusion coefficient of the 'world's fastest cell' is approximately 8000 µm 2 h −1 . Moreover, the cell speed of 5.2 µm min −1 is approximately equivalent to the speed of actin polymerization, 0.21 µm s −1 (or 12.6 µm min −1 ) [54], and therefore is likely to be close to the upper limit for possible cell motility via actin-based motility. Consequently, any values of diffusion beyond this theoretical upper limit likely do not represent actual single cell behavior. While some of the estimates in table 2 are within this limit, many are extremely large, particularly for estimates of white matter migration. Therefore, we were motivated to reexamine the assumptions of the RD models of tumor spread used to estimate such large values.
One issue we noticed with adapting an RD framework to estimate tumor progression parameters was the assumption that the diffusing species (i.e. the cells) is present at low concentrations (i.e. is dilute). However, characteristically, tumor cell density is very high, with cells closely packed together. Therefore, we derived the diffusion coefficient for a concentrated tumor environment. At high concentrations, we assume that the gradient in chemical potential (i.e. cell-cell interactions) dominates cell movements instead of the concentration gradient. Estimating the chemical activity as increasing linearly with mole fraction with an intercept at 1 (to be consistent with the low concentration equation), we found that the apparent diffusion due to high cell concentration is at most twice the intrinsic diffusion (see full derivation in the supplement (stacks.iop.org/CSPO/4/015001/ mmedia)).
Since the estimates of diffusion were much greater than twice the measured single-cell values and could not be explained by the low concentration assumption, we next examined the RD assumption that cells can be treated as points in space that occupy no volume. To investigate the consequences of this at the cell and tissue scales, we used a BD approach to simulate how single-cell behavior relates to overall tumor progression with and without this assumption. Since the BD approach simulates single stochastic events, this method is computationally intensive compared to the RD approaches, and therefore one-dimensional (1D) simulations were used to explore basic dynamic features and regimes. Moreover, the primary issue with current diffusion estimates is the use of uncharacteristically large values for migration along white matter tracts, which can be described as a quasi-1D route of migration [55], therefore we believe the 1D model is potentially relevant to glioma disease progression along these routes, and possibly others such as in the perivascular space of blood vessels. This is not to say that 2D or 3D approaches are invalid for modeling glioma progression; on the contrary-in future models we aim to explore these dimensions as well, as they pertain to in vitro cell cultures and macro-scale glioma spreading.
Starting from a single initiating cell, in silico 1D tumors were grown for a duration of time long enough to generate the sixth generation of cells (assuming a constant growth rate and unimpeded division) and positions of all of the cells in the simulation were tracked over time. In total, three 1D models were analyzed. The first model, termed the simple overlapping model, captures the assumptions of the mass conservation equation (equation (1)) by allowing cells to migrate, grow, and divide without regard for the presence of neighboring cells occupying the same space. In other areas of literature the simple overlapping model is referred to as 'branching Brownian motion' [56]. The second model, the overlapping model with carrying capacity (overlapping + cc), exhibits similar rules to the simple overlapping model in terms of diffusive and growth behavior, but does not allow cell division to occur if the cell in question is overlapping with another cell ( figure 1(A)), thus enlisting a local carrying capacity akin to the logistic growth term in the Fisher-Kolmogorov equation (equation (2)). The final model, the non-overlapping model, does not allow cells to occupy the same space at any point in their life-cycle, as diagramed in figure 1(B). A full description of these models is available in the methods section.
A non-overlapping BD model predicts a 'jammed' configuration and a correspondingly abrupt rise in tumor expansion rate By directly comparing the outcomes of the three models, it is evident that the conservation of volume used in the non-overlapping model enhances tumor spreading, as shown in figure 2(A) (see also supplemental figures 1-3 and supplemental movies 1-3). Moreover, proliferation is limited to the outer edge in the non-overlapping simulations (indicated by red dots) which is consistent with histological observations [11,13]. The distribution of final cell positions can be approximated as Gaussian (normal) in the overlapping models, consistent with a diffusive model, whereas the final cell positions in the nonoverlapping model exhibit an approximate uniform distribution, indicating non-diffusive behavior (figure 2(B)). Moreover, the curvature in the cell MSD plots in figure 2(C) is approximately linear in the overlapping models, again indicating diffusive cell movements, while the non-overlapping model MSD plot has upward curvature indicating nearly ballistic behavior. Plots of tumor radius over time (figure 2(D)) show that the overlapping models advance linearly with time. By contrast, the non-overlapping model tumor radius also increases linearly with time until it reaches a jammed configuration. After jamming, the speed of radial growth is increased (i.e. the slope is greater), consistent with previous AB modeling approaches of glioma migration and proliferation [19], however this effect is captured here without the need to assume a 'go or grow dichotomy' and with physically-grounded rules for cellular movements. Together these data indicate that the non-overlapping model enhances tumor spreading compared to the overlapping models with the same input parameters.

Parameter sensitivity analysis reveals synergy between proliferation and migration in the non-overlapping BD model
In each of the three models, the input parameters affect the progression characteristics differently. As shown in figure 3(A), in both of the cell-cell interaction models (overlapping + cc and non-overlapping) the ratio p app /p 0 of the apparent proliferation rate, p app , to nominal proliferation, p 0 , is decreased as p 0 is increased. As expected, the simple overlapping model shows no correlation with increased proliferation, which serves as a kind of positive control simulation. We note that the simple overlapping situation is clearly unrealistic as arbitrarily large numbers of cells can exist within a finite volume in this model; nevertheless, it is interesting to compare to the cell-cell interaction cases for reference. In the cell-cell interaction models, increasing input proliferation rate results in a relative decrease in the apparent proliferation, with the effect being stronger in the overlapping + cc case. Interestingly, the slope p app /p 0 in the non-overlapping model is similar between different diffusion coefficient values, but this is not the case for the overlapping + cc model, especially at the lowest diffusivity values. This is likely because the small diffusivity values do not allow cells to move away from each other to generate When comparing the effects on apparent diffusion, as expected, the overlapping models allow cells to dif-fuse freely, so the apparent diffusion, D app , is equal to the nominal diffusion coefficient D 0 at large values of D 0 ( figure 3(B)). For small values of D 0 , apparent diffusion is increased as p 0 increases because the generation of a new cell artifactually moves the cell centroid the distance of one cell radius (10 µm) during one cell cycle, which is a significant contributor to the movement in The non-overlapping model progresses linearly, similar to the overlapping models, then the slope of the linear progression suddenly increases after the simulation reaches a jammed configuration, as indicated by the yellow arrow. Thus, the non-overlapping model predicts fundamentally distinct tumor progression dynamics, and faster growth for the same input parameter values, as compared to the overlapping models. this regime. What we observe in the non-overlapping model is that central cells at high concentrations are unable to move because they are prevented from overlapping. This results in a 'jammed' configuration where interior cells nearly stop diffusing. At the same time, exterior cells have freedom to move and are biased toward outward movements away from the center of the tumor, making apparent diffusion larger for those cells. In this model, the diffusion ratio (D app /D 0 ) is not dependent on p 0 , but is strongly dependent on D 0 . At low D 0 , especially, the prevention of cell overlap can generate apparent diffusivity, D app , that reaches values up to 65 times that of D 0 , showing that the effect of jamming-induced ballistic behavior can largely account for the large discrepancy between single cell diffusion measurements and estimates of diffusion in literature.
Analysis of the tumor margin velocity of the three models shows a dependence of velocity on both prolif-eration and diffusion ( figure 3(C)). However, the overlapping models show a reduced tumor margin velocity, especially at lower diffusion coefficients. Compariso ns of the plots for the different models shows how a smaller input diffusion coefficient in the non-overlapping model can generate a tumor velocity similar to the overlapping + cc model at larger input diffusion (e.g.  show dependence on D 0 and p 0 , particularly at low diffusion coefficient and high proliferation rate. This is due to the growth and division of cells at early time points, rather than diffusive movement, which generates an artifactual apparent movement, which appears as a significant boost in motility when the nominal diffusion coefficient is small. The non-overlapping model shows little dependence of apparent diffusion on proliferation, but does have a true boost in apparent diffusion due to jamming effects, particularly at low D 0 . (C) Tumor velocity as a function of both input parameters. In most cases, the overlapping + cc model progresses much slower than the non-overlapping model, reflecting the abrupt rise in velocity due to jamming in the non-overlapping case only. Thus, the analysis reveals synergy between proliferation and migration in the non-overlapping model to strongly enhance tumor growth rate, particularly at the lower diffusion coefficients consistent with single cell experimental measurements.
velocity of Fisher's Approximation was also conducted and is presented in supplemental figure 6.

The non-overlapping model predicts intratumoral gradients in apparent cell diffusion rate and clinically-relevant tumor margin sizes and jamming times
Using an AB approach allows motion statistics to be computed for individual cells as a function of position in the tumor, relative to the tumor edge, which makes testable predictions about the spatial dependence of cell trajectories. As shown in figure 5(A), the overlapping models show no pattern in instantaneous velocity, while the non-overlapping model shows a distinct quasi-ballistic front with a zero velocity core ( figure 5(A)). Similarly, single cells in the overlapping models all have normalized instantaneous diffusion measurements centered around unity that do not vary with spatial position. In the non-overlapping model, cells on the ballistic front display superdiffusive/ convective behavior (with values over 1), while cells at the core display subdiffusive/confined diffusion behavior ( figure 5(B)). The length of the convective front in the non-overlapping model varied with input parameters, as shown in figure 5(C). Ballistic front length was largest for low values of p 0 and high values of D 0 , because as p 0 increases, more cells become jammed. In general, the non-overlapping model predicts reasonable edge lengths that are of order 0.1-1 mm, corresponding to the approximate length scale of a diffuse GBM tumor margin [57]. Likewise, the time to reach a jammed configuration depends on the input parameters, as shown in figure 5(D). Jamming occurs more quickly with high proliferation inputs, and depends only slightly on diffusion inputs. In the p 0 = 10/month situation, some simulations did not reach a jammed configuration in the time allotted for simulation, but would likely jam given sufficient time.
Overall, we find that the non-overlapping model requires much smaller diffusion coefficient values to reach the same rate of tumor spreading as the overlapping models because cells follow physical rules and are not permitted to occupy the same space, which in some cases results in the phenomenon of jamming and a ballistic advancing tumor front. Thus, this suggests that the discrepancy of diffusion coefficient values between bulk tumor estimates (table 2) and single-cell measurements from fluorescence microscopy (table 3) could possibly be explained by the fact that the current RD continuum modeling approaches do not account for the effect of incompressible tumor cell volume on spreading. Therefore, the RD models can fit the MRI data, but may require artificially large diffusion coefficients to do so.

Discussion
Several models have been developed to describe tumor dispersion. However, many of these models, particularly those based on RD approaches, require unrealistically large values of individual cellular diffusion coefficients which are estimated from apparent spreading based on CT or MRI (compare tables 2 and 3). In fact, singlecell values for animal glioma migration observed by fluorescence microscopy and modeled as a random diffusive process are 1-2 orders of magnitude smaller than estimates used in many of these simulations. To understand the source of this discrepancy in diffusion coefficients between estimates and measurements, we evaluated the assumptions of the RD approach. We first tested the RD assumption that particles are present at low concentrations and showed that a highconcentration correction based on thermodynamic activity coefficients only resulted in a 2-fold larger apparent diffusion coefficient. Therefore, we also evaluated the assumption that cells can be modeled as volume-less points in space. Using a 1D BD modeling approach that accounts for single-cell stochastic movements, we found that the inclusion of volume has a significant impact on overall tumor spreading dynamics. Specifically, imposing a simple constraint that cells cannot occupy the same space led to the development of a 'jammed' configuration that generated a non-motile tumor core and a quasiballistic tumor periphery. This is consistent with other AB models that have investigated the effects of volume exclusion and found that it gives rise to a convectivelike term [14,30]. Additionally, the non-overlapping correction allows realistic diffusion input parameters to generate faster tumor spreading velocities, and thus we hypothesize that cell jamming contributes to overall tumor spreading in vivo.
What is interesting about the non-overlapping model is that it generates increased spreading behavior without the inclusion of a chemotaxis term [42] or the 'go-or-grow' dichotomy [37], in which cells are either proliferative or migratory. All cells in the nonoverlapping model have the capacity to both proliferate and migrate, and a ballistic tumor front develops naturally from physical volume-exclusion principles. The advantage of this approach is that it is unnecessary to invoke a third parameter to represent the probabilistic switching between the proliferative ('grow') and migratory ('go') phenotypes, which has frequently been added empirically to generate realistic tumor growth [13].
In future studies, it should be possible to test the prediction that jamming has a significant impact on tumor spreading by comparing animal model behavior to computational model predictions. If jamming does significantly contribute to GBM spreading, the non-overlapping model predicts that there would be a difference in motility between cells at the tumor edges and in the tumor core, which can be tested by measuring single-cell trajectories as well as cell location relative to the central tumor. Furthermore, the non-overlapping model predicts that the migration of cells at the edges would be biased in the outward radial direction. These predictions can potentially be assessed in ex vivo animal systems such as those used to measure single-cell diffusion coefficients [56,58].
One potential limitation of the current model is that we characterize the effects of conserved volume in 1D, while real tumors progress in 3D. When exploring the behaviors of our model in 3D, we expect that jamming will contribute to tumor progression in that it will bias migration radially outward; however, since cells are able to move in multiple directions to escape jamming, we expect that it will take more time to achieve a jammed configuration. As the tumor grows radially in 3D, the radius of curvature of the tumor surface, relative to the cell size, will become large and the surface will appear locally to be flat, similar to the 1D Cartesian model presented here. Interestingly, our 1D model may have application to GBM tumor dissemination because GBM cells have been observed to migrate preferentially along linear tracts within the brain such as perivascular spaces [59] or white matter tracts [55,58,[60][61][62]. Therefore 1D invasion could be a clinically relevant mode of tumor spreading in GBM that contributes to its extensive invasiveness, and the current model could describe tumor progression at this mesoscale, and provide a mechanistic link between the single cell growth/movement scale (micrometers/hours) and the clinical scale (millimeters/months).
In terms of simulating GBM tumors, we also recognize that our model lacks the characteristic necrotic regions that are typical of GBM tumors. This could easily be incorporated into the model by the addition of oxygen consumption where more dense regions have to compete for resources and thus are more likely to become necrotic [3]. Furthermore, other models have demonstrated how the inclusion of oxygen consumption generates in silico tumors with more invasive, finger-like morphology [3,34], frequently observed for GBM tumors. Even so, the presence of necrotic regions are unlikely to affect the advance of the tumor margin, as the proximal necrotic cells are presumably in an environment that does not permit cell proliferation or migration.
In comparison to RD modeling approaches to GBM, the non-overlapping model presented here has an advantage in that it can simulate single-cell events. This is important because predictive models of tumor behavior should ideally be physically and biologically consistent across multiple spatio-temporal scales, and RD approaches have been previously shown to fail to capture cell-cell interaction and clustering effects, particularly when the proliferation rate is significant compared to diffusion coefficient [18,63], which is also when we would expect jamming to occur.
Although the BD approach used here is based on single cell events, it is perhaps not the optimal approach for generating patient predictions. A major disadvantage of BD approach is that it is computationally intensive, which may partly explain why the field has previously relied largely on RD approaches. For future studies we propose implementation of the concepts discovered here into more high-throughput computational models, perhaps by the inclusion of a convective term to an RD framework to simulate jamming effects, or by the addition of a concentration-dependent diffusion coefficient such as that used by Harko et al [7]. However, the model by Harko et al used a concentration-dependent D such that, as cellular concentration approaches c max (a carrying capacity), D approaches infinity. In our model we see the opposite-as local concentration reaches its maximum, cells become jammed and can no longer move. Therefore, rather than applying the change empirically, it may be possible to derive an analytical solution [30,64,65] to the jamming dynamics observed in our simulations so that the model retains its relevance at both the cell and tissue scales. In addition, our simulations reveal an early onset of the jammedballistic behavior, which can be regarded essentially as a steady state that will remain locked-in indefinitely, potentially making longer simulations unnecessary (supplemental figure 4). In the future it will be interest- ing to integrate our modeling predictions of jammed movements, which allows for voids between cells and breaking off of cells at the edge, with recent modeling of jamming of epithelial cell monolayers, which models the subcellular interactions as sets of vertices that define shared edges between cells [68].

Conclusions
Many different approaches have been taken to the problem of predicting tumor growth outcomes. One common method has been to adapt the RD framework, used to solve molecular transport problems, to the realm of cancer cell proliferation and migration via a diffusion-like process. However, since in vivo single-cell measurements of diffusion and proliferation have not been reported for human GBM, the parameters for this type of model applied to GBM have been estimated using bulk tumor characteristics of tumor progression from in vivo imaging modalities such as MRI. Within animal models of the disease, however, direct measurements of single cell diffusivity are now available. Despite the different species, we would expect that cell-level measurements between animals and humans would be similar, as cell migration generally employs conserved actomyosin-adhesion machinery [66], and yet there is a large discrepancy between single cell fluorescence microscopy measurements and estimation from MRI clinical data. Therefore, we re-evaluated the underlying assumptions of the RD framework and found a possible origin of this discrepancy, namely that the RD framework assumes individual particles do not occupy volume and therefore are permitted to overlap, while real cells clearly do occupy volume, and therefore overlapping configurations are physically unattainable.
In implementing this adjustment to the RD framework, we found that non-overlapping cells rapidly enter a jammed regime that produces slow-growing interior cells with subdiffusive migration characteristics and fast-growing peripheral cells with superdiffusive, ballistic, migration characteristics. The increased tumor spread due to jamming-induced convection at the periphery allows cells with nominally low diffusive characteristics to generate a fast progressing tumor that still depends strongly on diffusive cell movement on the local micrometer scale. Thus, the non-overlapping model helps to explain the apparent discrepancy between measured and estimated diffusion coefficients, and provides a more realistic in silico tumor simulator based on measurable parameters. In summary, we suggest that non-overlapping behavior should be implemented in future models of tumor progression to generate physically realistic and quantitatively accurate models that are predictive at multiple spatio-temporal scales.

1D BD simulation
Three distinct BD models were created using a variable time step Monte Carlo framework, which we refer to as the simple overlapping model, the overlapping plus carrying capacity (cc) model, and the non-overlapping model. Simulations were seeded with one initial 'cell' modeled as two completely overlapping spheres. Time (t) was incremented based on the nominal time increment (τ = 1 min) divided by the current number of cells, N cells (t) as shown in equation (3). This was to account for the fact that the probability of selecting any one particular cell is decreased proportionally by the current number of cells.
At every iteration (i) of the simulation, a uniformly distributed random number was used to select a single cell j for movement.
A second normally distributed random number was generated to determine the magnitude and direction of the movement of cell j based on the input diffusion coefficient, D, where x j (i) is the cell centroid position.
In the non-overlapping model, movements were rejected if the location of the centroid at time t(i + 1) caused the cell volume to overlap with any neighboring cell volumes. The simulation would then attempt up to 100 times to find a suitable location to move, and if it failed to do so, would return to its previous position at time t. The number of attempts is an adjustable parameter based on the trade-off between accuracy and speed of simulation; in our simulations we chose to be relatively conservative in this regard, favoring accuracy over simulation speed. In the overlapping model, however, this overlap check was not performed and cells were permitted to occupy the same space when diffusing. Another uniformly distributed random number was then used to select a cell for growth using equation (5) (above). To simulate cell growth, the two overlapping spheres, defining the volume of a single cell, gradually pushed apart such that, when unimpeded by other cells, the total cell volume increased linearly and volume was exactly doubled at time 1/p 0 , where p 0 is the nominal proliferation rate. At each iteration, the growth volume was proportional to the time step determined in equation (4) (above).
For every time point, the distance between two cell halves, or the centers of the two spheres which together represent a single cell (δ), was determined based on the total cell volume (as determined in equation (7)) and the volume of the spheres' overlapped region, known geometrically as the spherical cap (V cap ), where h is the cap height.
In the non-overlapping model, cells were only able to grow if the growth step did not result in overlap with any neighboring cells. In both the non-overlapping and overlapping + cc models, after the cell reached a volume of 2(4/3)πr 3 , where r = 10 µm, a check was performed to ensure the resulting daughter cell would not overlap any neighboring cells, and if the check was passed, the cell would divide and form 2 individual sets of completely overlapping spheres. This essentially imposed a physically-based, local carrying capacity limit where the probability of cell division decreases with increasing local cell density.

Data analysis
A minimum of 4 simulations per input parameter set was used for data analysis. Simulations were carried out to the number of iterations necessary for the initial cell to generate its sixth division event based on the nominal proliferation rate (i.e. the time it would have taken for cells to reach the 6th division if divisions were not impeded due to cell overlap). Apparent proliferation rate was calculated based on the actual number of division events that occurred within the given time. Apparent diffusion was calculated based on the MSD of individual cells over time, using equation (3), and MSD was calculated using the overlapping method [67]. Tumor velocities were calculated based on the distance of the most distal cells to the tumor center over time. Single-cell velocities were calculated as (x(t) − x(t − Δt))/Δt. The Δt value used to calculate instantaneous velocities was 8.3 h, consistent with typical brain slice imaging experiments. (See supplemental figure 5 for the effects of varying Δt.) Similarly, instantaneous diffusion coefficients were calculated as (x(t) − x(t − Δt)) 2 / (2Δt). The length of the ballistic front in the nonoverlapping model was estimated based on the linear fitting of the instantaneous velocity plots. Values less than 10% of the maximum were eliminated from the fit to remove jammed cells and the ballistic front lengths were taken as the x-intercept of the fit. Time to jamming was estimated based on when the center cell reached an instantaneous velocity of 0, with a tolerance of 10 nm/month, in the non-overlapping model.