A computational method for the covariance matrix associated with extracellular diffusivity on disordered models of cylindrical brain axons

1.   Introduction

Diffusion in complex media is a perennial research topic, whose underlying model is a partial differential equation (PDE) of parabolic type. Our motivation comes from the simulation of water diffusion inside the brain tissue (white and grey matter) problem in medical sciences ^[1].

Simulation of the water diffusion inside the brain tissue is motivated by the analysis of the diffusion magnetic resonance images (dMRI). The diffusion simulation is performed using Monte Carlo methods ^[2]. This study provides insightful novel theories on the diffusion properties with analytical representations ^[3] and generates useful data for the validation of model fitting ^[4]. In this problem, the domain is composed of cell bodies (soma, axons, neurites, glial cells, vessels, etc). Despite some of the cell structures can be approximated with simple geometries, the domains are geometrically complex. In particular, the extracellular spaces present arbitrary shapes (similar to a "gruyere-cheese shape"), hence numerical approximations are required to simulate the molecular diffusivity.

The approximation of the extracellular diffusivity profile on a disordered model of cylindrical brain axons, with the diameters computed from a Gamma distribution, is the focus of this study ^[2]. The distribution of diameters has been empirically observed in histological studies on ex vivo samples of the brain ^[5].

The extracellular diffusion process is characterized by the corresponding 3D covariance matrix $\Sigma$ ^[6] of a Gaussian Ensemble Average Propagator, i.e., the Spatio-temporal displacement distribution of the water molecules $P(x, t)$ ^[7]. The $\Sigma$ eigenvalues (matched with the corresponding eigenvectors) indicate the magnitude and orientational dependency, such that, it is possible to infer extracellular microstructure features like the main orientation of the axon bundles, percentage of the volume occupied by neurons (intracellular signal fraction), amount of diffusion anisotropy of the tissue (fractional anisotropy), and axial and radial diffusivity. Those descriptors computed from $\Sigma$ have been correlated with several brain damage and diseases ^[1].

In the Diffusion Weighted Magnetic Resonance (DWMR) medical literature, the covariance matrix is computed using the zeppellin model from the extracellular MR signal, i.e., a limited diffusion tensor model with two equal eigenvalues ^[6,8]. The estimated $\Sigma$ is also used to compute synthetic MR signals when the Short-Gradient-Pulse approximation is assumed ^[7].

The approximation of the covariance matrix is essential to correlate its properties to the microstructural features and their change, respectively, due to, tissue damages associated with brain diseases ^[1]. The covariance matrix can be obtained from a Gaussian density, which is the average density resulting from the solution of diffusion equations. We contribute by solving these diffusion equations using an ad hoc implementation of the Discontinuous Galerkin (DG) method ^[9].

In applications, the numerical solution of the PDE is part of the research process. This is performed in tandem, a model is first proposed, then a numerical method for simulation. This work applies the DG method to accomplish both objectives simultaneously.

Our implementation takes advantage of the underlying physical and geometric properties of the problem as posed by the clinicians. The so-called substrates are squared pixel domains that allow a uniform mesh. Also, we suppose that no diffusivity occurs between the axons and the extracellular region (impermeability) ^[6]. The usual approximation solves the diffusion equation in the extracellular region imposing a zero Neumann condition. The numerical fluxes then model the interaction between the axon walls and the extracellular region. This interaction corresponds to solving the diffusion equation in the whole domain. Thus, null computations are performed in axons. This redundant strategy not only circumvents the computationally expensive task, in Galerkin-type methods, of handling the boundary conditions, but more importantly, it allows the full strength of the DG method, to conduct 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-CUDA implementation is appropriate.

The actual problem of water diffusion inside the brain tissue is in 3D. We highlight the modeling and computational features of the DG method. We present computational results only in 2D, but our theoretical framework will be in 3D. The complex extracellular 2D domains used for our simulations are from the assumption that intracellular spaces (axons) can be modeled as cylinders such that the diffusion parallel to the main axon bundle orientation is unhindered, while the study focuses on the estimation of the diffusion properties perpendicular to the axons ^[10]. The aforementioned assumption is a well-known approximation when the axon pack density is high, as in the brain's corpus callosum region, and employed in various studies ^[3,4]. The outline is as follows:

Section 2 introduces the initial boundary value problem (IVBP) of the diffusion equation associated with the phenomenon. Next, the basics of the DG method are discussed to present the numerical fluxes modeling the interaction between the axons and the intercellular region. Then, the computational model for GPU-CUDA implementation is described. We complete the methodology with the scheme for computing the covariance matrix. Section 3 provides two numerical examples of the proposed method. First, the analytical solution of a free diffusion problem is considered to gauge the introduced DG-CUDA method and a classical Monte Carlo approach. Second, a benchmark type substrate associated with an ex vivo experiment is introduced to test the computational methodology. A computational comparison using the execution time is also described. We close with conclusions and a brief discussion on future work in Section 4.

2.   Materials and methods

2.1.   Water diffusion inside the brain tissue

We model a substrate occupying a domain (medium) comprised of two regions, the axons, and extracellular complement. The former is nondiffusive and the latter a diffusive region. We assume that the boundary between both regions is reflecting. Initial pulses are prescribed in the extracellular region, far away from the outer boundary which is modeled 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$ representing the axon and extracellular regions, respectively.

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

given Cauchy data

and boundary values

Here $\mathbf{n}$ is the outer unit normal to $\Omega_e$.

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

The discontinuous diffusion is given by

for a positive constant $k_0$.

2.2.   Numerical Flux and the Discontinuous Galerkin method

2.2.1.   The Discontinuous Galerkin method

Let $\Omega_e$ be partitioned into nonoverlapping polygonal elements, e.g. a triangulation. Denote by $K\equiv K^-$, one of such elements, see Figure 1.

A striking characteristic of the DG method is that it reduces the PDE to a first-order system. Consequently, consider

Multiplying (2.6) by the test functions $\mathbf{v} = (v_1, v_2)$, and $v$, we obtain after integrating by parts,

Continuity is not enforced at the boundary of adjacent elements. Therefore, the boundary terms $u$, $k \mathbf{q}$ are replaced with boundary numerical fluxes $u^*\equiv u^*(u^-, u^+)$, $(k\mathbf{q})^*\equiv(k\mathbf{q})^*(k^-, \mathbf{q}^-, k+, \mathbf{q}^+)$. As customary, the $-$ superscript denotes limits from the interior of $K^-$, and the $+$ superscript, limits from the exterior.

This yields in element $K^-$

These numerical fluxes chosen below model the underlying physical phenomena.

Another characteristic of the DG method, is the element-wise approximation of the unknown $\mathbf{q}$ and $u$.

Suppose $w$ is an approximation of any scalar functions $q_1$, $q_2$, $u$. Then, the approximation within $K^-$ is given by

For a triangulation, choosing the functions $N_j$ follows the Finite Element Method with Lagrangian interpolation.

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

Remark. The solution of the $p+1$ differential-algebraic system (2.7), is solved independently for each element. In practice $p\leq 3$ suffices. Thus, we have small systems, that do not exchange information on each time step, to solve.

2.2.2.   The numerical fluxes for water diffusion

No preferred direction of propagation in diffusion phenomena exists, thus for $u$, a central flux is considered, namely,

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

This is coined for the problem under consideration. 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 required.

2.3.   CUDA implementation

The meshing is time-consuming in Galerkin-type methods, as it is boundary conditions handling when assembling the local systems. Here, the domain is divided into square pixels. For meshing, we consider the simplest triangular mesh by separating these squares through the diagonals.

Also in this MRI application, the heat equation is solved in the extracellular region $\Omega_e$ with Cauchy conditions. No consideration is given to 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 balance computations on every element. From the remark above, the calculations in (2.7) are element-independent when updating the time. Consequently, the main ingredients for parallel processing using GPUs are met and include small balanced computations with no exchange of information between processors.

2.4.   The covariance matrix scheme

Let $(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 solution of the IVBP (2.1)$-$(2.4), where the Cauchy condition is the Dirac delta function supported in $(x_i, y_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 given as

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

The covariance matrix is given by,

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.

3.   Numerical results

3.1.   Free diffusion

In this example, we solve the heat equation with the Cauchy condition as the Dirac delta supported at the origin. The uniform diffusion coefficient is $k_0 = 450 \mu m/s^2$.

We compare the free diffusion approximation (without axons), versus the analytical solution. The latter is a bivariate Gaussian function $f(x, y;\Sigma)$ defined in (2.4), and the covariance matrix is

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

Graphical results are shown in Figure 2. The DG-CUDA Gaussian matrix fit for $400\times400$ is

An alternative construction of the covariance matrix is using Monte Carlo diffusion simulation ^[2]. We give the list of 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 using the relation

The obtained MC Gaussian matrix is

Remark. To gauge the approximations we compute the least squares residual

In both cases, the approximation of the Gaussian function is accurate. The least squares residual (3.1) is of the order $O(10^{-8})$. In practice, the Gaussian density functions coincide. However, the DG-CUDA approximation is structurally more consistent. The matrix is symmetric, the diagonal values coincide and the other terms are two orders of magnitude smaller than those obtained from MC.

3.2.   A case study

We consider a substrate consisting of 1901 nonoverlapping circles, which represent the axons, and consider the regions only within a square, the domain $\Omega$, that measures $50 \mu m$ on the side (Figure 3). The radius of the circles ranges from $0.150 \; \mu m$ 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 substance of the medium ^[11]. The smaller substrate is used.

The substrate under study occupies a square domain of side $50 \; \mu m$. Cauchy condition is given within a centered box of side $20 \; \mu m$ and observation time $t = 0.036s$. 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.

The axon region $\Omega_a$ is defined by 1901 axons. The diffusion coefficient as in (2.5). The scheme above is applied to a mixture of $m = 37$ densities. We show one PDE solution in Figure 4.

The covariance matrix by the DG-CUDA algorithm is

and the corresponden one obtained MC Gaussian matrix is

where the relative error between entries is smaller than 0.5%.

Remark. The choice of the sample size, $m$, for the Gaussian mixture is heuristic. From our experiments, $m$ need not be large. In two dimensions $O(10)$ will suffice.

3.3.   Computational efficiency

We summarize the execution time in the following table

Regardless of the diffusion coefficient, the GPU processing is more than 1000 times faster for the solution of the heat equation with a single delta as Cauchy condition.

For the prototype computation, we implement parallel processing using GPUs in the CUDA language of a DELL laptop with hardware:

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

GPU: NVIDIA Corporation GM107GLM [Quadro M1000M]

4.   Conclusions

This study proposed a methodology for the efficient computation of the covariance matrix associated with the extracellular diffusivity profile on a disordered model of cylindrical brain axons. The methodology relies on two strategies; introducing physical numerical fluxes in the Discontinuous Galerkin method and null computations on the axon region for straightforward GPU implementation.

We present satisfactory results on two-dimensional examples that show the proposed methodology outperforms sequential implementations and is more structurally consistent than a Monte Carlo approach.

Presently, we are working toward extending our methodology to the three-dimensional case. We contend that our methodology can include the corresponding complexities in a 3D extension. We mention that the choice of the sample size, $m$, for the Gaussian mixture is heuristic and as such, an in-depth analysis is required. These issues make up our current and future work.

Acknowledgments

M. A. Moreles acknowledges partial support from the project CONACYT A1-S-17634. The authors would like to thank Enago (www.enago.com) for the English language review.

Conflict of interest

All authors declare no conflicts of interest in this paper.