Estimating the extent of glioblastoma invasion
NNoname manuscript No. (will be inserted by the editor)
Estimating the extent of glioblastoma invasion
Approximate stationalisation of anisotropic advection-diffusion-reactionequations in the context of glioblastoma invasion
Christian Engwer · Michael Wenske
Received: date / Accepted: date
Abstract
Glioblastoma Multiforme is a malignant brain tumor with poor prognosis.There have been numerous attempts to model the invasion of tumorous glioma cellsvia partial differential equations in the form of advection-diffusion-reaction equa-tions. The patient-wise parametrisation of these models, and their validation via ex-perimental data has been found to be difficult, as time sequence measurements aregenerally missing. Also the clinical interest lies in the actual (invisible) tumor extentfor a particular MRI/DTI scan and not in a predictive estimate. Therefore we proposea stationalised approach to estimate the extent of glioblastoma (GBM) invasion atthe time of a given MRI/DTI scan. The underlying dynamics can be derived from aninstationary GBM model, falling into the wide class of advection-diffusion-reactionequations. The stationalisation is introduced via an analytical solution of the Fisher-KPP equation, the simplest model in the considered model class. We investigate theapplicability in 1D and 2D, in the presence of inhomogeneous diffusion coefficientsand on a real 3D DTI-dataset.
Keywords glioblastoma modelling, stationalisation, reaction-diffusion
Mathematics Subject Classification (2010) · · · · Treatment of glioblastoma multiforme (GBM) turmors usually consists of a combina-tion of tumor resection (operation), radio- and chemotherapy (Sathornsumetee et al.,
C. EngwerInstitut f¨ur Numerische und Angewandte Mathematik, WWU M¨unster, M¨unster, GermanyE-mail: [email protected]. Wenske (cid:0) (Corresponding author)Institut f¨ur Numerische und Angewandte Mathematik, WWU M¨unster, M¨unster, GermanyE-mail: m [email protected] a r X i v : . [ q - b i o . CB ] J a n Christian Engwer, Michael Wenske (Swanson et al., 2008; Patel and Hathout, 2017), sothat tissues are segmented as healty, although they still contain a significant numberof tumor cells.In clinical practice an average extent of this invisible infiltration of 2 cm normalto the visible tumor is assumed. Using mathematical modelling the aim is to estimatethe extent of resection and radiotherapy to be applied outside of the tumorous regionsvisible on the medical images.
Modelling:
Many efforts have been made to mathematically model the behaviour ofthe tumorous glia cells within the brain. The mathematical approaches are numerousand include a wide range of effects possibly influencing the behaviour of the glioblas-toma (GBM) cells. One prominent mathematical approach is to describe the prolifer-ation and the movement of the tumor by macroscopic partial differential equations.Most of these models take the form of diffusion-reaction-advection equations.Harpold et. al wrote a review article in 2007 about the evolution of mathemat-ical modelling of GBM (Harpold et al., 2007). The modelling started from simplereaction-diffusion equations with exponential growth ( ∂c∂t = ∇ ( D ∇ c ) + ρc ) (Tracquiet al., 1995). From here, the extensions from homogeneous diffusion to distinguisheddiffusivities in grey- and white matter were performed by (Swanson et al., 2000).With the advent of diffusion tensor imaging (DTI) it was also possible to make useof directional information of water diffusion within the brain. Jbabdi et al. (2005)estimated the tumor diffusion from the water diffusion matrices available from DTIand simulated the tumor invasion making full use of the medical imaging informa-tion. It was also possible to link the predicted spreading velocity of the Fisher model( v ≈ √ Dρ ) to experimental data for low-grade gliomas. They effectively found thatthe increase in radius was linear. Thereby, giving further indication that the tumorgrowth may be modelled by reaction-diffusion equations (Mandonnet et al., 2003).This may indicate that the rate of advance can be estimated with the knowledge ofthe diffusive properties of the surrounding tissue, and the reproductive rate of the tu-morous glia cells. There are also quite rigorous approaches to link knowledge on mi-croscopic cell behaviour to macroscopic partial differential equations via upscaling.The mathematical framework was given by (Othmer and Hillen, 2000, 2002; Hillenand Painter, 2013). The underlying idea is to identify the most important influenceson the cell’s behaviour on the microscopic and mesoscopic scale and mathematicallyderive the influence on the macroscopic population dynamic. With this powerful ap-proach it has been possible to microscopically motivate haptotaxis, chemotaxis andproliferation (Engwer et al., 2016, 2015; Hunt and Surulescu, 2017; Jan Kelkel, 2011;Corbin et al., 2018).Woodward et al. additionally modelled the effect of resection on the tumor inva-sion, finding some correlation with real survival times of patients (Woodward et al., stimating the extent of glioblastoma invasion 3 Parameters:
There have been approaches to assess the growth characteristics and pa-rameters in in-vitro experiments, e.g. (Oraiopoulou et al., 2018). It is an open ques-tion whether the information gathered in these experiments is transferable to in-vivosituations. Caragher et al. investigated treatments using novel 3D cell culturing meth-ods in the context of GBM therapy developement (Caragher et al., 2019). Even ifin-vitro experiments may prove essential in understanding the involved biochemistryand qualitative effects, their use for the validation problem is limited. In order to im-prove our ability to accurately describe GBM invasion the interlinked problems ofparametrisation and availability of data have to be addressed. The reason for this isthat the large number of free parameters can not be met with experimental data toestimate them with reasonable accuracy.
Validation:
Although there are a number of mathematical descriptions available, ithas proven difficult to derive their parametrisation from medical data or experiments.In medical practice the diagnosis of GBM is often rapidly followed by a combinationof tumor resection, radio- and chemotherapy, thereby severely altering the growthcharacteristics of the tumor (Sathornsumetee et al., 2007). One problem is, that inorder to validate any given forward GBM-model given in the form of advection-diffusion-reaction equations, one would need a time-series of DTI/MR/CT scans ofthe patient, where the tumor growth had been left without cessation. That setup wouldallow for direct comparison of the in-silico experiments and the medical images ofthe progressing tumor. One also preferably had these datasets from a large number ofrepresentative patients. In order to retrieve a dataset like that, one needed to deprive ahigh number of patients of live prolonging treatment while undergoing regular med-ical scans. The ethical impossibility is obvious. It is also questionable whether anysuch datasets will be available in the near future. Even in more accessible subjects
Fig. 1
Schematic time line of ideal datasets for the validation of GBM forward models. Christian Engwer, Michael Wenske
Fig. 2
Schematic time line of a realistic sequence of measurements. The initial DTI scan after diagnosisis rapidly followed by a combination of gross tumor resection (OP), radio- and chemotherapy (RT, CT).During treatment, there may be follow-up MRI’s . like rodents, it proved difficult to determine a good parametrisation for numericalmodels. Rutter et al. studied tumor growth in five mice which were injected with tu-morous glia cells under controlled conditions. Even with a reportedly careful experi-mental setup, the resulting tumor sizes varied significantly and fitting parameters of asimple Fisher-KPP tumor model proved difficult (Rutter et al., 2017). Given the goalof improving the ability to accurately describe the tumor invasion in real patients, theproblem of validation/falsification has to be addressed.
Contribution:
We present a stationalisation approach, that opens the possibility toestimate the invisible tumor extent at measurement time, building upon existing in-stationary tumor growth models. It may also alleviate part of the parametrisationproblem by only being sensitive to the relative strength of physical effects and noton their absolute quantitative parametrisation. We will first state the class of consid-ered partial differential equations in section 2.1. In section 2.2 we present a method-ology to find a stationalised analytical expression for the 1D Fisher-KPP equation,determining the front shape of a travelling wave solution. We also explain how theanalytically derived stationalisation term may be used to approximate the tumor den-sity. After numerical verification of the derived gradient distribution, we investigatewhether this stationalisation approach is viable in the presence inhomogeneous mate-rial properties in section 4. Finally the advantages and shortcomings of the proposedprocedure will be critically discussed in 5. u : Ω × T (cid:55)→ R , withthe d-dimensional domain Ω ⊂ R d and a time range T = [ t , t e ] . The solution u describes the volume percentage of cancerous cells and therefore ≤ u ( x ) ≤ .The dynamic of the density profile is given in the form of a macroscopic partialdifferential equation. We consider the time-dependent fully anisotropic advection-diffusion-reaction equation in the following form: ∂u∂t − ∇ ( D t ( x ) ∇ u ) − ∇ (( ∇ · D t ( x )) u ) = ρu (1 − u ) in Ω × T, (1a) stimating the extent of glioblastoma invasion 5 with boundary and initial conditions u ( x ,
0) = g ( x ) in Ω, (1b) ∇ u ( x , t ) = 0 on ∂Ω × T, (1c)with the diffusion parameter D t ( x ) being a symmetric positive definite matrix D t ( x ) ∈ R d × d , i.e. the diffusivity can be inhomogeneous and anisotropic. Equation (1) is aprototype model for the invasion of GBM in the sense that many models differ fromit merely in the reconstruction of the tumor diffusion matrix D t ( x ) from DTI data, orthe addition of other chemotaxis terms (Hunt, 2018).2.2 Stationalisation of the Fisher-KPP equationIn one dimension and for isotropic diffusive properties with ( D t = 1 , ρ = 1 ), equa-tion (1) degenerates to the classical Fisher-KPP equation (Kolmogoroff et al., 1988;Fischer, 1937); ∂u ( x (cid:48) ) ∂t = ∆u ( x (cid:48) ) + u ( x (cid:48) )(1 − u ( x (cid:48) )) , (2)with x (cid:48) ∈ R as a spatial coordinate. The first term may be interpreted as the passivediffusive spread of the cells due to a random walk, the rightmost term as a logisticproliferation term, as often encountered in biological contexts. For the initial condi-tions: lim x (cid:48) →−∞ u (cid:48) ( x ) = 1 , lim x (cid:48) →∞ u (cid:48) ( x ) = 0 , the equation allows travelling wave solutions (Kolmogoroff et al., 1988). These solu-tions converge over asymptotically long times cales to a wave-front which is constantin shape, moving laterally with a globally constant velocity v ∈ R . With that infor-mation, equation (2) can be stated in equivalent form by introducing a co-movingframe x = x (cid:48) − vt . ∆u ( x ) + u ( x )(1 − u ( x )) − v · ∇ u ( x ) (cid:124) (cid:123)(cid:122) (cid:125) stationalisation (3)The advective term − v · ∇ u ( x ) results from the coordinate transformation into thecomoving frame and can be understood as simply counteracting the lateral movementof the propagating front (Ablowitz and Zeppetella, 1979; Kolmogoroff et al., 1988).We will refer to the advective term and its approximations as the stationalisationterm . Any laterally shifted profile u ( x − c ) , with a constant c ∈ R is equally asolution to (3). The minimum limit speed of the travelling wave solutions has beenstated to be v = 2 (cid:112) D t ρ, but Kolomogorov also states that there are solutions for speeds higher than v . Thisfact is due to the reaction term acting independently of the local environment. For thelimit of diminishing diffusion, the reaction term will induce very high travelling wavespeeds for initial conditions which are close to horizontal. Contrarily, the wave speed Christian Engwer, Michael Wenske will be determined to a higher degree by the diffusive flux for initial conditions thatfall from 1 to 0 in a very localized region (Kolmogoroff et al., 1988). In the context ofglioma invasion only the very rapid spatial decay of tumor density is relevant. Eventhough the sigmoidal shape of the advancing front profile is formed quite rapidly, itmay take asymtotically long to reach a state where the shape of the advancing frontis no longer changing. For long time scales and no boundary conditions, a stationarywave form can be given in analytical form for the special wave-speed of v = √ u ( x ) = 1(1 − r exp( x √ )) , (4)where r = − (Ablowitz and Zeppetella, 1979). We will call that profile the equili-brated wave-form. In the fully equlibrated case, we can state an analytical expressionfor ∇ u ( x ) . Since the analytical solution is invertible in the relevant range u ∈ [0 , ,we may also express the gradient in terms of u . Assuming an analytical expressionfor p ( u ) = v ∇ u ( u ) , we may express the comoving 1D Fisher-KKP equation (3) inequivalent form: ∆u + u (1 − u ) − p ( u ) (cid:124) (cid:123)(cid:122) (cid:125) stationalisation . The gradient of the analytical solution, can be given as ∇ u = − (cid:113) exp( x √ )(1 + exp( x √ )) . (5)The inverse of the analytical solution is given by x ( u ) = √ (cid:16) ± √ u − (cid:17) . (6)and represents the mapping from a given amplitude u to the corresponding posi-tion x within the wave form profile. Substituting the inverse (6) into the gradientexpression(5) yields an analytical term for ∇ u ( u ) and thereby a closed form for thestationalisation term p ( u ) = | v | (cid:114)
23 (1 − √ u ) u. Note that the strict equivalence between p ( u ) and − v ∇ u only holds for the limitcase, in which the analytical solution is available. Also, the expression only holdsfor the case without boundary conditions, i.e. Ω ≡ R . Although this formulationis only strictly equivalent in case of a completely equilibrated wave form, we willuse the penalty term p ( u ) as an approximation of v ∇ u ( u ) leading to a stationalisedwave-pinning type model: ∆u + u (1 − u ) (cid:124) (cid:123)(cid:122) (cid:125) tumor model prototype − | v | (cid:114)
23 (1 − √ u ) u (cid:124) (cid:123)(cid:122) (cid:125) stationalisation in Ω, (7a) ∇ u ( x ) = 0 at ∂Ω, (7b) u ( x ) = c at ∂Ω i . (7c) stimating the extent of glioblastoma invasion 7 Within the encompassing domain Ω we define a smaller enclosed domain Ω i ⊂ Ω representing the tumor visible on the medical images. The additional internal bound-ary condition on ∂Ω i is used to localize the stationalised solutions.The Fisher equation is one example out of a family of KPP-type equations whichcombine a diffusive term with a nonlinear reaction term f ( u ) . The reaction term is of-ten chosen in a manner so that it dynamically connects two fixed points of amplitude: f (0) = 0 , f (1) = 0 , f (0 < u < > . Although an exact analytical description ofthe stable wave fronts proves difficult, the rough characteristic of propagating frontsfound in nature (combustion, bacteria growth etc.) is often similar to a sigmoid func-tion. The gradient distribution of any sigmoid-like travelling wave front will have asimilar shape like p ( u ) . For any sigmoidal wave-form |∇ u | will be zero for u = 0 and u = 1 and of higher amplitude for < u < .The underlying idea of the stationalisation procedure is, that the gradient distri-bution, and therefore the penalty term necessary to calculate the travelling wave formstationally, may possibly be approximated to a good degree of accuracy and evenfor those cases where the analytical solution is not at hand. Considering the generalmodel in the form of (1), we define the corresponding stationalised problem as: Find u s such that ∇ ( D t ( x ) ∇ u s ) + ∇ (( ∇ · D t ( x )) u s )+ ρu s (1 − u s ) − | v | (cid:114)
23 (1 − √ u s ) u s , in Ω, (8a) ∇ u ( x ) = 0 at ∂Ω, (8b) u ( x ) = c at ∂Ω i . (8c)The above problem closely resembles the general model (1a) augmented with thepenalty term p ( u ) and internal boundary conditions. The matrix D t ( x ) is a recon-struction of the tumor diffusivity from the DTI datasets, ρ is a growth parameter and v the penalty parameter.It is known that in higher spatial dimensions, the Fisher-KPP equation has relatedsigmoid-like travelling wave solutions. There are more additional stable wave frontpatterns in higher dimensions. One of them describes a v-shaped waveform propagat-ing through the two-dimensional medium, which can be interpreted as two straightwave fronts collapsing into each other at a certain angle. Observed in the directionof a bisecting line, this combined wave front is indeed stationary at certain speeds.There also exist spatially oscillating front shapes, but these profiles are only possiblefor the extension of u ( x ) out of the relevant range u ∈ [0 , (Brazhnik and Tyson,1999). We focus on cases where the wave propagation occurs as a radial expansionfrom a centred mass. In the direction of the propagation, equation (4) provides goodestimates on the wave fronts profile. In this section we discuss the methods used in our numerical experiments to mea-sure the modelling error of the proposed stationalisation. In section 4 we will use
Christian Engwer, Michael Wenske the methods to investigate the validity of our approach in a series of problems withgrowing complexity.3.1 Numerical schemeWe follow the method of lines approach to split temporal and spatial operators. Thetemporal discretisation is an implicit-euler scheme, which is unconditinally stable.For the spatial discretisation we use for simplicity a standard finite element dis-cretisation on cubic grids with multi-linear trial- and test-functions.The matrix divergences are pre-calculated by a first order finite difference approx-imation within each grid cell. The inhomogeneous diffusion matrices at the quadra-ture point are evaluated by nearest neighbour interpolation.Positivity of our solution might be violated in finite precision calculations. Wher-ever the numerical iteration leaves the physically sensible range of u ∈ [0 , , wedisregard the reaction term of the given forward- or stationalized model and insteademploy the following artificial numerical penalty term n ( u ) = (cid:40) − ρu if u < ρ (1 − u ) if u > . (9)This is done, because a logistic reaction term ∝ ρu (1 − u ) may otherwise amplifynumerical fluctuations which produce a slightly negative amplitude in u .The stationalisation procedure should be largely independent of the chosen nu-merical discretisation. The largest source of error does not lie in the numerical treat-ment, but in the approximations made in the parametrisation and in the estimation ofthe internal Dirichlet constraints.3.2 ImplementationThe implementation was realized within DUNE software framework (Bastian et al.,2008b,a; Blatt et al., 2016). The finite element discretisation was implemented withindune-pdelab (Bastian et al., 2010). The non-linear system is solved with a classicalNewton-Krylov method, using linear search. The linear systems are solved with anAMG- preconditioned BiCGSTAB solver, using the dune-istl module (van der Vorst,1992; Blatt and Bastian, 2007). The release version for all DUNE modules was 2.6.For the realistic head model the diffusion matrices were reconstructed by thecamino software package (Cook et al., 2006).3.3 Error measureIn section 4.2 we will investigate the impact of the stationalisation error on the ob-served tumor front. In this course we will compare the reconstructed tumor front of afully instationary simulation with the reconstructed tumor front using the stationali-sation approach. stimating the extent of glioblastoma invasion 9 Fig. 3
Schematic representation of A ⊕ B : grey regions Given a reference solution u a , an approximation u b and a threshold value θ wedefine two domains A and B as A = { x | u a ( x ) ≥ θ } , B = { x | u b ( x ) ≥ θ } . The medically relevant information is the spatial discrepancy between two level-sets( ∂A, ∂B ) of these density profiles. An absolute measure for this error is the symmet-ric difference A ⊕ B , as depicted in Figure 3. It measures those sub-volumes whichare either included A but not in B , or vice versa. That way, both over- and underes-timations of the approximation u b are represented. The most expressive informationin the medical context might be the average distance between the two level-sets. Wetherefore introduce the characteristic level-set distance. Definition 1 (Characteristic level-set distance)
For a given level-set value θ , wedefine the characteristic level-set distance between ∂A and ∂B as L B := | A ⊕ B || ∂A | . (10)It quantifies the average distance between theAssuming a spherical reference geometry for A we can simplify this expressionand avoid evaluating | ∂A | . Given the radius r A (or an approximate average radius of A ), L B simplifies to L B = | A ⊕ B | , L B = | A ⊕ B | πr A and L B = | A ⊕ B | πr A . u ( x c ) = c within the domain and Neumann bound-aries at the border, we may find two solutions mirrored around the constrained pointwithin the profile. Since we want to use an internally constrained region given by seg-mentation information we may expect two possible solutions to the formulation (3).The first solution corresponds to the stationary approximation to the travelling wavesolution of the forward model, which moves outwards from the tumor center, while satisfying the Dirichlet constraint. This solution will have relatively low mass outsidethe constrained region. The second possible solution corresponds to a travelling wavemoving into the constrained region, and will have a high mass on the outside. Thetwo solution types are sketched in Figure 4. When considering the effective reaction . . . . u [ a r b ] x [arb]HMSLMSthreshold . . . . − .
02 0 0 . u [ a r b ] effective reaction [arb]eff. reactionthreshold Fig. 4 Left:
Schematic of high- and low mass solutions (HMS, LMS), both fulfilling the internal Dirichletconstraints (black dots).
Right:
Plot of effective reaction term consisting of the logistic growth and thepenalty term. The penalizing regime is indicated in grey. The visibility threshold is within the penalizingregime. term consisting of the logistic growth combined with the penalty term, we find that ithas a penalizing regime for < u < and a growth regime for < u < (Fig. 4).The diffusive process transports mass from high amplitudes to lower amplitudes. Thecombined reaction term counteracts this process. The visibility threshold, and there-fore the constraints, lie within the penalizing regime. In order to select the low masssolution branch, we use initial guesses u s ( x ) (cid:28) which are well within the purelypenalizing regime, so that the newton iteration converges to the low mass solutionreliably. We present the results of the numerical validation studies of the stationalisation pro-cedure in 1D and 2D. We first compare the gradient distributions of forward simu-lations with the analytical expression found in equation (2.2) in section 2. After thatwe compare forward simulations with their stationary approximations and assess theimpact of our modelling error on the reconstructed tumor front. We present an easilyreproducible 2D example in 4.2.2. Finally we present the applicability to a realistic3D DTI dataset in 4.3.4.1 Investigation of modelling errorWe will numerically investigate the effect of inhomogeneous diffusion on the gradientdistribution, and by this the applicability of the chosen stationalisation term. As theanalytic formulation of ∇ u is only valid in the homgeneous 1D case, we expect amodelling error, which we will assess in different test cases. stimating the extent of glioblastoma invasion 11 . . . . . . | ∇ u ( u ) | u1D gradient distribution |∇ u ( u ) | inhomogeneous with β = 0 . , t (cid:29) inhomogeneous with β = 0 . , t (cid:29) homogeneous, t (cid:29) ∇ u ( u ) analytic Fig. 5
Scatter plot of 1D gradient distribution: magnitude of gradient of 1D wave-front profile. The numer-ical solution approaches the analytic expression only at asymptotically large time-scales. The underlyingcharacteristic is not destroyed by the introduced inhomogeneities.
We introduce a test cases with randomly perturbed diffusion coefficients. D β = d ( δ ) d , (11)with δ being a uniformly distributed random value with a spread of β around anaverage of 1. Here, d is the spatial dimension. The exponent of δ is chosen to allowcomparisons between the isotropic homogeneous case and the randomized case, byassuring that the average of the determinants of the diffusive medium are close to 1.0for every realisation of the random field: | ¯ D β | = ¯( δ ) ≈ . (12)We only compare the homogeneous and inhomogeneous case with equal grid resolu-tion. We evaluate the random inhomogeneous diffusion on the dual grid and assumeit to be piecewise constant therein. The diffusivity within one dual grid cell is sta-tistically independent from any neighbour. The effect of statistical scattering of thediffusive properties on the macroscopic front speed is non-trivial. Since the globalfront speed appears as a linear factor in the analytic derivation of the stationalisa-tion term, we may not expect perfect convergence to the analytically derived gradientdistribution. The results may however illustrate the effect of realistic datasets. We consider a one dimensional forward simulation of equation (1) starting from aGaussian initial condition in the center of a one dimensional domain x ∈ [0 , .We first simulate the homogeneous case with D t = 1 where we expect perfect con-vergence of the gradient distribution to the analytical expression. In the homogeneouscase there are no advective terms active, as ∇· D t is zero. Upon start of the simulation,two travelling waves form and move away from the center of the domain. After a short initial phase, the wave-fronts in the homogeneous medium asymp-totically approach the front shape of the analytical solution (4) and its symmetriccounterpart. The front speeds within an inhomogeneous material are slightly per-turbed. Similarly the relation ∇ u ( u ) is disturbed.In Figure 5 we investigate how the gradient distributions approach the analyticalexpression for the homogeneous and inhomogeneous case and how strongly it differsfrom the analytical expression in the case of heterogeous coefficients. We consider a square domain ( L x = L y = 200 ) in 2D with a Gaussian initial con-dition in the center. After the start of the simulation, the wave-front circularly prop-agates outward from the central point. Slow convergence of the gradient distributionto the analytic expression (2.2) can be replicated in the 2D homogeneous case. Con- . . . . . . | ∇ u ( u ) | u2D gradient distribution |∇ u ( u ) | inhom. β = 0 . , t (cid:29) inhom. β = 0 . , t (cid:29) homogeneous, t (cid:29) ∇ u ( u ) analytic Fig. 6
Gradient distribution of a 2D simulation with different spreads in diffusive properties ( β =0 . , . , . ). trary to the one dimensional case, there is an effect of curvature in higher dimensionsslowing the convergence towards the 1D gradient distribution. However in the limitcase the wave propagation still reduces to a 1D dynamic in the propagation direction(J.D.Murray, 2007). It is obvious that the introduction of random material propertiesbreaks the strict applicability of the stationalisation term. However Figure 6 suggestthat also in 2D the underlying polynomial relation between the wave-fronts ampli-tude and its gradient is merely perturbed by the material properties. Compared tothe usual parameter uncertainties, we consider this modelling error small and thusa stationalised solution, making use of the analytical gradient distribution, shouldstill provide reasonable estimates on the density profile. Since the global propagationspeed v appears as a linear factor, any process that alters the propagation speed awayfrom the analytical value should have a linear effect on the gradient distribution. stimating the extent of glioblastoma invasion 13 ∇ u . The numerical observation in section 4.1 suggests that the average behavior isstill well described by our analytic reformulation 7a. In the following we will inves-tigate the impact of this stationalisation error on the actual tumor front. We will usetest cases with growing complexity. . . . .
81 0 10 20 30 40 50 60 70 80 90 10000 . . . .
81 0 10 20 30 40 50 60 70 80 90 100 t u m o r d e n s i t y Homogeneous case: D t = 1 t u m o r d e n s i t y xInhomogeneous case: β = 0 . . . . . .
25 0 2 4 6 8 10 12 1400 . . . . .
25 0 2 4 6 8 10 12 14
Zoom of homogeneous caseu(x,t= 0)u(x,t=10)u(x,t=20)u s ( x ) threshold xZoom of inhomogeneous case Fig. 7
Direct comparison of 1D forward simulations and their stationalizations. The black dots indicatethe begin of the internal constraint given by the simulated imaging threshold ( u = 0 . ). The domain wasdiscretised into 1000 equidistant elements. Top left.:
Initial condition, and states of the forward simulation(homogeneous).
Bottom left.:
Initial condition, and states of the forward simulation (inhomogeneous).
Top right.:
Zoom of the forward solution and the corresponding stationalization (homogeneous).
Bottomright.:
Zoom of the forward solution and the corresponding stationalization (inhomogeneous). Both theforward and the stationalised solution show slight deviations from a smooth decay in the inhomogeneouscase. Note that the front speeds are not perfectly identical.
We try to mimic the situation observed in the medical application and presentthe procedure of estimating the tumor extent at the time of diagnosis. To generateartificial datasets in controlled scenarios, we simulate the carciogenesis by first as-suming a small Gaussian initial condition. Secondly, we propagate the density profilefor a certain time simulating the uninterrupted tumor growth. Finally, at the time ofdiagnosis, we use a level-set of u = 0 . to represent the thresholded medical imagingmodalities. Other choices of threshold value are possible. We then use the thresholdedvolume as an internal Dirichlet constraint and solve the stationary problem (8a). Inthis numerical setup both the full forward density profile u and the stationary profile u s are known and compared. In any real world situation only the thresholded imageinformation would be accessible. We first show results for a simple 1D case with D t = 1 and with inhomogeneouscoefficients ( β = 0 . ). By introducing inhomogeneous diffusive coefficients, the ad-vective drift term will produce small contributions to the equation. We chose thepenalty parameter v = √ and ρ = 1 . The one dimensional setup is not very real-istic, but practical to illustrate the procedure. Figure 7 presents a direct comparisonof forward simulations with their stationalized counterparts. The snapshots at dif-ferent times of the forward solution suggest, that the rough form of the advancingfront is formed rapidly. If the forward solution converges rapidly towards the form ofthe analytical solution, with only diminishingly small corrections to the wave frontsshape at larger times, then correspondingly the approximation to the gradient statisticwill perform well even early in the simulation. The direct comparisons show that thestationalization produces a reasonable approximation to the full forward solutions.The inhomogeneous coefficients induce small deviations from a smooth front shape,which are present in both the forward simulations profile as well as in the stational-ized solution. It is to be expected that wherever the internal constraint results fromthresholding of an underlying smooth distribution, the stationalisation will performwell since the real density distribution is close to the equilibrated wave-form. In order to assess the viability of the stationalisation for more realistic tumor modelswe now move to two dimensions with inhomogeneous coefficients (eq. (8a)). We setup an inhomogeneous but isotropic field for the diffusion matrix by scaling the unitmatrix according to its x position D t ( x ) = (1 . (cid:16) πL x x (cid:17) . . (13)This effectively separates the domain in three distinct regions, with the left- and right-most third of the domain having higher diffusivity and the middle strip having re-duced diffusivity. The changes in diffusivity may represent grey and white matterregions in a primitive way. In this example there is more long-range deviation of thediffusive properties than in the 1D examples in section 4.2.1. We again chose thepenalty parameter v = √ and ρ = 1 . In a medical situation the task is to estimatethe region and intensity of radiotherapy to be applied and the area of resection fromonly the thresholded information visible at the time of diagnosis.Figure 8 compares predictions of the forward simulation of equation (1) with thestationary solutions of equation (8a). The stationalisation greatly benefits from theinternally constrained region. Since most of the tumor mass is above the visibilitythreshold, the stationalization only has to provide the estimation on the surroundingregion. Similar to the 1D constrained situation, the stationalisation captures the profilequite well, but slightly overestimates the invasion extent in low density regions. stimating the extent of glioblastoma invasion 15 Fig. 8
2D inhomogeneous isotropic testcase.
Top left: level-sets on the forward solution at t=20.
Topright: diffusivity as given in (13) and horizontal cut through the Gaussian initial condition at x = (50 , . Bottom left: level-sets on the solution of the stationary problem.
Bottom right:
Error field indicatingregions of over- and underestimation. software project (Cook et al., 2006). We relate thetumor diffusion to the water diffusion by a simple scalar factor: D t = α D w . (14)More advanced reconstructions are possible, but not central to this example. A vari-ety of different tumor diffusion models has been proposed in the literature, see e.g.(Hunt, 2018; Painter and Hillen, 2013; Conte et al., 2020). We use the forward model(1) and the corresponding stationary problem (8a), with the parameters in table 1 toscale the equations to a realistic range. Here α is a dimensionless parameter, v the http://camino.cs.ucl.ac.uk/index.php?n=Tutorials.DTI parameter value α ρ v Table 1
Parameters used to scale the terms in (1) and (8a) to realistic ranges. The penalty parameter v hasthe unit of /s since we approximated the gradient within the advective term with the analytical expression(2.2). Fig. 9
Horizontal slice of the 3D results and dataset. The black dot indicates the position of the smallGaussian initial condition at x = (0 . m, . m, . m ) . Top left: level-sets (0.2, 1e-3, 1e-4) on u ( x ) after 90 days. Top right:
Fractional Anisotropy of the reconstructed tumor diffusion matrix D t ( x ) fromthe camino dataset (Cook et al., 2006). Bottom left:
Identical level-sets on u s ( x ) . Bottom right:
Regionsof over- and underestimation by the stationary solution. penalty parameter, and ρ is a growthrate in /s . We again follow the procedure de-scribed in section 4.2 and start a forward simulation from a small Gaussian at t = 0 until t e = 90 d , and use a level-set ( u = 0 . ) as the constrained region for the station-alisation. Figure 9 shows the direct comparison of the forward simulation and thestationalization. All local extentions or reductions induced by the local increase or stimating the extent of glioblastoma invasion 17 L c [ m ] Level-set value determining the symmetric difference regionCharacteristic length error L c at t = 90 d . . . .
004 0 0 . .
001 0 . . Fig. 10
Characteristic length errors L B for given level-set values for the 3D example camino dataset. Thegrey region indicates the errors between the level-sets used in Figure 9. decrease of the underlying diffusivity are present in both the forward and the station-alized solution. For this particular example, we measured the characteristic levelsetdistance L B as given in (10) for a series of small levelset values. Figure 10 indicatesthat the distance error is between − mm for the level-sets chosen in our numeri-cal example. These levelsets were chosen to replicate the current treatment radius ofabout 2cm. We presented a stationalisation approach for the estimation of the glioblastoma inva-sion extent. The stationalisation approach partially addresses the problems of para-metrization and data availability. The stationary simulations do not depend on the complete knowledge of the initial condition to produce reasonable tumor invasion es-timates. The thresholded information provided by the medical images might be fullyutilized with the limited information it provides. The stationalisation only requiresdatasets from one point in time, i.e. one DTI scan and a medical segmentation ofthe tumorous region. This is the data that is routinely gathered in medical practice,as it is used for planning the radiation therapy and resection. Stationary simulations,as presented here, may provide an additional tool in this regard without altering theimaging practices.The problem of quantifying model parameters is also partly alleviated by the factthat in the stationary formulation the solution depends merely on the strength of theparameters with respect to each other, and not their absolute values. It may be easierto experimentally determine ratio values between the three important parameterizedterms: D t , v and ρ than the full set of parameters needed for a forward simulation.Forward models which fit into the form of equation (1) may only produce reliable re-sults if the correct initial condition g ( x ) is known and the parameters are determinedto a reasonable degree of accuracy. In the medical setting the only information athand is the medical imaging at time of diagnosis. For the forward models to producereasonable estimates on the tumor invasion profile we would firstly need informa-tion on the location and time of the carciogenesis, secondly the information on thenon-degraded diffusive properties of the tissue surrounding it and lastly the correct parametrisation on a per-patient basis. We want to emphasise that the temporal dy-namics of the tumor growth, although scientifically interesting, are not necessarilyrelevant in the medical treatment planning. The information needed for treatmentwith the current techniques is an accurate description of the tumor density field attime of diagnosis .The derivation of the stationalisation term for the one dimensional Fisher equationis not strictly transferable to tumor models which incorporate medical data or haveadditional advection terms, but the results from section 4.3 indicate that they maystill produce reasonable predictions. If the underlying diffusion matrix field includedstrong inhomogeneities inducing strong advective terms, the procedure might lose itsvalidity, however an investigation of the given 3D DTI dataset shows that the peclet-number relating the drift and diffusion strength as given in (8a), τ = | v L || D | , staysmostly below 0.3 when D t is derived from the DTI data by simple scaling (14). Mosttumor models show travelling wave characteristics, where the main physical effectsinclude diffusion and nonlinear growth. We recommend close inspection of both thegradient distributions and the peclet numbers if the stationalised model should beextended. The stationalisation procedure should retain its applicability as long as thegradient distribution retains its underlying characteristic, and the distance betweenmedical segmentation and the border induced by the level-sets is not chosen too large.In the case that the presented approximation fails, it may also be possible to introducemore elaborate numerical ways to fit the necessary penalty term locally.We compared the level-sets of the forward simulation with those of the station-ary simulation and found characteristic distances between them of about L B ≈ v for a given set of forward model parameters. A medicalpractitioner might choose the actual level-set value to replicate the current practiceof treating a 2 cm radius around the bulk tumor and then use qualitative informa-tion in the form of locally recommending an extension or retraction of the treatmentradius. The stationary model will correctly capture the effect of the material proper-ties, as presented in Figure 8. Where the tissue is more diffusive, the level-set on u s will overextend the 2 cm radius and where the diffusivity is small, more brain tissuemight be left untreated. If the underlying tumor model should be extended, all effectsincreasing- or diminishing the spread of glia cells will be reflected accordingly in thestationalised solution.Should time-series datasets become available, then the stationalisation may beused to estimate the initial condition g ( x ) for a forward simulation from the firstdataset. An initial condition calculated in this way should be closer to a presumablysmooth real density profile than using a stepped profile with steep gradients.Naturally the computation of solutions to the stationary problem take less timethan a full forward simulation. While compute servers are usually available in aca-demic institutions, the possibility to calculate the results on a regular consumer pcwith only short computation times is important in order to transfer such methodsinto clinical application. The stationary simulation allows for the computation of sets stimating the extent of glioblastoma invasion 19 of solutions for varying parameters within a short timeframe. We present exemplaryruntimes for the camino dataset on recent hardware for the two cases in Table 2. simulation type hardware runtime Table 2
Runtimes for the testcases in two- and three dimensions on different hardware. v ≥ √ Dρ ) in the one dimensionalFisher-KPP equation may be an indicator for how a localized penalty parameter couldbe improved. Instead of choosing a constant v globally, it might be possible to set apenalty factor linearly combined with local information. Thereby incorporating thelocal increases and decreases in diffusivity into the stationalisation.In section 4.3 we used a real dataset, but a comparatively primitive tumor model.It should be possible to extend the stationalisation procedure to incorporate additionaleffects like chemo- or haptotaxis as long as the dynamic of producing travelling wavesolutions of sigmoidal shape is not altered by the additions. In the example in section4.3, we used a level-set on the forward simulation as the internal Dirichlet constraintsfor the stationalisation. In a medical setting one would directly use the medical seg-mentation from the DTI/MRI modalities. There are promising advances in generatingtumor segmentations in an automated fashion, e.g. BraTumIA (Porz et al., 2014).Automating the process of the segmentation opens up the possibility to use a fullyautomated process to advise the treatment planning in real patients.We showed how well a stationalised formulation would perform compared to aforward simulation if all the necessary information were present. The fact that time-series datasets are largely unavailable and therefore no parametrizations can be de-rived from them, makes direct comparisons between existing tumor models difficult.It is, however possible to compare the stationalized versions of existing models withonly the datasets from the time of diagnosis. One could then compare these levelsetsto the clinical target volume (CTV) regularly produced in medical practice. This of-fers an attractive approach to perform a model comparison for a wide range of tumormodels proposed in the literature. Acknowledgements
The authors would like to express gratitude towards the BMBF (Bundesministerium f¨ur Bildung undForschung) for funding the GlioMath project (funding code: 05M16PMA), and all collaborating partnerswithin.
References
Ablowitz MJ, Zeppetella A (1979) Explicit solutions of fisher’s equation for a spe-cial wave speed. Bulletin of Mathematical Biology 41(6):835–840, DOI 10.1007/BF02462380Bastian P, Blatt M, Dedner A, Engwer C, Kl¨ofkorn R, Kornhuber R, Ohlberger M,Sander O (2008a) A Generic Grid Interface for Parallel and Adaptive ScientificComputing. Part II: Implementation and Tests in DUNE. Computing 82(2–3):121–138, DOI 10.1007/s00607-008-0004-9Bastian P, Blatt M, Dedner A, Engwer C, Kl¨ofkorn R, Ohlberger M, Sander O(2008b) A Generic Grid Interface for Parallel and Adaptive Scientific Comput-ing. Part I: Abstract Framework. Computing 82(2–3):103–119, DOI 10.1007/s00607-008-0003-xBastian P, Heimann F, Marnach S (2010) Generic implementation of finite elementmethods in the Distributed and Unified Numerics Environment (DUNE). Kyber-netika 46:294–315Blatt M, Bastian P (2007) The iterative solver template library. In: Kagstr¨om B, Elm-roth E, Dongarra J, Wasniewski J (eds) Applied Parallel Computing – State of theArt in Scientific Computing, Springer, Berlin/Heidelberg, pp 666–675Blatt M, Burchardt A, Dedner A, Engwer C, Fahlke J, Flemisch B, Gersbacher C,Gr¨aser C, Gruber F, Gr¨uninger C, Kempf D, Kl¨ofkorn R, Malkmus T, Mthing S,Nolte M, Piatkowski M, Sander O (2016) The Distributed and Unified NumericsEnvironment, Version 2.4. Archive of Numerical Software 4(100):13–29, DOI10.11588/ans.2016.100.26526Brazhnik PK, Tyson JJ (1999) On traveling wave solutions of fisher’s equation in twospatial dimensions. SIAM Journal on Applied Mathematics 60(2):371–391Caragher S, Chalmers AJ, Gomez-Roman N (2019) Glioblastomas next top model:Novel culture systems for brain cancer radiotherapy research. Cancers 11(1), DOI10.3390/cancers11010044Chang EL, Akyurek S, Avalos T, Rebueno N, Spicer C, Garcia J, Famiglietti R,Allen PK, Chao KC, Mahajan A, Woo SY, Maor MH (2007) Evaluation of peritu-moral edema in the delineation of radiotherapy clinical target volumes for glioblas-toma. International Journal of Radiation Oncology,Biology,Physics 68(1):144 –150, DOI 10.1016/j.ijrobp.2006.12.009Conte M, Gerardo-Giorda L, Groppi M (2020) Glioma invasion and itsinterplay with nervous tissue and therapy: A multiscale model. Journalof Theoretical Biology 486:110088, DOI https://doi.org/10.1016/j.jtbi.2019.110088, URL
Cook PA, Bai Y, Nedjati-Gilani S, Seunarine KK, Hall MG, Parker GJ, Alexan-der DC (2006) Camino: Open-source diffusion-mri reconstruction and process- stimating the extent of glioblastoma invasion 21 ing. URL
Corbin G, Hunt A, Klar A, Schneider F, Surulescu C (2018) Higher-order models forglioa invasion: from a two-scale description to effective equations for mass densityand momentum. Math Models Meth Appl SciEngwer C, Hillen T, Knappitsch M, Surulescu C (2015) Glioma follow white mattertracts: a multiscale dti-based model. J Math Biol 71:551–582Engwer C, Hunt A, Surulescu C (2016) Effective equations for anisotropic gliomaspread with proliferation: a multiscale approach. J Math Med Biol 33:435–459Fischer RA (1937) The wave of advance of advantageous genes. Annals of Eugenics7(4):355–369, DOI 10.1111/j.1469-1809.1937.tb02153.xHarpold HL, Alvord EC Jr, Swanson KR (2007) The evolution of mathematicalmodeling of glioma proliferation and invasion. Journal of Neuropathology Experi-mental Neurology 66(1):1–9, DOI 10.1097/nen.0b013e31802d9000, URL http://dx.doi.org/10.1097/nen.0b013e31802d9000
Hillen T, Painter KJ (2013) Transport and Anisotropic Diffusion Models for Move-ment in Oriented Habitats, Springer Berlin Heidelberg, Berlin, Heidelberg, pp177–222. DOI 10.1007/978-3-642-35497-7 7, URL https://doi.org/10.1007/978-3-642-35497-7_7
Hunt A (2018) Dti-based multiscale models for glioma invasion. doctoralthesis, Tech-nische Universit¨at Kaiserslautern, URL http://nbn-resolving.de/urn:nbn:de:hbz:386-kluedo-53575
Hunt A, Surulescu C (2017) A multiscale modeling approach to glioma in-vasion with therapy. Vietnam Journal of Mathematics 45(1):221–240,DOI 10.1007/s10013-016-0223-x, URL https://doi.org/10.1007/s10013-016-0223-x
Jan Kelkel CS (2011) On some models for cancer cell migration throughtis-sue networks. Mathematical Biosciences & Engineering 8:575, DOI 10.3934/mbe.2011.8.575, URL http://aimsciences.org//article/id/0fac622a-15b3-4a5d-9a14-4414a389a19f
Jbabdi S, Mandonnet E, Duffau H, Capelle L, Swanson KR, PlgriniIssac M, GuillevinR, Benali H (2005) Simulation of anisotropic growth of lowgrade gliomas usingdiffusion tensor imaging. Magnetic Resonance in Medicine 54(3):616–624, DOI10.1002/mrm.20625, URL https://doi.org/10.1002/mrm.20625
JDMurray (2007) Mathematical Biology 1: An Introduction. SpringerKolmogoroff A, Petrovsky I, Piscounoff N (1988) Study of the diffu-sion equation with growth of the quantity of matter and its applica-tion to a biology problem. In: Pelc P (ed) Dynamics of Curved Fronts,Academic Press, San Diego, pp 105 – 130, DOI https://doi.org/10.1016/B978-0-08-092523-3.50014-9, URL
Konukoglu E, Sermesant M, Clatz O, Peyrat JM, Delingette H, Ayache N (2007) Arecursive anisotropic fast marching approach to reaction diffusion equation: Appli-cation to tumor growth modeling. Information Processing in Medical Imaging pp687–699
Konukolu E, Clatz O, Bondiau PY, Delingette H, Ayache N (2006) Extrapolatingtumor invasion margins for physiologically determined radiotherapy regions. In:Larsen R, Nielsen M, Sporring J (eds) Medical Image Computing and Computer-Assisted Intervention – MICCAI 2006, Springer Berlin Heidelberg, Berlin, Hei-delberg, pp 338–346Mandonnet E, Delattre JY, Tanguy ML, Swanson KR, Carpentier AF, Duffau H,Cornu P, Van Effenterre R, Alvord Jr EC, Capelle L (2003) Continuous growthof mean tumor diameter in a subset of grade ii gliomas. Annals of Neurology53(4):524–528, DOI 10.1002/ana.10528Oraiopoulou ME, Tzamali E, Tzedakis G, Liapis E, Zacharakis G, Vakis A, Papa-matheakis J, Sakkalis V (2018) Integrating in vitro experiments with in silico ap-proaches for glioblastoma invasion: the role of cell-to-cell adhesion heterogeneity.Scientific reports 8(1):16200–16200, URL
Othmer HG, Hillen T (2000) The diffusion limit of transportequations derived from velocity-jump processes read more:https://epubs.siam.org/doi/10.1137/s0036139999358167. SIAM Journal onApplied Mathematics SIAM J. Appl. Math., 61:751775. (25 pages)Othmer HG, Hillen T (2002) The diffusion limit of transport equations ii: Chemo-taxis equations. SIAM Journal on Applied Mathematics 62:12221250, DOI https://doi.org/10.1137/S0036139900382772, URL https://doi.org/10.1137/S0036139900382772
Painter K, Hillen T (2013) Mathematical modelling of glioma growth: The useof diffusion tensor imaging (dti) data to predict the anisotropic pathways ofcancer invasion. Journal of Theoretical Biology 323:25 – 39, DOI https://doi.org/10.1016/j.jtbi.2013.01.014, URL
Patel V, Hathout L (2017) Image-driven modeling of the proliferation and necro-sis of glioblastoma multiforme. Theoretical Biology and Medical Modelling14(1):10, DOI 10.1186/s12976-017-0056-7, URL https://doi.org/10.1186/s12976-017-0056-7
Porz N, Bauer S, Pica A, Schucht P, Beck J, Verma RK, Slotboom J, ReyesM, Wiest R (2014) Multi-modal glioblastoma segmentation: Man versus ma-chine. DOI 10.1371/journal.pone.0096873, URL https://doi.org/10.1371/journal.pone.0096873
Rutter EM, Stepien TL, Anderies BJ, Plasencia JD, Woolf EC, Scheck AC,Turner GH, Liu Q, Frakes D, Kodibagkar V, Kuang Y, Preul MC, KostelichEJ (2017) Mathematical analysis of glioma growth in a murine model. Sci-entific reports 7(1):2508–2508, URL
Sathornsumetee S, Reardon DA, Desjardins A, Quinn JA, Vredenburgh JJ, Rich JN(2007) Molecularly targeted therapy for malignant glioma. Cancer 110(1):13–24,DOI 10.1002/cncr.22741Silbergeld DL, Chicoine MR (1997) Isolation and characterization of humanmalignant glioma cells from histologically normal brain. Journal of Neuro-surgery 86(3):525 – 531, URL https://thejns.org/view/journals/ stimating the extent of glioblastoma invasion 23 j-neurosurg/86/3/article-p525.xml
Swanson KR, Alvord EC, Murray JD (2000) A quantitative model for differentialmotility of gliomas in grey and white matter. Cell Proliferation 33(5):317–329,DOI 10.1046/j.1365-2184.2000.00177.xSwanson KR, Rostomily RC, Alvord EC (2008) A mathematical modelling toolfor predicting survival of individual patients following resection of glioblastoma:a proof of principle. British Journal of Cancer 98(1):113–119, URL https://doi.org/10.1038/sj.bjc.6604125
Tracqui P, Cruywagen GC, Woodward DE, Bartoo GT, Murray JD, Alvord JrEC (1995) A mathematical model of glioma growth: the effect of chemother-apy on spatio-temporal growth. Cell Proliferation 28(1):17–31, DOI 10.1111/j.1365-2184.1995.tb00036.xvan der Vorst HA (1992) Bi-cgstab: A fast and smoothly converging vari-ant of bi-cg for the solution of nonsymmetric linear systems. SIAM Journalon Scientific and Statistical Computing 13(2):631–644, DOI 10.1137/0913035,URL https://doi.org/10.1137/0913035 , https://doi.org/10.1137/0913035https://doi.org/10.1137/0913035