A GPU implementation of the Discontinuous Galerkin method for simulation of diffusion in brain tissue

by   Daniel Cervantes, et al.

In this work we develop a methodology to approximate the covariance matrix associated to the simulation of water diffusion inside the brain tissue. The computation is based on an implementation of the Discontinuous Galerkin method of the diffusion equation, in accord with the physical phenomenon. The implementation in in parallel using GPUs in the CUDA language. Numerical results are presented in 2D problems.



There are no comments yet.


page 5

page 9

page 10


Implementation of high-order, discontinuous Galerkin time stepping for fractional diffusion problems

The discontinuous Galerkin dG method provides a robust and flexible tech...

A hybridizable discontinuous Galerkin method for electromagnetics with a view on subsurface applications

Two Hybridizable Discontinuous Galerkin (HDG) schemes for the solution o...

Practical computation of the diffusion MRI signal of realistic neurons based on Laplace eigenfunctions

The complex transverse water proton magnetization subject to diffusion-e...

Molecular dynamic simulation of water vapor interaction with blind pore of dead-end and saccate type

One of the varieties of pores, often found in natural or artificial buil...

A comparison of matrix-free isogeometric Galerkin and collocation methods for Karhunen–Loève expansion

Numerical computation of the Karhunen–Loève expansion is computationally...

Moving Mesh with Streamline Upwind Petrov-Galerkin (MM-SUPG) Method for Convection-Diffusion Problems

We investigate the effect of the streamline upwind Petrov-Galerkin metho...

Parameter-robust Multiphysics Algorithms for Biot Model with Application in Brain Edema Simulation

In this paper, we develop two parameter-robust numerical algorithms for ...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1. Introduction

The study of the diffusion phenomenon in porous media is an active line of research with many applications in the applied sciences. The lack of analytical solutions for arbitrary shaped domains requires the use of numerical solvers to compute reliable solutions. The analytical solutions are restricted to simple–shaped geometries as spheres, parallel walls and cylinders [14]

. However, real–life applications require to compute estimations of the diffusion on complex domains.

An important research area in the medical sciences, and the focus of this work, is the simulation of water diffusion inside the brain tissue (white and grey matter). The problem is to characterize the hydrogen molecular displacement due to brownian motion. Literature on the subject is vast, let us present a short review.

The probability of the displacement

(in micrometers) of an ensemble of water molecules in a given time (in miliseconds) is summarized by the Ensemble Average Propagator (EAP) [18]. The Diffusion Weighted modality of the Magnetic Resonance allows to capture data to infer the molecular displacement on in vivo and ex vivo experiments [2, 3]. The quantification of molecular displacement of water molecules on in vivo patientes allows to infer properties of the microstructure of brain tissue [20, 16]

. The information above is used to detect tissue disruption associated with deseases: for instance, a premature reduction of the cellular volume is an indicator of neuron death. The numerical simulation of the water diffusion allows to generate useful data for the validation of model fitting

[10], generation of novel theories about the diffusion properties with analytical representations [4], model improvement [11], etc. In this problem the domain is composed of cell bodies (soma, axons, neurites, glial cells, vessels, etc). Despite the fact that some of the cell structures can be approximated with simple geometries, in general the domains are much more geometrically complex, in particular, the extracellular spaces presents arbitrary shapes (similar to a ”gruyere–cheese shape”) hence numerical approximations are required to simulate the molecular diffusivity [12].

Numerical solvers to compute the EAP of the hydrogen displacement are varied. For methods based on Monte Carlo diffusion simulators see [12, 10, 11, 12]

and for partial differential–equation based solvers

[15, 8]. However, the need to produce massive simulation data for: a) simulating experiments for different machine parameters (magnitude of magnetic gradientes, experimental times, etc) [9], b) validating complex models [16], c) training automatic learning algorithms [7], etc., requires to produce accurate simulations in optimal computational times. Nowadays, the diffusion simulators in the state–of–the–art [12] requiere from minutes to days to estimate the diffusion phenomena on complex domains [19].

This leads to the aim of this work, the approximation of the extracellular diffusivity profile on a disordered medium. On a disordered model of cylindrical brain axons, with the diameters computed from a Gamma distribution

[1, 12] the whole process can be characterized by the corresponding 3D covariance matrix . The eigenvalues (matched with the corresponding eigenvectors) indicate the magnitude and orientational dependency, such that, it is possible to infer extracellular microstructure features as: the main orientation of the axon bundles, the percentage of the volume occupied by neurons (intra cellular signal fraction), the amount of diffusion anisotropy of the tissue (fractional anisotropy), among other descriptors. In the DWMR medical literature, this covariance matrix is computed by the diffusion tensor (DT) from the MR signal [17, 16]. Those descriptors computed from have been correlated with several brain damages and diseases [13].

Consequently, the approximation of the referred covariance matrix is of great interest. We shall introduce a numerical methodology for approximation. As a first step we develop a 2D version of the problem.

The covariance matrix is obtained from a gaussian density, formed by averaging densities which result from the solution of diffusion equations. Our main contribution is to solve these diffusion equations using and ah hoc implementation of the Discontinuous Galerkin (DG) method, see [5] for a thorough discussion.

Our implementation takes advantage of the underlying physical and geometric properties of the problem as posed by the clinicians. In practice, the so called substrates are squared pixel domains allowing for a uniform mesh. Also there is an assumption of no diffusivity between the axons and the extracellular region. The common approximation solves the diffusion equation only in the extracellular region imposing a zero Neumann condition. An alternative is to use the numerical fluxes, a main feature of the DG method, to propose an interaction between the axons and the intercellular region. With this interaction, we are solving the diffusion equation in the whole domain. Thus null computations are carried out in the axons. This, apparently redundant strategy, free us of boundary handling, a computationally expensive task in Galerkin type methods. But more importantly, it allows for the full strength of the DG method, to carry out the time update of the solution in parallel for all elements. A small ODE problem is solved in each element with no communication. Consequently, a GPU implementation is appropriate.

The outline is as follows.

In Section 2, the Initial Boundary Value Problem (IVBP) of the diffusion equation associated to the phenomenon is introduced. The a substrate for study is described for a simulated ex vivo experiment.

The basics of the DG method are presented in Section 3. Therein, the physical and geometric properties of the case study are used to tune our DG-GPU implementation.

In Section 4 a scheme to approximate the diffusion encoding covariance matrix is introduced. Performance is illustrated with free diffusion and the case study associated to an ex vivo experiment. Conclusions and a brief discussion on future work close our exposition.

2. Water diffusion inside the brain tissue

In this section we introduce a substrate with such a realistic shape and properties. First, we discuss the Initial Boundary Value Problem (IBVP) for the diffusion equation to be used throughout.

2.1. The diffusion problem

In practice, one models a substrate occupying a domain, which is a medium comprised of two regions, the axons and the extracellular complement. The former is regarded as non diffusive and the latter a region with constant diffusion. It is assumed that the boundary between both regions is reflecting. Initial pulses are prescribed in the extracellular region, far away from the outer boundary which is modelled as a perfect absorber.

Let be the domain occupied by the substrate. This domain is the union of two intertwined adjacent regions, and . These are respectively, the the axon and extracellular regions.

The IVBP consists on finding that solves the diffusion equation


given Cauchy data


and boundary values


Here is the outer unit normal to . Notice that .

The Neumann boundary condition (3) corresponds to a reflecting boundary, whereas the Dirichlet boundary condition (4) to that of a perfect barrier.

The discontinuous diffusion is given by


for a positive constant .

2.2. A case study

The substrate under consideration consists of 1901 non-overlapping circles which represent the axons and and we only take into account the regions that are within a square, the domain , that measures on the side (Figure 1). The radius of the circles is within the range of up to . The extracellular region is the exterior of the circles and the diffusion coefficient is set to .

The ex vivo coefficient oscillates between 450 and 600, depending on the temperature and the substances of the medium, [6]. We use the smaller one to be able to use small substrates.

Figure 1. The substrate consists of 1901 non-overlapping circles.

It is apparent that numerical approximations are required to simulate the molecular diffusivity is this rather complex porous medium.

3. DG-CUDA solution of the heat equation

3.1. The Discontinuous Galerkin Method

Let be partitioned into non overlapping polygonal elements, e.g. a triangulation. Let us denote by one of such elements, see Figure 2.

Figure 2. Integration element and its normal.

A defining feature of the DG method is to reduce the PDE to a first order system. Consequently, let us consider


Multiplying (6) by the test functions , and , we obtain after integrating by parts,

A second feature of the DG method is the element-wise approximation of the unknown and . Continuity is not enforced at the boundary of adjacent elements. Consequently, the boundary terms , are replaced by boundary fluxes , . As customary, the superscript denotes limits from the interior of , and the superscript limits from the exterior.

This yields in element


Let be an approximation of any of the scalar funcions , ,. Within the approximation is given by


For a triangulation, the functions

are chosen as in the Finite Element Method with lagrangian interpolation.

It is assumed that and are known in (7). Hence, we are led to solve an differential-algebraic system for and . The solution in time is advanced by a Runge-Kutta method.

Remark. We stress that the solution of the differential-algebraic system (7), is solved independently for each element. In practice suffices. Thus, we have small systems to solve that do not exchange information in each time step.

3.2. Numerical flux

There is no preferred direction of propagation in the heat equation, thus for a central flux is considered, namely

A physical assumption is that there is no flow between axons and the extracellular region. Consequently, for , we propose the numerical flow

This is coined for the problem under consideration. For instance, if and

, the harmonic mean forces a zero Neumann condition, hence there is no flow trough the boundary of the element

as expected.

3.3. Cuda implementation

Meshing is a time consuming task in Galerkin type methods, as it is boundary conditions handling when assembling the local systems. In our case, the domain is divided in square pixels which we use to our advantage. More precisely, these squares are separated by the diagonal in two triangles. A triangle constructed in this fashion, is the basic element for discretization. See Figure 3

Figure 3. Meshing.

Also in this MRI application, the heat equation is solved in the extracellular region where Cauchy data is prescribed. Apparently, there is no need to consider the axon region . Nevertheless, we solve the heat equation in elements contained in where the contribution to the solution is null.

We are led to balanced computations on every element. As pointed out in the remark above, the calculations in (7) are element independent when updating time. Consequently the main ingredients for parallel processing using GPUs are met. Namely, small balanced computations with no exchange of information between processors.

Parallel processing using GPUs is implemented in the CUDA language in a DELL laptop with hardware:

CPU: Intel(R) Core(TM) i7-6820HQ CPU @ 2.70 GHz

GPU: NVIDIA Corporation GM107GLM [Quadro M1000M]

4. The Gaussian profile of the extra-cellular diffusivity

In this section we introduce a scheme to approximate the covariance matrix of the EAP. Hence describing the Gaussian profile of the extra-cellular diffusivity. Numerical results are also shown for the case of free diffusion, and a comparison with MCMC.

4.1. The scheme

Let us choose , randomly and uniformly in . Set a final time .

  • For , let be the DG-Cuda solution of the IVBP (1)(4), where the Cauchy data is the Dirac’s delta function supported in .

  • For , let be obtained from by centering and normalization to yield a density function

  • Construct the mixture model

  • Fit a Gaussian density to , that is, determine a covariance matrix , such that

The covariance matrix is given by,

It is computed by quadrature rules using , the values of the numerical solution at the nodes in the mesh.

4.2. Numerical results

The substrate under study occupies a square domain of side . Cauchy data is given within a centered box of side . The observation time is seconds. At this given time the outer boundary is not reached by diffusion. In all cases a square grid of mesh is used.

Free diffusion. In this simple example, We solve the Heat Equation with Cauchy data a Dirac’s delta supported at the origin. The uniform diffusion coefficient is . Graphical results in Figure 4.

Figure 4. Left: 2D free diffusion. Right: 1D view.

Let us compare the free diffusion (without axons) approximation, versus the analytical solution. The latter is a bivariate Gaussian function defined in (4.1), where the covariance matrix is

Taking and , non zero coefficients in the covariance matrix are equal to .

The DG-CUDA Gaussian matrix fit for is

An alternative construction of the covariance matrix is by means of Monte Carlo Diffusion Simulation. See [12] for details. We just list the corresponding data in their notation.

Diffusion constant , duration of the diffusion simulation. is the number of time steps in the simulation. The step length is obtained from the relation

The obtained MC Gaussian matrix is

To gauge the approximations we compute the least squares residual


In both cases the approximation of the Gaussian function is highly accurate. The least squares residual (9) is of the order . For practical purposes, the Gaussian density functions coincide. But the DG-CUDA approximation is structurally more consistent. The matrix is symmetric, the values in the diagonal coincide and the other terms are near zero.

Case study. The axon region is defined by 1901 axons. The diffusion coefficient as in (5). The scheme above is applied to a mixture of densities. An instance of one PDE solution is shown in Figure 5. The solution of the full scheme is in Figure 6.

Figure 5. Left: One PDE solution. Right: Normalized and centered solution for one PDE.
Figure 6. Left: 2D free diffusion. Right: 1D view.

The covarince matrix by the DG-CUDA algorithm is

Finally, let us summarize execution time in the following table

Time stps No. Deltas CPU GPU
400400 5184 1 6585.47 sec. 5.13649 sec.
400 400 5184 37 21877.252 sec. 188.597 sec.

We remark that regardless of the diffusion coefficient, the GPU process is more than 1000 times faster for the solution of one instance of the heat equation with a single delta as Cauchy data.

5. Conclusions and future work

In this work we have provided a methodology for the efficient computation of the diffusivity properties in porous media. These results can be used for the simulation of the DWMR signal given the numerical estimation of the EAP. Such a computation can be analytically performed by the application of the Fourier transform on the EAP, however, simplistic assumptions about the machine model have to be made

[18]. In order to provide a useful tool for the medical researchers, the numerical simulator of the signal should take into account the signal changes associated to a realistic MR machine parameters, as the non squared (but trapezoidal) and the finite duration of the magnetic pulses. It is our contention that our methodology is versatile to include this complexities in a 3D extension. The latter is part of our current and future work.


D. Cervantes and M. A. Moreles thank the support of ECOS-NORD through the project: 000000000263116/M15M01. M. A. Moreles also acknowledges partial support from the project CIB CONACYT 180723.


  • [1] F. Aboitiz, A.B. Scheibel, R.S.Fisher, and E.Zaidel. Fiber composition of the human corpus callosum. Brain Res., 143(53), 1992.
  • [2] D.C. Alexander. An introduction to computational diffusion mri: the diffusion tensor and beyond. In Hagen H. Weickert J., editor, Visualization and Processing of Tensor Fields, chapter 10, pages 83–106. Springer Verlag, 2006.
  • [3] D.C. Alexander. Modelling, fitting and sampling in diffusion mri. Magnetic Resonance, 103:255–260, 2009.
  • [4] L. M. Burcaw, E. Fieremans, and D. S. Novikov. Mesoscopic structure of neuronal tracts from time-dependent diffusion. NeuroImage, 114:18–37, 2015.
  • [5] Bernardo Cockburn, George E Karniadakis, and Chi-Wang Shu. Discontinuous Galerkin methods: theory, computation and applications, volume 11. Springer Science & Business Media, 2012.
  • [6] Tim B Dyrby, Matt G Hall, Maurice Ptito, Daniel Alexander, et al. Contrast and stability of the axon diameter index from microstructure imaging with diffusion MRI. Magnetic Resonance in Medicine, 70(3):711–721, 2013.
  • [7] Nedjati-Gilani Gemma L. et. al. Machine learning based compartment models with permeability for white matter microstructure imaging. Neuroimage, 150:119?135, 2017.
  • [8] Moroney B. F., Stait-Gardner T., Ghadirian B., Yadav N. N., and Price WS. Numerical analysis of nmr diffusion measurements in the short gradient pulse limit. Journal of Magnetic Resonance, 234:165–175, 2013.
  • [9] Uran Ferizi, Torben Schneider, Thomas Witzel, Lawrence L. Wald, Hui Zhang, Claudia A.M. Wheeler-Kingshott, and Daniel C. Alexander. White matter compartment models for in vivo diffusion MRI at 300 mt/m. NeuroImage, 118:468–483, 2015.
  • [10] Els Fieremans, Dmitry S. Novikov, Jens H. Jensen, and Joseph A. Helpern. Monte carlo study of a two-compartment exchange model of diffusion. NMR in biomedicine, 23 7:711–24, 2010.
  • [11] Nima Gilani, Paul N. Malcolm, and Glyn Johnson. An improved model for prostate diffusion incorporating the results of monte carlo simulations of diffusion in the cellular compartment. NMR in biomedicine, 30 12, 2017.
  • [12] M.G. Hall and D.C. Alexander. Convergence and parameter choice for monte-carlo simulations of diffusion mri. IEEE TRANSACTIONS ON MEDICAL IMAGING, 28(9), 2009.
  • [13] Derek Jones. Diffusion MRI theory methods and applications. Oxford University Press, 2011.
  • [14] C.H. Neuman. Spin echo of spins diffusing in a bounded medium. Chemical Physics, 60(11):4508?4511, 1974.
  • [15] Dang Van Nguyen, Jing-Rebecca Li, Denis Grebenkov, and Denis Le Bihan. A finite elements method to solve the bloch-torrey equation applied to diffusion magnetic resonance imaging. J. Comput. Physics, 263:283–302, 2014.
  • [16] E. Panagiotaki, T. Schneider, B. Siow, M.G. Hall, M.F. Lythgoe, and D.C. Alexander. Compartment models of the diffusion mr signal in brain white matter: A taxonomy and comparison. NeuroImage, 2012.
  • [17] C Pierpaoli, P Jezzard, P J Basser, A Barnett, and G Di Chiro. Diffusion tensor MR imaging of the human brain. Radiology, 201(3):637–648, 1996.
  • [18] W.S. Price. Pulsed-field gradient nuclear magnetic resonance as a tool for studying translational diffusion: Part 1. basic theory. G-Animal’s Journal, 1997.
  • [19] Alonso Ramirez-Manzanares, Mario Ocampo-Pineda, Jonathan Rafael-Patino, Giorgio Innocenti, Jean-Philippe Thiran, and Alessandro Daducci. Quantifying diameter overestimation of undulating axons from synthetic dw-mri. In Procc. In Annual Meeting of the International Society of Magnetic Resonance in Medicine, page 748. Annual Meeting of the SMRM, París, France., 2018.
  • [20] P. van Gelderen, D. DesPers, C.M. van Zijl, and C.T.W. Moonen. Evaluation of restricted diffusion in cylinders. phosphocreatine in rabbit leg muscle. Magnetic Resonance, 103:255–260, 1994.