### Model construction

We propose a model that combines the diffusion tensor driven migration of cells along white matter bundles (Eq. 2) with logistic tumor cell proliferation and with a novel necrosis term which is activated once cell density has surpassed tissue supporting capacity:

$$ \frac{\partial c\left(\mathbf{x},t\right)}{\partial t} = \left\{ \begin{array}{lr} \nabla \cdot \left(\mathbf{D}\left(\mathbf{x} \right)\nabla c \right) + \rho c \left(1-{c}/{c_{m}}\right) & : \max \left[c\left(\mathbf{x}\right)\right] < \tau \\ \eta c & : \max \left[c\left(\mathbf{x}\right)\right] \geq \tau \end{array} \right. $$

(3)

We introduce the necrosis rate, *η*∈(0,1), which produces an exponential decay of tumor cells once cell concentration exceeds a threshold, *τ*. Here, max[*c*(**x**)] denotes the maximum cell concentration at position **x** across time; in keeping with empiric observations, our model presumes a physiologic change that maintains the necrotic state even as cell death causes the tumor cell density to fall below *τ*.

We utilize MR imaging features of each GBM lesion to estimate model parameters, thus producing a customized, tumor-specific predictor of future growth and proliferation. Specifically, we measure tumor radii on serial T_{2}-weighted and contrast-enhanced T_{1}-weighted images, which have been hypothesized to correspond to cell densities of 0.16*c*
_{
m
} and 0.80*c*
_{
m
}, respectively [18, 19]. This relationship is illustrated schematically in Fig. 1. Given the tumor radial velocity (*v*) and cell density gradients thus estimated, Fisher’s relationship (\(v=2\sqrt []{\rho D}\)) enables the calculation of the proliferation parameter, *ρ*, and cell diffusivity, *D* [5]. We note the physiological implications of this choice of sequences. The enhancing tumor corresponds to areas of angioneogenesis with capillary leak. The T_{2} margin represents infiltrating nonenhancing tumor and associated vasogenic edema; the notion of subthreshold tumor beyond the area of T_{2} abnormality has also been described and is important in radiation therapy planning [20].

As has been previously described, to account for the observed strong preferential migration of GBM along white matter tracts, we consider **D** as the tumor cell diffusion tensor, derived via a differential scaling of the eigenvalues of the symmetric, positive definite spin diffusion tensor calculated from diffusion-weighted MRI [13, 21]. Furthermore, we scale **D** such that the mean cell diffusivity matches the *D* estimated from the serial MR images as described above.

We fix *c*
_{
m
}=10^{5}cells/mm^{3} as in previous work [13]. The necrosis threshold, *τ*, is chosen as a fraction of *c*
_{
m
}, and is customized for each tumor based upon the earliest detection of necrosis on contrast-enhanced T_{1}-weighted imaging; lesions which exhibit earlier necrosis have lower values for *τ*. Finally, *η* is also empirically selected for each tumor based on the width of the enhancing rim on T_{1}-weighted imaging such that lesions with narrower regions of enhancement have lower *η* (more rapid necrosis).

### Model validation

To validate our model, we simulate tumor growth in the brain of a healthy adult volunteer using parameters derived from the imaging data of a rare patient with GBM and serial preoperative imaging. We computationally seed tumor cells in a voxel within the brain of our volunteer corresponding to the location of the observed tumor isocenter, and we study the simulated GBM progression over time. Allowing for intersubject anatomical variability—which we attempt to minimize by using an age and gender matched volunteer—we expect that an ideal model would reproduce the observed tumor growth very closely. Thus, the ability of our simulation to recreate the lesion provides insight into the validity of our technique. We qualitatively compare the observed tumor profile to those produced by the proposed and prior models. We further perform a quantitative comparison by nonlinearly registering the two individuals’ brains and computing the mutual information between the simulated tumor profiles and the contrast-enhanced T_{1}-weighted image of the observed tumor. In this context, the mutual information (MI) is a measure that quantifies the amount of information that the simulated cell density shares with the post-contrast T_{1} signal intensity (*s*). If *p*(·) denotes a marginal probability density function and *p*(·,·) indicates a joint probability density function, the mutual information is given by:

$$ \text{MI}\left(c;s\right) = \sum_{C \in c} \sum_{S \in s} p\left(C,S\right) \log_{2} \left[\frac{p\left(C,S\right)} {p(C) p (S)} \right] $$

(4)

In addition, to examine the sensitivity of our model to the position of the initial seed voxel, we repeat the above experiment with a fixed *ρ* and *D*, but we displace the initialization voxel in each of four cardinal directions from the original seed. We qualitatively describe the resulting tumor progression profiles in relation to the original, and we also quantify the differences between profiles initialized with slightly different seeds using the mutual information metric as described above.

All images were obtained using a 3 T scanner with typical acquisition parameters for T_{1}-weighted (T_{R} = 487 ms, T_{E} = 16 ms, 5 mm slice thickness), T_{2}-weighted (T_{R} = 3670 ms, T_{E} = 93 ms, 5 mm slice thickness), and diffusion-weighted (T_{R} = 4100 ms, T_{E} = 95 ms, 4 mm slice thickness) sequences. Diffusion-weighted images were acquired over 64 distinct gradient directions at a *b*-value of 1000mm/s^{2} with 3 repeat acquisitions which were averaged to maximize the signal-to-noise ratio. Following correction for eddy current distortions and head motion, diffusion tensors were estimated from diffusion-weighted MR images using the Diffusion Imaging Reconstruction and Analysis Collection (Laboratory of Neuro Imaging, Los Angeles) [22]. Modeling of tumor progression was performed using custom implementations written in MATLAB®; (MathWorks®;, Natick). Equation 3 was solved using a finite difference scheme with a small discretization timestep (*Δ*
*t*=1day) to maintain numerical stability. Registration was performed using the FNIRT utility from FSL (Centre for Functional MRI of the Brain, Oxford) [23].