A Multiscale Model for Glioma Spread including Cell-tissue Interactions and Proliferation

Glioma is a broad class of brain and spinal cord tumors arising from glia cells, which are the main brain cells that can develop into neoplasms. Since they are highly invasive they are hard to remove by surgery, as the tumor margin it most often not precisely enough identifiable. The understanding of glioma spread patterns is hence essential for both radiological therapy as well as surgical treatment. In this paper we propose a multiscale model for glioma growth including interactions of the cells with the underlying tissue network, along with proliferative effects. As in [12], we assume that cancer cells use neuronal fibre tracts as invasive pathways. Hence, the individual structure of brain tissue seems to be decisive for the tumor spread. Diffusion tensor imaging (DTI) is able to provide such information, thus opening the way for patient specific modelling of glioma invasion. Starting from a multiscale model involving subcellular (microscopic) and individual (mesoscale) cell dynamics, we do a parabolic scaling to obtain an approximating reaction-diffusion-transport equation on the macroscale of the tumor cell population and perform numerical simulations based on DTI data in order to assess the performance of our modeling approach.


Introduction
Gliomas are malignant primary tumors in the human brain with a poor prognosis [37].The common treatment approach involves surgical resection, usually followed by radiotherapy.Thereby, the assessment of the tumor margin is, however, a challenging issue, due to the high invasiveness of gliomas and their finger-like, often disconnected patterns.Mathematical models aim at providing enhanced information about position, shape, and extent of the tumor by relying on medical imaging techniques.In particular, diffusion tensor imaging (DTI) is one of the most common radiological methods in the detection of brain tumors [32].On T1-and T2-weighted images neoplastic and healthy tissue differ in intensity, thus making the core part of the tumor visible.
DTI techniques are non-invasive and they permit to assess the specific brain architecture.Modern imaging methods rely on evaluating the diffusion behaviour of water molecules within structured tissue.This reveals the orientation and anisotropy (meaning the extent to which fibers align together) of neuronal fibers making up white matter tracts.This information is particularly interesting for surgical and radiation therapy.
Several key features have been identified to be crucial for tumor progression, see e.g., [17].In this work we focus on the interaction of tumor cells with the underlying tissue and on proliferative effects.The model we propose here has a multiscale character, with the purpose of connecting the various levels on which the relevant biochemical processes triggering glioma spread take place.For modeling tumor growth we rely on the go-or-grow hypothesis, which states that cancer cells can either move or proliferate [15,21].We assume that cancer cells use neuronal fibre tracts as invasive pathways, which makes the individual brain structure interesting for patient specific modelling of glioma.Diffusion tensor imaging data are able to provide this information [34,35].Starting in Section 2 with a kinetic model on the mesoscale and an ODE on the subcellular level, a parabolic scaling yields an approximating reaction-diffusion-transport equation on the macroscale.Thereby, the precise form of the coefficients is determined by using DTI data.Section 3 is concerned with the numerical simulation of the macroscopic equation.Eventually, in Section 4 we comment on the performance of this modeling approach.

Multiscale modelling of glioma growth
The invasion of tumor cells is a highly complex process involving a plethora of phenomena on different spatial and temporal scales [17].Multiscale models can help to map this high degree of complexity into mathematical structures.Such mathematical settings involve two or several scales and accordingly more or less details of the processes they describe.Micro-macro models for tumor cell migration have been proposed and analysed e.g., in [28], which accounts for the effect of heat shock proteins (microlevel) on cancer invasion (on the macrolevel) and extends the model in [33], or more recently for the migration of cancer cells performing haptotaxis and chemotaxis (macroscale) under the influence of integrin binding to soluble and unsoluble ligands (microscale) [27,30] or for the acid-mediated tumor invasion with pH taxis [31].Different, but related approaches are hybrids between discrete and continuous modeling.These use cellular Potts or (lattice gas) cellular automata to provide lattice based simulations of individual cells involving more or less details on the subcellular level, from which information about the behavior on the population level can be extracted, see e.g., [8,18] and the references therein.More complex continuum models also involving mesoscale dynamics (individual cell-tissue and/or cell-cell interactions) were considered by Bellomo et al. [4,5,6] or -fitting in the same framework-in e.g., [12,22,23,26].In our model for glioma growth and spread we distinguish between three different scales.The microscale refers to processes happening on the subcellular level.In our setting these reduce to the binding of cell surface receptors to unsoluble ligands (proteins on the myelinated axon fibers) and employs an ordinary differential equation for the corresponding mass action kinetics.The mesoscale accounts for the behavior of individual cells and their interactions with their surroundings, in this case the (anisotropic) tissue network.The dynamics are characterized with the aid of kinetic transport equations for the densities of moving and resting cells (see Section 2.3).At this point, we rely on the above mentioned go or grow hypothesis.Finally, the macroscale refers to the population level on which the tumor as a whole can be observed.Unlike the previous two scales, it is not contained in the original model, but deduced from it.More precisely, it is our goal in this paper do derive an equation for the macroscopic density N (t, x) of tumor cells at time t at position x, which also contains the information of the descriptions of the microand mesoscales.

Diffusion tensor imaging
At present, medical imaging methods like computed tomography (CT) along with magnetic resonance imaging based methods like diffusion tensor imaging (DTI) play a central role in the diagnosis of brain tumors in general.Cancerous tissue and healthy tissue differ in intensity in the representation, and hence malignant tissue alterations can be observed.One significant shortcoming is that not the whole tumor is visible, but only its core part.The actual tumor margin cannot be assessed.This is a large drawback, since after resection of the malignant region visible on DTI or CT scans the tumor can start growing again due to the incomplete removal.However, in particular the diffusion tensor imaging (DTI) method is of high value for glioma prognosis as it provides information about the local tissue strcuture within the brain [35].DTI relies on the measurement of the diffusion behaviour of water molecules within structured tissue, relating it to the anisotopy of neural fibre tracts [9,34].The water molecules are assumed to diffuse faster along fibre structures than orthogonally to them, which allows to capture the mentioned anisotropy.Technically, DTI measures the aparrent diffusivity of water molecules per volume element (voxel).This can be characterized with a symmetric tensor, the water diffusion tensor Thereby, the diffusion tensor can be visualized e.g., with the aid of the fractional anisotropy index [7] or with tensor glyphs.While the former is a scalar value between zero and one, computed by using the eigenvalues of the apparent diffusion tensor1 , the latter convey tensor variables by mapping the tensor eigenvectors and eigenvalues to the orientation and shape of e.g., a cuboid or an ellipsoid.Glyphs commonly used in this context are ellipsoids and peanuts.For more details on this topic we refer to [12] and the references therein.Tractography aims to reconstruct the local tissue orientation out of DTI data [10,11].This is a quite complex process and still part of ongoing research.As we are primarily interested in one main direction of tissue fibres per voxel, we will use a simple mathematical expression to connect the water diffusion tensor and the local fibre distribution in Section 2.4.

Transport equations with resting phases
We aim at providing a multiscale model for glioma spread wich includes proliferation and interactions of the cells with the extracellular matrix.It has been proposed that the go-or-grow hypothesis is a good starting point in modeling proliferative effects on the cell scale [8,18].Here we start from the approach in [19] and propose a setting which is expanded by a model for the cell receptor dynamics on the microscale as in [12] and which also accounts for the interaction between tumor cells and tissue fibres.We aim at observing the density functions p(t, x, v, y) of moving cells at time t, position x ∈ R n , velocity v ∈ V ⊂ R n , and internal state y ∈ Y ⊂ R d and of resting (and hence proliferating) cells r(t, x, y).This leads to a system of equations of the form Here dv is the turning operator modeling the cell velocity innovations due to contact guidance, environmental cues etc.These influences are contained in the turning kernel K.As in [12], we choose K(x, v) = q(x,v) ω , where v is the normalized velocity, q(x, v) is the directional distribution of tissue fibers, and ω = V q(v)dv = s n−1 is a scaling constant.The function λ(y) denotes the turning rate of cells.The vector y stands for the internal state of a cell: its components can be, for instance, concentrations of proteins involved in some intracellular signaling network or -as will be the case throughout this work-the concentration of integrins bound to tissue fibres in the cell environment.Its dynamics is characterized by a system of ordinary differential equations where the right hand side is involving the volume fraction of tissue A(t, x) (including ECM and brain fibers) as input, see [12].For simplicity we account in this work only for binding of cell surface receptors to unsoluble components on the tissue, modeled by A. This assumption reduces the above ODE system to just one equation for the receptor binding dynamics, see ( 6) below.The function N in equations ( 2) and ( 3) denotes the total cell density and is given by and g(N ) and l(N ) are functions representing gain and loss due to cell death and proliferation.α and β are the rates with which cells stop and proliferate, respectively start moving after a resting (proliferating) phase.Thereby, the cells which exit the proliferative phase and start to move are doing this by interacting with the tissue. 2Both α and β depend on the position x.

Cell surface receptor dynamics
As in [12], we focus on the cell surface receptor dynamics of the cells.Transmembrane adhesion proteins bind to unsoluble proteins in the brain tissue.In this work we want to study the impact of this microscopic process on the macroscopic movement behaviour on the tissue scale.Starting from simple mass action kinetics, the receptor binding dynamics on the cell surface is described by the following ordinary differential equation: where R 0 denotes the total number of receptors on the cell (we assume it is conserved), whereas y(t) gives the density of integrins bound to ECM fibres.The constants k + and k − denote the reaction rates for the reversible binding of integrins to the tissue fibres.The steady state of equation ( 6) is given by y * = k + AR0 k + A+k − .Following [12] we introduce a new internal variable z := y * − y measuring deviations from the steady state.Next consider the path of a single cell starting in x 0 and moving with velocity v through a time-invariant density field A(x).Then with the notations x = x 0 + vt and it follows that for any t and hence z satisfies the equation This equation can be solved explicitly for z to yield where y 0 = y(0).Hence, z is bounded as long as ∇A it is and its sign depends on the current orientation of the cell w.r.t. the gradient of the fibre volume fraction.We choose the turning rate to be of the form λ(z) = λ 0 − λ 1 z ≥ 0, where λ 0 and λ 1 are some positive constants with ) and K r := k + /k − .Note that this is equivalent to considering λ(y) = λ 0 − λ 1 y * + λ 1 y, which becomes indeed larger when many receptors are bound to fibres.This choice corresponds to the one proposed in [13] for the turning rate of bacteria and -after a linearizationalso to that chosen in [16].Further, observe that z involves the spatial position x of the cell as a parameter and there exists ensures that λ(z) ≥ 0.

The mesoscopic transport equation and its scaling
With the above analysis for the dynamics on the microscale, we obtain now a system of evolution equations for the density p(t, x, v, z) of moving cells and the density r(t, x, z) of resting cells.The system reads Introduce the notations where Z ⊆ [y * − R 0 , y * ] is our new domain for the internal dynamics.In the following we assume the data to be compactly supported in the (x, v, z)-space, which allows to perform the subsequent calculations.Integration of the equations ( 9) and ( 10) with respect to z gives Multiply equations ( 9) and ( 10) by z and integrate again with resp. to z.As in [12] we neglect moments of higher order in z and obtain Now, we use a parabolic scaling t → ε 2 t, x → εx, then drop the hats for notational simplification.In addition, we scale the birth-death dynamics as in [19] by This scaling relies on the assumption that the time scale on which birth and death events occur is much slower than the random walk process.From equations ( 14), ( 15), ( 16) and ( 17) we obtain Now, we use Hilbert expansions for m, m z , w, w z , M, M z of the form and compare terms of equal order in ε.
ε 2 : It is our goal now to derive an evolution equation for the macroscopic cell density We start with the equations for the coefficients of ε 0 .First of all, we obtain from ( 27) that w z 0 = α β M z 0 .Integration of ( 26) with respect to v yields which implies M z 0 = 0, m z 0 = 0 and w z 0 = 0. Equation (25) gives Finally, we obtain from ( 24) Now let us analyze the equations for the coefficients of ε 1 .Starting with eq. ( 31) we obtain From equation ( 30) we obtain after integration with respect to v that M z 1 = 0 and hence also w z 1 = 0. Here, we used the assumption vq dv = 0 for undirected tissue.Then, again from (30) we obtain the following expression for m z 1 : Equation ( 28) leads to the expression Now use the properties of the operator L[λ 0 + α] defined on the weighted L 2 -space L 2 q (V ).Thereby, the weight function is q −1 (v) and L[λ 0 + α] is a compact Hilbert-Schmidt operator (see [20]) whose kernel is given by the linear space q , denoting the subspace of L 2 q (V ) spanned by q.The pseudoinverse L[λ 0 + α] −1 |<q> ⊥ of this operator is (see again [20] or [12]) is just the multiplication by the factor − 1 λ0+α .This leads to Using equations (32) and (33) we can derive an evolution equation for the macroscopic cell density (34).
From (33) we have Plugging this term into equation ( 32) we obtain (42) After integrating again with respect to v this yields Rearranging gives With the defined macroscopic quantity we obtain Now evaluate V vm 1 dv.Recall It follows Finally, we obtain with M 0 = β α+β N 0 and the evolution equation Notice that we still have to determine the form of q in order to obtain explicit expressions for the tumor diffusion tensor D T (x).

Determination of the coefficients
Now we are interested in the form of the diffusion tensor D T (x) in equation (52).To determine its specific form we have to choose a concrete form for the fibre distribution q.As shown in [12], a possible choice can be q The tensor D W (x) is the water diffusion tensor and can be obtained for each voxel x by DTI measurements.This formula allows for computing the fibre distribution function for each voxel.In [12] it has been shown that for the choice (53) for the distribution q the tensor (51) has the form The rates α and β are assumed not to vary drastically within the brain, and are hence constant in space.

The full macroscopic model
For the sake of clarity we collect the results of our calculations and formulate the whole macroscopic model again.We obtained from our scaling an approximating partial differential equation of advectiondiffusion-reaction type for the macroscopic cell density N 0 .We assume that the rates α and β do not vary in space and furthermore assume that the proliferation behaves logistic.Specifying g and l to be of the form with appropriate positive constants c g , c l ∈ R, we obtain the equation

Numerical simulations
We present 2D simulations of the resulting macroscopic advection-diffusion-reaction equation (55).The tumor diffusion tensor, the tumor drift velocity and the local growth rates are precalculated using MAT-LAB 3 .The simulations of the macroscopic evolution equation are implemented using the DUNE framework [2].The macroscopic quantities D T (x) and v(x) are spatially dependent, and we expect regions in space where the system is diffusion dominated and regions where it is drift dominated.This can be expressed by the spatially dependent Péclet number Pe = v(x)L D T (x) , where the norm in the numerator is the L 2 norm, while we take in the denominator the Frobenius norm of D T .L is a macroscopic characteristic length scale (i.e. a mesh width of 2 mm, in our present setting and simulations).For Pe 1 we have a diffusion dominated case, and classical methods for parabolic equations will apply, while for Pe 1 the equation is drift dominated, and hyperbolic numerical methods need to be used.Note that the equation is still parabolic, but the numerical method should be chosen based on its characteristics scales.The numerical scheme needs to be able to handle non-linear degenerated anisotropic parabolic equations and full tensors.Furthermore it should be locally mass conservative.For these reasons we decided to employ a first order discontinuous Galerkin (dG) scheme in space and an implicit Euler scheme in time.The non-linearities are handled using an outer Newton scheme for each time step.This method is overall first order accurate.

Spatial Discretization
For the spatial discretization we use a structured mesh M(Ω) = {E i }, which is a subset of the voxel mesh of the DTI data.Since our mesh is Cartesian with mesh-width h, the sets E i are simply the grid cells that belong to the brain tissue.We define the skeleton of M as the boundary between those grid cells, We use a symmetric interior penalty discontinuous Galerkin method (SIPG) [36] as implemented in DUNE [2].We denote with V h the dG trial-and test space, which is in our case the space of piecewise bilinear polynomials where Q 1 (E) denotes the set of bilinear functions on E. We test the above equation (55) with these ansatz functions, use a weak formulation, and express coupling conditions along the skeleton through jumps • and averages { • } in a similar notation as in [1].We use weighted averages, as proposed in [14].Omitting the computational details, which can be found for example in [1], the resulting semi-discrete problem reads: with a nonlinear reaction term R, the bilinear form a h , and a penalty term J h given by and where η > 0 denotes the penalty factor.For a good choice of η we refer to [14].
Note that we require the normal velocity to be continuous, i.e. v| ∂En n n = −v| ∂Em n m , where n n denotes the outer normal vector of grid cell E n .To ensure this, we use a Raviart-Thomas RT0 approximation of the velocity field.On an internal edge with two adjacent cells E n and E m we define x = x| ∂En n n +x| ∂Em n m and { x } = ω n x| ∂En + ω m x| ∂Em , with weights ω n , ω m and the unit outer normal vectors n n , n m .To be robust with respect to heterogeneous diffusion coefficients, we use the weights ω n =

Temporal Discretization
For the time discretization of N 0h we use an implicit Euler scheme and choose the time step width τ such that the discrete CFL (Courant-Friedrichs-Lewy) condition of 1 is satisfied.For each time step we obtain a non-linear system, due to the growth term.The modeling of the growth term requires that the tumor cell density is positive.The employed dG scheme is not monotone, small undershoots are possible and will occur in practice.Therefore we cannot directly compute R(N 0 ), but employ a clipping at zero.The non-differentiability of the clipped operator at zero leads to severe problems in the Newton method.Thus, we additionally regularize the clipping with a parameter , which yields the clipping operator For a time step k + 1, k ∈ [1, N ), and a discrete time step width of τ , the fully discrete problem now reads: Find N k+1 0h ∈ V h such that After linearization of R we can employ a standard Newton scheme.To solve the linearized system a conjugate gradient method with an ILU precoditioner with zero-fill-in is applied.Table 1 shows the parameter values used in the simulations.
Table 1: Choice of parameters used for simulations.

Simulation results
The choice of appropriate terms for the gain and loss functions g and l is still an issue.When assuming -as in the most existing works on the macroscale-that cellular proliferation can be described with a logistic growth, the choice to make is to let g(N 0 ) = r be constant for the proliferation rate r g ∈ R + and l(N 0 ) = r l N 0 for r l ∈ R + .Altogether this leads to the reaction term Simulations of this model have been perfomed on a horizontal brain slice with a binary inital tumor placed in the left hemisphere.Figures 1 and 2 show simulations for the pure diffusion case and for the full system including haptotactic drift.As expected, proliferation causes sharper gradients compared to the advection-diffusion model without proliferative effects.The corresponding tumor growth patterns evolve faster than the patterns without proliferation.However, for the pure diffusion case the patterns do not differ significantly from the models without proliferation and, again, significant anisotropic behavior caused merely by diffusion cannot be observed.
Adding the haptotactic drift terms leads to branched structures, as Figure 2 shows.Obviously, the tumor spreads predominatly in brain white matter; from the clinical point of view these finger-like patterns are more realistic than an isotropic spread.

Discussion
The DTI-based multiscale model for glioma invasion proposed and analyzed in this work is an extension of the one introduced in [12].The novelty is here that we account for proliferation of tumor cells.Thereby, we made use of the go and grow hypothesis and identified the resting cells with those proliferating.The choice of gain and loss functions g(N ) and l(N ) due to cell proliferation and death, as well as of the stopping (growth) rates can be done in a less restrictive way, leading -after the scaling-on the macroscale to rather general proliferation terms, of which the logistic growth (which is commonly chosen when writing some partial differential equations for tumor invasion directly on the population level) is merely a particular case.The numerical simulations of the macroscopic equation show a highly anisotropic spread of the tumor cells, according to the underlying brain structure.As before in [12] this was achieved by the multiscality of our model (in particular the inclusion of subcellular receptor binding processes and their involvement in the cell turning rates), which led to the dominant haptotactic term in equation (55).When simulating only the pure diffusion model the corresponding invasive behavior is not as known from practice, where branched, fingering patterns like those in Figure 2 are observed.When compared with our previous multiscale model without proliferation in [12] the present setting predicts -as expected-a faster invasion, with larger cell clusters.The nontrivial issue of modeling cell proliferation in a way which is less artificial than just choosing it to be a logistic growth on the macrolevel was handled here by accounting for the resting (and hence dividing) cells.This, however, does not seem to be the only way to characterize it.Bellomo et al. [4,5,6]  proposed to model proliferation on the mesolevel by letting the cells with different internal states4 interact with each other, which is arguable, as cancer cells can divide without needing to stay in touch with their neighbours.The modeling of cell proliferation on the mesolevel and the deduction of the appropriate mesoscopic equation are thus an aim to follow.
T n +D T m )nn and ω m = n t m D T m nm n t m (D T n +D T m )nm respectively, where D T n denotes the tumor diffusion tensor computed in E n .Note that due to n n = −n m the relation ω n + ω m = 1 holds.The advective term is stabilized using an upwind formulation, with the upwind reconstruction N 0 ↑ h of the cell density.