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 DGGPU 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
(1) 
given Cauchy data
(2) 
and boundary values
(3) 
(4) 
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
(5) 
for a positive constant .
2.2. A case study
The substrate under consideration consists of 1901 nonoverlapping 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.
It is apparent that numerical approximations are required to simulate the molecular diffusivity is this rather complex porous medium.
3. DGCUDA 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.
A defining feature of the DG method is to reduce the PDE to a first order system. Consequently, let us consider
(6) 
Multiplying (6) by the test functions , and , we obtain after integrating by parts,
A second feature of the DG method is the elementwise 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
(7) 
Let be an approximation of any of the scalar funcions , ,. Within the approximation is given by
(8) 
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 differentialalgebraic system for and . The solution in time is advanced by a RungeKutta method.
Remark. We stress that the solution of the differentialalgebraic 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
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) i76820HQ CPU @ 2.70 GHz
GPU: NVIDIA Corporation GM107GLM [Quadro M1000M]
4. The Gaussian profile of the extracellular diffusivity
In this section we introduce a scheme to approximate the covariance matrix of the EAP. Hence describing the Gaussian profile of the extracellular 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 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.
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 DGCUDA 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
(9) 
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 DGCUDA 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.
The covarince matrix by the DGCUDA 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.Acknowledgements
D. Cervantes and M. A. Moreles thank the support of ECOSNORD through the project: 000000000263116/M15M01. M. A. Moreles also acknowledges partial support from the project CIB CONACYT 180723.
References
 [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 timedependent diffusion. NeuroImage, 114:18–37, 2015.
 [5] Bernardo Cockburn, George E Karniadakis, and ChiWang 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] NedjatiGilani 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., StaitGardner 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. WheelerKingshott, 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 twocompartment 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 montecarlo 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, JingRebecca Li, Denis Grebenkov, and Denis Le Bihan. A finite elements method to solve the blochtorrey 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. Pulsedfield gradient nuclear magnetic resonance as a tool for studying translational diffusion: Part 1. basic theory. GAnimal’s Journal, 1997.
 [19] Alonso RamirezManzanares, Mario OcampoPineda, Jonathan RafaelPatino, Giorgio Innocenti, JeanPhilippe Thiran, and Alessandro Daducci. Quantifying diameter overestimation of undulating axons from synthetic dwmri. 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.
Comments
There are no comments yet.