In Silico Neuro-Oncology: Brownian Motion-Based Mathematical Treatment as a Potential Platform for Modeling the Infiltration of Glioma Cells into Normal Brain Tissue.

Intensive glioma tumor infiltration into the surrounding normal brain tissues is one of the most critical causes of glioma treatment failure. To quantitatively understand and mathematically simulate this phenomenon, several diffusion-based mathematical models have appeared in the literature. The majority of them ignore the anisotropic character of diffusion of glioma cells since availability of pertinent truly exploitable tomographic imaging data is limited. Aiming at enriching the anisotropy-enhanced glioma model weaponry so as to increase the potential of exploiting available tomographic imaging data, we propose a Brownian motion-based mathematical analysis that could serve as the basis for a simulation model estimating the infiltration of glioblastoma cells into the surrounding brain tissue. The analysis is based on clinical observations and exploits diffusion tensor imaging (DTI) data. Numerical simulations and suggestions for further elaboration are provided.

Introduction

Glioblastoma multiforme (GBM) is the most malignant of all brain tumors. Apart from heterogeneity and highly invasive behavior, GBM cells tend to infiltrate the surrounding tissue by solely leaving the main tumor mass and traveling long distances inside the brain.1 This diffusive behavior of GBM cells is one of the main causes of tumor relapse after resection. Since infiltrating cells are generally not visible by an magnetic resonance imaging (MRI) and an extensive resection may damage surrounding tissue, small populations of such cells are almost always left at the resection margin of GBMs and are widely believed to drive tumor relapse. Therefore, in order to provide quantitative insight into the nonimageable phenomenon of tumor cell invasion, a number of mathematical models have been developed.2–8 Most of the published mathematical models ignore the anisotropic character of diffusion of glioma cells since truly exploitable tomographic imaging data to this end are rare. In order to enrich the anisotropy-enhanced glioma model weaponry so as to increase the potential of exploiting available tomographic imaging data, we propose a Brownian motion-based mathematical analysis that could serve as the basis for a model estimating the infiltration of glioblastoma cells into the surrounding normal brain tissue. The analysis is based on clinical observations and exploits diffusion tensor imaging (DTI) data. Numerical simulations and suggestions for further elaboration are provided. A realistic model based on the analysis presented could be useful for diagnosis as well as resection and radiotherapy planning.

Mathematical Modeling of Diffusion

The reaction–diffusion equation has been widely used for modeling the diffusion of tumor cells. In Refs. 2 and 5, the authors assume that tumor cells move from regions of higher to lower densities and exploit the reaction–diffusion equation using an additional term corresponding to the net proliferation of tumor cells: where c(x, t) is the concentration of glioma cells at time t and location x, and ρ reflects the net proliferation of glioma cells. D is the diffusion coefficient, assuming different scalar values on regions of white (D) and gray (D) matter, where D. In this model, diffusion of tumor cells is mathematically described by the equation where D = D is the diffusion coefficient in white matter regions, D = D is the diffusion coefficient in gray matter regions of the brain, and x, i = 1, 2, 3, denote the three Cartesian coordinates. This model has been widely used throughout the literature. For a review, see Ref. 9. Numerical treatments of Equation (1) can be found in Ref. 7. In Refs. 10 and 11, the authors have expanded this model, introducing in Equation (1) the diffusion coefficient D in the tensorial form, thereby including anisotropy in their models. Measurements of D in the tensorial form were acquired through DTI techniques, which are described in the following section.

The starting point for developing our model is the observation that Equation (1a) is the Fokker–Planck equation corresponding to the stochastic differential equation where B denotes the standard Brownian motion in R3.12 Given the initial position x of a glioma cell, the distribution of the random variable x (ie, the solution of Equation (1b) at time t) provides the probability distribution over all possible locations of this cell at time t. Given the initial position of a glioma cell in terms of a probability distribution p(x, 0), the probability distribution p(x, t) can be found by the following two equivalent ways: (i) by solving Equation (1a) as a partial differential equation in order to find the evolution of this distribution through time, or equivalently, (ii) by solving the stochastic differential Equation (1b) with initial distribution p(x, 0), in order to find the probability distribution of the random variable x.

In case we desire to study the movement of many glioma cells constituting a tumor, this notion of distribution is interpreted as follows: integration of p(x, t) over an area A of R3 yields the fraction of the total tumor cell population that is located in A.

In the following sections of this article, it will be shown that this point of view provides an intuitive and natural way of modeling anisotropy in the movement of tumor cells.

Diffusion Tensor Imaging

DTI is an MRI technique measuring the diffusion properties of water molecules along specific directions. This is done by defining an ellipsoid in every voxel of a three-dimensional space which mathematically corresponds to a 3 × 3 positive-definite symmetric matrix:

This matrix can be decomposed into the following form: where λ1, λ2, λ3 are the eigenvalues of D (positive, since D is positive definite) and u1, u2, u3 are the corresponding orthonormal eigenvectors. The eigenvalues and eigenvectors of D define an ellipsoid with principal axes lying on the directions of u1, u2, u3 and having lengths , , , respectively (Fig. 1).

Figure 1

The triaxial diffusion ellipsoid. See main text for details.

This ellipsoid, called the diffusion ellipsoid, describes the anisotropic diffusion of dyed water molecules in the specific voxel it refers to. If at the beginning of the observation period, a droplet of water molecules is placed at the center of the ellipsoid, after some time, the front of the diffusing water molecules will form an ellipsoid like the one in Figure 1. This reflects the fact that at a certain location, water molecules do not move toward all directions at equal rates. Diffusion is faster in the larger axis (ie, along the eigenvector corresponding to the larger eigenvalue). For each eigenvector of D, the larger the corresponding eigenvalue, the larger the diffusion along its direction. DTI measurements provide this ellipsoid (actually, the principal axes directions and lengths) for each voxel. For a visual representation, see http://commons.wikimedia.org/wiki/File:DTI-axial-ellipsoids.jpg.

Diffusion tensor MRI has been used extensively for tractography in vivo. It is widely assumed that water molecules tend to move more easily along white matter tracts; therefore, the fiber bundle of local white matter tracts is assumed to be aligned with the largest axis of the diffusion ellipsoid.13 Furthermore, measuring the surface of the front of the diffusing water molecules (ie, the surface area of the diffusion ellipsoid) provides a quantity known as apparent diffusion coefficient (ADC), which has been reported to inversely correlate with local cell density.14,15 Taking into account that brain tumor cells are generally assumed to diffuse toward regions of lower cellular density2,5 and invade the surrounding tissue by moving along white matter tracts,1 the diffusion tensor acquired by DTI has been used to describe the diffusion of tumor cells as well.10,11

Derivation of the Model

Our model will use the diffusion tensor measurements provided by DTI to describe the stochastic movement of tumor cells within the brain tissue. Let us assume that we are given a DTI atlas, consisting of the diffusion tensor in each voxel like in Ref. 10. Since the voxel size is of the order 0.5–1 mm3, we will assume that the diffusion tensor is piecewise constant: the diffusion tensor is constant within a voxel and is defined by the measured diffusion ellipsoid of the particular voxel. Thus, supposing that the observation period of the tensor measurement (see before) is Δt, the movement of a particular water molecule within a voxel can be described by the equation where p(x, x, Δt) is the probability of a molecule starting at x to be at x after time Δt. U is the matrix whose columns are the orthonormal vectors u1, u2, u3 of the diffusion tensor of the voxel, and L is a diagonal matrix with the respective eigenvalues λ1, λ2, λ3 as elements of the main diagonal. Thus, given x, the probability density function of x is a multivariate Gaussian with mean x and covariance matrix ULU. The movement of tumor cells in the same interval Δt can be described by the equation where

Using Equations (4) and (5), Equation (3) can be written equivalently as

The parameter α is positive and rescales the eigenvalues of the tensor, thereby rescaling conformally the axes of the diffusion ellipsoid. This reflects the fact that tumor cells may tend to move along the axes of the ellipsoid, but do so with a different velocity than water molecules. Since tumor cells move slower than water molecules, we typically expect α <1. Equations (3) and (6) again are Gaussian probability densities for the random variable x, ie, the position of a tumor cell after time Δt given that the initial position of the cell is x. Using a standard linear transformation, the random variable x can be equivalently written in the form where Z is a normally distributed random vector with mean x, and whose covariance matrix is the identity matrix in R3. Equivalently, we can write for the random variable x – x, where Z′ is a normally distributed random vector with mean (0, 0, 0) and covariance matrix the identity matrix in R3. This leads us to model the movement in continuous time by the equation where b is a three-dimensional random vector, ie, the distribution of b is normal with zero mean and covariance matrix the identity matrix times Δτ. This is equivalent to the stochastic differential equation where B denotes standard Brownian motion in R3. The matrix U(x) depends on x, and its columns are the orthonormal eigenvectors u1(x), u2(x), u3(x) of the diffusion tensor at x. The matrix L1/2(x) is diagonal, and its main diagonal entries are the square roots of the eigenvalues λ1(x), λ2(x), λ3(x) of the diffusion tensor at x. The parameter a is to be estimated from data. Equation (7) is the model we propose for describing the movement of tumor cells in the brain.

By assuming a twice differentiable diffusion tensor and denoting the probability density function of x by p(x, t), the corresponding Fokker–Planck equation describing the evolution of p(x, t) through time is where β(x) are elements of the matrix where M is the diagonal matrix with main diagonal entries ,· ,· ·12 Denoting by p(x,0) the probability distribution of the position of a cell at time 0, we can estimate the probability that the particular cell will lie at a ball of center x at time t, by integrating the function p(x, t) on that ball. Sampling from the distribution p(x, t) can give us an estimate on how the cells have spread at the end of the time interval [0, t]. Assuming that the diffusion tensor is twice differentiable, solution of Equation (8) could be approximated by numerical methods like finite differences. But DTI provides piecewise constant values for the diffusion tensor, so one should first approximate DTI measurements by a twice differentiable function.

We note that in the case of isotropic diffusion, U(x) is a constant and equals the identity matrix. The eigenvalues λ1(x), λ2(x), λ3(x) are also constants and λ1(x) = λ2(x) = λ3(x) = λ. In this case, equation reduces to which is the form of the diffusion term in Equation (1a) for isotropic diffusion. A comparison of these equations leads to

There are well-known estimates of the diffusion coefficients throughout the literature.2,5,7 The last equation provides means for estimating the parameter a. As seen from this equation, an increased value of the parameter a means a higher diffusion coefficient. Biologically, this corresponds to higher cell motility and a more extended infiltration of glioma cells throughout the brain.

Equation (7) can be written in the form where , , are the one-dimensional components of the three-dimensional vector dB.

Comparing Equations (1b) and (7a) and taking into account equation 2D = aλ reveals the connection of our approach with previous ones. Equation (1b) implies that the diffusion ellipsoid is a sphere and that given its initial position, a cell can move along any unit direction with the same probability. By taking into account the axes of the diffusion ellipsoid as provided by DTI and by varying the quantities , , , our model is able to quantify the preferential movement of tumor cells along specific unit directions.

Previous approaches involve numeric integration of Equation (1a). We will not try to solve Equation (8) numerically. Rather, we are going to assume an initial distribution of tumor cells and numerically simulate sample paths from Equation (7) on the interval [0, t]. These paths are seen as possible cell trajectories, along which tumor cells infiltrate the surrounding brain tissue. We are going to keep track of the random variable x, ie, the position of the simulated cell at time t resulting after each simulation. The value of x is an approximate sample from the distribution p(x, t).

Through our simulations, we will assume that the matrices U(x) and L(x) are piecewise constants. This assumption is imposed by the minimum voxel size of the DTI technique. Interpolation methods could be used to obtain a smooth version of DTI measurements. Owing to absence of an actual DTI atlas of the human brain, we will not do this in this work.

Numerical Simulations

To produce sample paths of Equation (7), we will use the Euler–Maruyama16 numerical scheme: where Z is an R3-valued, normally distributed random variable Z ~ N(0, I3), and Δt is the discretization step. In each iteration, the matrices U(x), L(x) are defined by the eigenvalues and eigenvectors of the diffusion tensor of the voxel in which x lies.

Scheme (9) has order of strong convergence 1/2. This means that if is the solution of Equation (7) at time T as calculated from Equation (9) and x is the actual solution of Equation (7) at time T, then for some positive constant C.16

Rewriting Equation (9) in view of Equation (7a) yields where Z1, Z2, Z3 ~ N(0, 1) that is, Z1, Z2, Z3 are one-dimensional, independent, and normally distributed random variables and Δt is the discretization step.

We will use Equation (9a) for our simulations. It is worth noting that, in view of the previously derived equation, 2D = aλ, since the diffusion coefficient unit is (surface units)/(time units), eg, mm2/h, the quantities , i = 1, 2, 3 in Equation (9a) are length units (eg, mm).

Owing to absence of an actual DTI atlas of the brain, simulations were performed using various synthetic diffusion tensors.

Estimates of D were taken from Refs. 2 and 5. A 10-fold variation of D(D = 10D) was used; we have estimated the quantities , i = 1, 2, 3, using the equation 2D = aλ (Table 1). As for the orthonormal eigenvectors u1(x), u2(x), u3(x) of the local diffusion tensor, we used the following examples:

Table 1

Diffusion coefficient and , i = 1, 2, 3, min/max values.

	MIN	MAX (10×)
D(mm^2h)	0.0054	0.054
Δtαλ_i(x), i = 1, 2, 3, (mm)	0.1	0.32

,

, and

.

In each simulation, we have used a population of 40,000 cells, initially distributed according to a standard (non-skewed) Gaussian distribution with small variance (≤1). We have performed simulations for 400, 800, and 1500 hours using a simulation step equal to 1/60 hours. Since each path can be simulated independently, we have parallelized the algorithm on a quad-core computer. The algorithm took about 170 seconds to perform a 1500-hour simulation. Usage of more processors is expected to improve this time.

Simulation results are shown in Figure 2.

Figure 2

Probable (sample) paths of a cell, with starting and ending points. A cell initially located at the start position will perform random, Brownian-like movements described by Equation (7). Sample paths have been obtained by numerical integration of Equation (7). Performing several such simulations and keeping track of the respective end points, one obtains an image of the distribution of cells throughout the brain, as those depicted in Figures 3–5.

Figure 3A and 3B shows isotropic diffusion. The simulation time is 800 hours. For both figures, the diffusion tensor is constant and its eigenvectors consist of triad (I). For Figure 3A, , i = 1, 2, 3. For Figure 3B, , i = 1, 2, 3. As can be seen, diffusion of cancer cells takes place symmetrically along all axes. In Figure 3B, the larger diffusion coefficient results in cells infiltrating a larger area.

Figure 3

Visualization of isotropic diffusion. See text for details.

Figure 4A and 4B shows anisotropic diffusion. For both simulations, the following synthetic diffusion tensor has been used: For x1 ≥0, the diffusion tensor is defined by , u2(x)=(0,1,0), and , , . For x1 <0, the diffusion tensor is defined by u1(x) = (1,0,0), u2(x) = (0,1,0), u3(x) = (0,0,1) and , , . The simulation time is 800 hours and 1500 hours, respectively. Inspection of Figure 4A and 4B reveals the tendency of the cells to move parallel to when they lie on the x1 ≥0 side and parallel to (1,0,0) when they lie on the x1 <0 side. Figure 4(C) depicts the projection of Figure 4(B) on the x1x2 plane. There is a slight wider spread along the x2 axis on the x1 ≥0 side compared to the one on the x1<0 side. This is due to the higher (0.15 to 0.1) value of on the x1 ≥0 side.

Figure 4

Visualization of anisotropic diffusion. See text for details.

Figure 5(A) and (B) also show anisotropic diffusion. The simulation times are 800 hours and 1500 hours, respectively. The following synthetic diffusion tensor has been used: For x1 ≥0, x3 ≥0, the diffusion tensor is defined by , u2(x)=(0,1,0), and , , For x1 ≥0, x3 <0, the diffusion tensor is defined by , u2(x)=(0,1,0), and , , For x1 <0, the diffusion tensor is defined by u1(x) = (1,0,0), u2(x) = (0,1,0), u3(x) = (0,0,1) and the quantities , , . In view of Equation (9a), there is a clear tendency of the cells to move parallel to eigenvectors with greater coefficients. In Figure 5(C), we show the projection on the x1x2 plane. Cells on the x1 < 0 side appear with a smaller density and a wider spread along the x1 axis than those on the x1 ≥ 0 side. Again, this is partially explained by the higher value of (0.2 to 0.15 mm) on the x1 <0 side. Cells on the x1 <0 side tend to move more easily parallel to the x2 axis. Additionally, by comparing the projections on the x1 axis of the vectors of both x1 ≥0, x3 ≥0 and x1 ≥0, x3 <0 regions with the value of on the x1 <0 region, we see that in the x1 <0 region, cells tend to move more easily along the x1 axis as well.

Figure 5

Visualization of anisotropic diffusion for a different diffusion tensor. See text for details.

Discussion

We have modeled the movement of glioma cells in the brain using a stochastic differential equation. Our method relies on the assumption that glioma cells perform a Brownian-like motion while infiltrating the surrounding brain tissue. Our starting observation has been that many popular models describing the same phenomenon do actually imply the same hypothesis; the well known reaction–diffusion Equation (1a) constitutes the Fokker–Planck equation for a simple stochastic differential equation with zero drift and a constant diffusion term. However, it is perceivable that Brownian motion-based models may not be adequate for describing cell movement. Cell movement is affected by several mechanisms, including directed cell motility along gradients of chemoattractant molecules (eg, growth factors, nutrients, cytokines, chemokines) or along gradients of extracellular matrix molecules.17 These mechanisms are known as chemotaxis and haptotaxis, respectively.18,19 To the best of the authors’ knowledge, in vivo measurements concerning such phenomena are not yet feasible. Our approach provides a basis for modeling the random motion of tumor cells, which can be extended and readjusted if such measurements become available. For example, knowledge that tumor cell motility is affected by certain chemotactic or haptotactic gradients can be readily modeled by introducing appropriate nonzero drift terms to Equation (7), each one corresponding to a specific gradient.

We have not included proliferation of cells in the analysis as yet. The latter can be introduced by allowing the sample paths to branch during the time of the simulation in the following way: if a cell located at position y divides, sample paths can be produced from Equation (7) with initial position y, thus producing probable paths for the cells resulting after mitosis. Undoubtedly, this approach implies larger computational burden. However, as previously stated, our model and its simulation algorithm are suitable for parallelization. We leave this for future work.

Since the model is hypothesis oriented, we provide some discussion on how its performance can be assessed, provided that suitable MRI and DTI are available. As described in the previous section, given an initial distribution over possible locations of tumor cells, the model actually produces approximate samples from the distribution p(x, t) of possible locations x of tumor cells at time t. We denote these samples by x, i = 1, …, N. The actual distribution of tumor cells q(x, t) can be deduced by a three-dimensional MRI. We let each voxel correspond to a cube and consider the σ-algebra Σ produced by these cubes in R3. Cellular densities of each voxel are normalized to sum to 1.

We use the samples x, i= 1, …, N to build an approximation of the distribution p(x, t) as follows: where 1 is the indicator function of the set A.

The total variation distance between distributions and q(x, t) is

Small variation distance between measured and simulated distributions implies good model performance. Minimization of this distance may also be used for parameter estimation.

Conclusions – Future Work

We have proposed a novel mathematical treatment for simulating the invasion of glioma tumor cells into the brain. The method consists of simulating probable paths that the tumor cells may follow, using a stochastic differential equation. An inspection of the results obtained so far shows that our approach satisfies basic qualitative characteristics of anisotropic glioma cell diffusion. Therefore, it will be further elaborated – including, ie, the explicit modeling of cell proliferation – in order to finally serve as a tomographic data adaptable diffusion model of glioma invasion and thus eventually improve the simulation of the corresponding biological phenomena.

The proposed model provides a distribution of the glioma cell population throughout the brain. The latter can lead to an estimate of local cell density in any position within the brain. By discretizing the anatomic region of interest through the utilization of a three-dimensional cubic mesh, we can have an estimate of glioma cell density in each cube.

The glioma cells residing within each elementary cube of the mesh can be clustered into subcategories depending on their cell cycle phase status and metabolic activity. At this point, the discrete-entity, discrete-event approach described in Ref. 17 can be recruited in order to simulate the response of glioma to treatment, such as radiotherapy and/or chemotherapy.