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

TL;DR

This paper presents a GPU-accelerated Discontinuous Galerkin method for simulating water diffusion in brain tissue, enabling efficient approximation of the covariance matrix in diffusion studies.

Contribution

It introduces a parallel GPU implementation of the Discontinuous Galerkin method tailored for brain tissue diffusion simulation, improving computational efficiency.

Findings

GPU implementation achieves significant speedup

Numerical results validate the method in 2D cases

Method accurately approximates covariance matrices

Abstract

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.

Equations

$$\frac{\partial u}{\partial t}=\nabla \cdot \left(k(x)\nabla u\right),(x,t)\in {\Omega }_e\times (0,T),$$

$$u(x,0)=u_0(x),x\in {\Omega }_e,$$

$$\frac{\partial u}{\partial \mathbf{n}}=0,(x,t)\in \partial {\Omega }_e\times (0,T),$$

$$u(x,t)=0,(x,t)\in \partial \Omega \times (0,T).$$

$$k(x)=\left\{\begin{array}[]{l}0,\quad x\in\Omega_{a}\\ \\ k_{0},\quad x\in\Omega_{e}\end{array}\right.$$

$$\{\begin{array}{ll}\mathbf{q}& =\nabla u\\ \frac{\partial u}{\partial t}& =\nabla \cdot \left(k\mathbf{q}\right).\end{array}$$

$${\int }_{K^{-}}\mathbf{q}\cdot \mathbf{v}d\mathbf{x}$$

$$={\int }_{\partial K^{-}}(u\mathbf{v}\cdot \mathbf{n})d\mathbf{s}-{\int }_{K^{-}}u\cdot \nabla \mathbf{v}d\mathbf{x}$$

$${\int }_{K^{-}}{u}_{t}vd\mathbf{x}$$

$$={\int }_{\partial K^{-}}v(k\mathbf{q})\cdot \mathbf{n}d\mathbf{s}-{\int }_{K^{-}}k\mathbf{q}\nabla \cdot vd\mathbf{x}$$

$$\begin{array}[]{rcl}\int_{K^{-}}\mathbf{q}\cdot\mathbf{v}d\mathbf{x}&=&\int_{\partial K^{-}}\mathbf{v}\cdot h_{u,K^{-}}(u^{-},u^{+},\mathbf{n}^{-})d\mathbf{s}-\int_{K^{-}}u\cdot\nabla\mathbf{v}d\mathbf{x}\\ &&\\ \int_{K^{-}}u_{t}vd\mathbf{x}&=&\int_{\partial K^{-}}vh_{\mathbf{q},K^{-}}(\mathbf{q}^{-},\mathbf{q}^{+},\mathbf{n}^{-})d\mathbf{s}-\int_{K^{-}}k\mathbf{q}\nabla\cdot vd\mathbf{x}\end{array}$$

$$w(x,t)=\sum _{j=0}^p{w}_j(t)N_j(x).$$

$$h_{u,K^{-}}(u^{-},u^{+},{\mathbf{n}}^{-})=\frac{u^{-}+u^{+}}{2}{\mathbf{n}}^{-}.$$

$$h_{\mathbf{q},K^{-}}({\mathbf{q}}^{-},{\mathbf{q}}^{+},{\mathbf{n}}^{-})=\frac{2k^{-}{k}^{+}}{k^{-}+k^{+}}\frac{1}{2}({\mathbf{q}}^{-}+{\mathbf{q}}^{+})\cdot {\mathbf{n}}^{-}.$$

$$u(x,y)=\frac{1}{m}\sum _{i=1}^m{u}_i(x,y,T).$$

$$u(x,y;\Sigma )\approx \frac{1}{2\pi \sqrt{|\Sigma }|}\exp \left(-\frac{1}{2}\left[((x,y)-({\mu }_x,{\mu }_y){)}^T{\Sigma }^{-1}((x,y)-({\mu }_x,{\mu }_y))\right]\right),$$

$$\Sigma =\left(\begin{array}{cc}E[(EX-{\mu }_x)(EX-{\mu }_x)]& E[(EX-{\mu }_x)(EY-{\mu }_y)]\\ E[(EY-{\mu }_y)(EX-{\mu }_x)]& E[(EY-{\mu }_y)(EY-{\mu }_y)]\end{array}\right).$$

$$\Sigma =\left(\begin{array}{cc}2Tk& 0\\ 0& 2Tk\end{array}\right).$$

$$\left(\begin{array}{cc}32.53& 0.00038\\ 0.00038& 32.53\end{array}\right)$$

$$l=\sqrt{4\frac{D{t}_s}{T}}\approx 0.11\mu m.$$

$$\left(\begin{array}{cc}32.487358329& -0.075889940\\ -0.075889940& 32.378282498\end{array}\right)$$

$$\sum _i\sum _j[u(x_i,y_i;\Sigma )-u_{ij}{]}^2.$$

$$\left(\begin{array}{cc}19.50& -0.0088\\ -0.0088& 19.50\end{array}\right)$$

Full text

1. 1 Introduction

2. 2 Water diffusion inside the brain tissue

3. 2.1 The diffusion problem

4. 2.2 A case study

5. 3 DG-CUDA solution of the heat equation

6. 3.1 The Discontinuous Galerkin Method

7. 3.2 Numerical flux

8. 3.3 Cuda implementation

9. 4 The Gaussian profile of the extra-cellular diffusivity

10. 4.1 The scheme

11. 4.2 Numerical results

12. 5 Conclusions and future work

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 $x$ (in micrometers) of an ensemble of water molecules in a given time $\Delta$ (in miliseconds) is summarized by the Ensemble Average Propagator (EAP) $P(x,\Delta)$ [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 $\Sigma$. The $\Sigma$ 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 $\Sigma$ 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 $\Omega$ be the domain occupied by the substrate. This domain is the union of two intertwined adjacent regions, $\Omega_{a}$ and $\Omega_{e}$. These are respectively, the the axon and extracellular regions.

The IVBP consists on finding $u$ that solves the diffusion equation

[TABLE]

given Cauchy data

[TABLE]

and boundary values

[TABLE]

[TABLE]

Here $\mathbf{n}$ is the outer unit normal to $\Omega_{e}$. Notice that $\partial\Omega\subset\partial\Omega_{e}$.

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

[TABLE]

for a positive constant $k_{0}$.

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 $\Omega$, that measures $50\mu m$ on the side (Figure 1). The radius of the circles is within the range of $0.150\;\mu m$ up to $1.141\;\mu m$. The extracellular region is the exterior of the circles and the diffusion coefficient is set to $k_{0}=450\mu m/s^{2}$.

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.

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 $\Omega_{e}$ be partitioned into non overlapping polygonal elements, e.g. a triangulation. Let us denote by $K\equiv K^{-}$ one of such elements, see Figure 2.

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

[TABLE]

Multiplying (6) by the test functions $\mathbf{v}=(v_{1},v_{2})$, and $v$, we obtain after integrating by parts,

[TABLE]

[TABLE]

A second feature of the DG method is the element-wise approximation of the unknown $\mathbf{q}$ and $u$. Continuity is not enforced at the boundary of adjacent elements. Consequently, the boundary terms $u\mathbf{n}$, $k\mathbf{q}\cdot\mathbf{n}$ are replaced by boundary fluxes $h_{u,K^{-}}(u^{-},u^{+},\mathbf{n}^{-})$, $h_{\mathbf{q},K^{-}}(\mathbf{q}^{-},\mathbf{q}^{+},\mathbf{n}^{-})$. As customary, the $-$ superscript denotes limits from the interior of $K^{-}$, and the $+$ superscript limits from the exterior.

This yields in element $K^{-}$

[TABLE]

Let $w$ be an approximation of any of the scalar funcions $q_{1}$, $q_{2}$ ,$u$. Within $K^{-}$ the approximation is given by

[TABLE]

For a triangulation, the functions $N_{j}$ are chosen as in the Finite Element Method with lagrangian interpolation.

It is assumed that $u^{+}$ and $\mathbf{q}^{+}$ are known in (7). Hence, we are led to solve an differential-algebraic system for $u^{-}$ and $\mathbf{q}^{-}$. The solution in time is advanced by a Runge-Kutta method.

**Remark. **We stress that the solution of the $p+1$ differential-algebraic system (7), is solved independently for each element. In practice $p\leq 3$ 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 $u$ a central flux is considered, namely

[TABLE]

A physical assumption is that there is no flow between axons and the extracellular region. Consequently, for $\mathbf{q}$, we propose the numerical flow

[TABLE]

This is coined for the problem under consideration. For instance, if $k^{-}>0$ and $k^{+}=0$, the harmonic mean forces a zero Neumann condition, hence there is no flow trough the boundary of the element $K^{-}$ 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

Also in this MRI application, the heat equation is solved in the extracellular region $\Omega_{e}$ where Cauchy data is prescribed. Apparently, there is no need to consider the axon region $\Omega_{a}$. Nevertheless, we solve the heat equation in elements contained in $\Omega_{a}$ 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 $(x_{i},y_{i})\in\Omega_{e}$, $i=1,2,\ldots,m$ randomly and uniformly in $\Omega_{e}$. Set a final time $T$.

• •

For $i=1,2,\ldots,m$, let $U_{i}(x,y,T)$ be the DG-Cuda solution of the IVBP (1)$-$(4), where the Cauchy data is the Dirac’s delta function supported in $x_{i}$.

• •

For $i=1,2,\ldots,m$, let $u_{i}(x,y,T)$ be obtained from $U_{i}(x,y,T)$ by centering and normalization to yield a density function

• •

Construct the mixture model

[TABLE]

• •

Fit a Gaussian density to $u$, that is, determine a covariance matrix $\Sigma$, such that

[TABLE]

The covariance matrix is given by,

[TABLE]

It is computed by quadrature rules using $u_{ij}$, the values of the numerical solution at the nodes $(x_{i},y_{j})$ in the mesh.

4.2. Numerical results

The substrate under study occupies a square domain of side $50\;\mu m$. Cauchy data is given within a centered box of side $20\;\mu m$. The observation time is $t=0.036$ seconds. At this given time the outer boundary is not reached by diffusion. In all cases a square grid of $K\times L\equiv 400\times 400$ 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 $k_{0}=450\mu m/s^{2}$. Graphical results in Figure 4.

Let us compare the free diffusion (without axons) approximation, versus the analytical solution. The latter is a bivariate Gaussian function $f(x,y;\Sigma)$ defined in (• ‣ 4.1), where the covariance matrix is

[TABLE]

Taking $k=450\mu m/s^{2}$ and $T=0.036s$, non zero coefficients in the covariance matrix are equal to $32.4$.

The DG-CUDA Gaussian matrix fit for $400\times 400$ is

[TABLE]

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 $D=4.5e-10$, $t_{s}=0.036$ duration of the diffusion simulation. $T=5000$ is the number of time steps in the simulation. The step length $l$ is obtained from the relation

[TABLE]

The obtained MC Gaussian matrix is

[TABLE]

To gauge the approximations we compute the least squares residual

[TABLE]

In both cases the approximation of the Gaussian function is highly accurate. The least squares residual (9) is of the order $O(10^{-8})$. 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 $\Omega_{a}$ is defined by 1901 axons. The diffusion coefficient as in (5). The scheme above is applied to a mixture of $m=37$ densities. An instance of one PDE solution is shown in Figure 5. The solution of the full scheme is in Figure 6.

The covarince matrix by the DG-CUDA algorithm is

[TABLE]

Finally, let us summarize execution time in the following table

[TABLE]

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.

Acknowledgements

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.