1 Introduction

The glymphatic system (GS) is the structural entity whereby waste products are transported from the brain and into lymphatic vessels located outside, in the meninges and along the neck vasculature [8]. Importantly, the GS also flushes out of the brain soluble amyloid beta (A\(\beta \)) and tau proteins, the main culprits of Alzheimer’s disease (AD) in humans and animals [5]. Despite the potential implications of the GS for AD and other neurodegenerative conditions, there are significant gaps in our understanding of the waste clearance mechanisms and the physical forces controlling transport.

Glymphatic transport behavior can be observed with a temporal series of contrast enhanced MR images of the rodent brain. Briefly, the small molecular weight gadolinium (Gd) contrast agent (tracer) is infused into the cerebrospinal fluid (CSF) reservoir of the cisterna magna and its spatial distribution into the brain is captured by the successive acquisition of 3D T1-weighted images (234 \(\mu \)m resolution) every \(\sim \!\!4\) minutes, for a total of \(\sim \!\!3\) hours [6]. However, these MRIs do not provide directional information on the tracer movement between time points. Therefore, there is an urgent need for a mathematical framework that can capture and help visualize the dynamic tracer behavior in a manner aligning well with the biological understanding.

In this work, we present a novel visualization framework, GlymphVIS, for studying glymphatic transport pathways using regularized optimal transport (OT). The theory of OT seeks the most feasible way to redistribute mass from one given distribution to another while minimizing the associated cost of transportation ([9, 12]). OT has been used for registration and connectivity analysis of brain white matter [7], image morphing [4], and has recently been extended to the case of measures of different total mass [2].

Fig. 1.
figure 1

GlymphVIS Pipeline: (a) The generalized regularized OT procedure (GR-OT) takes initial and final ‘observed’ density images as input and returns the ‘clean’ or ‘believed true’ final density image along with the corresponding velocity vector field describing the deformation. (b) The output density images and velocity are then subsequently passed to the flow pattern analysis procedure (FPA) which returns pathway and streamline clustering visualization for the whole time domain.

Full size image

Ratner et al. [10] modeled the glymphatic flow using the traditional OT formulation. This approach yielded promising results, but at the expense of some unrealistic assumptions: (1) the movement of the contrast agent is not affected by diffusion, (2) the total mass of the tracer remains constant over time, and (3) the given MRIs represent the true density distribution at that time. From the implementation perspective, the mass conservation constraint requires normalizing the density distribution which can be drastically altered by the presence of noise in the data and additional noise interference is caused by taking the given images as fixed endpoints. Finally, the authors of [10] do not explicitly model time which means the resulting deformation field cannot reflect time-varying dynamics. In this work, we introduce a new and more physiologically relevant model inspired by the work of Benamou and Brenier [1] that relaxes the above-mentioned unrealistic constraints. Specifically, the contributions of this paper are enumerated as follows:

  1. 1.

    We replace the continuity equation with the advection-diffusion equation to more accurately model the flow behavior and smoothen the deformation field;

  2. 2.

    We no longer enforce total mass conservation, so normalization of the density distributions is not needed;

  3. 3.

    We treat the final time condition as a free endpoint which prevents overfitting to noise;

  4. 4.

    We explicitly model the time domain, which allows for a direct temporal analysis of the dynamic flow behavior.

This project was supported by AFOSR grant FA9550-17-1-0435), ARO grant (W911NF-17-1-049), grants from National Institutes of Health (1U24CA18092401A1, R01-AG048769), MSK Cancer Center Support Grant/Core Grant (P30 CA008748), and a grant from Breast Cancer Research Foundation (grant BCRF-17-193).

2 GlymphVIS

Benamou and Brenier [1] recast the OT problem in the context of fluid mechanics that explicitly yields a time-interpolant between the two densities. This naturally motivates an ideal framework for studying the glymphatic pathways because it allows for more direct control and variation in modeling its dynamic flow behavior. Here, we introduce the following two terms to the original Benamou and Brenier OT formulation: (1) a regularization term to alleviate the effect of noise and (2) a diffusion term in the standard continuity equation to better model both advection and diffusion in the glymphatic system. We then clusters the streamlines from the resulting velocity field in order to elucidate and visualize the conduits of glymphatic flow and efflux; see Fig. 1.

2.1 Regularized OT

In order to motivate our model formulation, we begin with our assumptions about the data:

Assumption 1

Image intensity is proportional to tracer mass, (and we therefore refer to the intensity as mass).

Assumption 2

Tracer is transported via glymphatic pathway, as supported by experimental findings [5].

Assumption 3

Apparent motion of glymphatic transport is governed by the advection-diffusion equation (ADE),

(1)

where \(0\le \rho :[0,T]\times \mathcal {D}\rightarrow \mathbb {R}\) is a density with compact support, \(\mathcal {D}\subset \mathbb {R}^d,\) \(v:[0,T]\times \mathcal {D}\rightarrow \mathbb {R}^d\) is velocity and \(0 \le \sigma \in \mathbb {R}\) is diffusivity.

Assumption 4

The MR images we are given are noisy observations of the tracer’s conditions at time \(t=T_i\),

$$\begin{aligned} \rho (T_i,x) + \epsilon = \rho _{T_1}^\mathrm{obs}(x), \quad i = 0,...,N, \end{aligned}$$
(2)

where \(\epsilon \) is a random Gaussian iid with covariance \(\varSigma \).

Given initial and final observations of tracer density \(\rho _0^\mathrm{obs}\) and \(\rho _T^\mathrm{obs}\) at times \(t=0\) and \(t=T\) respectively, our goal is to find the velocity field v and the ‘believed true’ or ‘clean’ image \(\rho _T\) such that the constraint (1) is satisfied. To this end, we propose to minimize the objective function

$$\begin{aligned} \mathcal{J}[\rho ,v] = \int _0^T\int _{\mathcal {D}}\, \frac{1}{2}\rho \Vert v\Vert ^2\, dxdt +\alpha \Vert \rho (T,x) - \rho _T^\mathrm{obs}(x)\Vert _{\varSigma }^2, \end{aligned}$$
(3)

subject to the ADE constraint (1) with the initial condition \(\rho (0,x) = \rho _0^\mathrm{obs}(x)\). Here, we have added the second term to the Benamou and Brenier energy functional [1] so that noise, which is inherent in all sensor-derived data, is explicitly taken into account. The parameter \(\alpha \) weighs the balance between fitting the data and minimizing the energy associated with transporting the mass. We refer to (3) as the generalized regularized OT problem (GR-OT) and note that supplemental regularization can easily be implemented by adding intermediate densities to help guide the optimization procedure toward more accurate results.

Remark:

In the original Benamou-Brenier formulation of OT, \(\sigma =0\) and endpoint distribution is specified in Eq. (1), and \(\alpha =0\) in Eq. (3).

Proposed Minimization Method. While it is possible to solve the optimization problem (3) as a constrained one, it is straightforward to eliminate the ADE (1) and solve the problem for v alone. Accordingly, consider solving the PDE for a given v, obtaining the smooth, differentiable map

$$\begin{aligned} F(v) = \rho (t,x), \quad t\in [0,1]. \end{aligned}$$
(4)

Discretization. Suppose the given images are \((n_1\times \ldots \times n_d)\) in size, let \(s=n_1 * \ldots * n_d\) denote the total number of voxels, and let m denote the number of time steps such that \(m*\delta t = T\). We will use bold font to denote linearized variables.

We use mimetic methods, designed to keep their properties when considering inner products, for the discretization of the problem. First, we use operator splitting to discretize the ADE (1) as an advection step and diffusion step, independently. In the first step, we consider the advection equation and solve the problem

$$ \frac{\partial \rho }{\partial t} + \nabla \cdot \,(\rho v) = 0 \quad \rho (t_n,x) = \rho _n. $$

Using a particle in cell (PIC) method, we obtain the discrete equivalent \({\varvec{\rho }}_{n+1}^* = \mathbf{S}(\mathbf{v}_n) {\varvec{\rho }}_n\) where \(\mathbf{S}\) is a linear interpolation matrix. The method is conservative which means that no mass is lost during this step. For the second step, we consider the diffusion equation and solve the problem

Using the backward Euler method, we obtain the discrete equivalent

$$\begin{aligned} (\mathbf{I}- \delta t \mathbf{A}) {\varvec{\rho }}_{n+1} = {\varvec{\rho }}_{n+1}^* \end{aligned}$$
(5)

where \(\mathbf{I}\) is the identity matrix and \(\mathbf{A}\) is a discretization of the diffusion operator on a cell centered grid. Combining these two steps, we obtain the corresponding discrete forward problem \((\mathbf{I}- \delta t \mathbf{A}) {\varvec{\rho }}_{n+1} = \mathbf{S}(\mathbf{v}_n) {\varvec{\rho }}_n, n=0,\ldots ,m.\) Clearly, the density at any time step depends only on the initial density \({\varvec{\rho }}_0\) and the velocity \(\mathbf{v}\), allowing us to define the discrete map \(F(\mathbf{v})\) (4) that maps the velocity to the density at all times. Next, defining \({\varvec{\rho }}= [{\varvec{\rho }}_1^{\top },\ldots ,{\varvec{\rho }}_{m+1}^{\top }]^{\top }\) and \(\mathbf{v}= [\mathbf{v}_0^{\top },\ldots ,\mathbf{v}_m^{\top }]^{\top }\), a straightforward discretization of the energy yields

$$\begin{aligned} \int _0^T \int _{\mathcal {D}} \rho \Vert v\Vert ^2 dx\, dt \approx h^d \delta t{\varvec{\rho }}^\top (\mathbf{I}_m \otimes \mathbf{A}_v) (\mathbf{v}\odot \mathbf{v}), \end{aligned}$$
(6)

where \(h^d\) is the volume of each cell, \(\mathbf{I}_k\) is the \(k \times k\) identity matrix, \(\mathbf{A}_v\) is a \(1 \times d\) block matrix of \(\mathbf{I}_s\), \(\otimes \) denotes the Kronecker product and \(\odot \) denotes the Hadamard product. We then solve the discrete optimization problem which now reads

$$\begin{aligned} \min&\phi ( \mathbf{v}) = \frac{1}{2}h^d\, \delta t {\varvec{\rho }}^\top (\mathbf{I}_m \otimes \mathbf{A}_v) (\mathbf{v}\odot \mathbf{v}) + \alpha \Vert {\varvec{\rho }}_n - {\varvec{\rho }}_n^\mathrm{obs}\Vert ^2 \\ \nonumber \mathrm{subject \ to}&{\left\{ \begin{array}{ll} (\mathbf{I}- \delta t \mathbf{A}) {\varvec{\rho }}_{n+1} - \mathbf{S}(\mathbf{v}_n) {\varvec{\rho }}_n = 0 \\ {\varvec{\rho }}_0 = {\varvec{\rho }}_0^\mathrm{obs}. \end{array}\right. } \end{aligned}$$
(7)

Note that the objective function is quadratic with respect to \(\mathbf{v}\) and the interpolation matrix \(\mathbf{S}\) is linear with respect to \(\mathbf{v}\) as it contains the weights on the linear interpolation. Following [11], one can use a Gauss-Newton like method to solve the problem.

2.2 Flow Pattern Analysis

The time-interpolant of density images and corresponding time-varying velocity vector field \(\mathbf{v}(t,x)\) directly output by the GR-OT procedure is then fed into our flow pattern analysis procedure (FPA). For each time step, we construct streamlines by integrating the velocity field \(\mathbf{v}\). By looking at the streamline density through each voxel, we get a global visualization of the GS ‘pathways’. In order to supplement this with local information, we cluster the streamlines using the QuickBundles algorithm [3]. Significant clusters provide more information regarding different flow trajectories within different pathways and fluid reservoirs. Both the pathways and clusters are converted to NIfTI files where they are analyzed by overlaying anatomical masks using Amira software specifically designed for visualization of data in 3D and 4D. We discuss these results in the following section.

3 Results

In order to quantitatively assess the performance of our model, we look at the registration error between the model returned ‘clean’ density and the target image density. Taking the mean square of the error and the infinity norm of the error, our model (with no diffusion, i.e. \(\sigma =0\)) yields up to 5 times smaller errors than the traditional OT model proposed in [10] (Fig. 2). This large improvement was possible due to the aforementioned adjustments made to account for noise in the data. We then introduce a little diffusion (\(\sigma = 0.002\)) and look at the root mean square error between the returned ‘clean’ densities with the ‘clean’ densities obtained by increasing the diffusion parameter by factors of 10 (\(\sigma = 0.02, 0.2\)). The robustness of the diffusion parameter is shown by these errors, given in Fig. 3c, as well as by the consistent pathways found with multiple values of \(\sigma \), see Figs. 3a and b.

Fig. 2.
figure 2

Registration error between model returned final density and target image density for the traditional OT model [10] (shown in blue) and our GR-OT model (shown in red).

Full size image
Fig. 3.
figure 3

Robustness of diffusion parameter. (a) Pathways obtained with \(\sigma = 0.002\) and (b) pathways obtained with \(\sigma = 0.2\). (c) Root mean square error between ‘clean’ densities obtained with \(\sigma = 0.002\) and \(\sigma = \tilde{\sigma }\). Top row: \(\tilde{\sigma } = 0.02\). Bottom row: \(\tilde{\sigma } = 0.2\).

Full size image
Fig. 4.
figure 4

GlymphVIS pathways. (a) Original contrast enhanced MRI highlighting the MCA area. (b) 3D volume rendering of the GlymphVIS pathways in relation to the whole rat brain (grey scale, volume rendered) demonstrating that GlymphVIS pathways track CSF transport along the MCA from the level of the Circle of Willis (CW) to where it crosses the olfactory tract (not shown) and proceeds dorsally onto the surface of the brain. (c) GlymphVIS pathways without the whole brain. Details of the pathways in other areas are now visible including pathway reservoirs associated with the basal cistern, the interpenduncular cistern (IpC) and cleft between the hippocampus and other brain nuclei.

Full size image
Fig. 5.
figure 5

GlymphVIS clusters. (A,B) Clusters shown in different colors (Olf=olfactory bulb and Cb=cerebellum). (C) Anatomical MRI from the ventral surface, where the MCA and internal carotid artery (ICA) can be visualized as single, vascular structures running along the surface of the brain. In addition, the acoustic nerve and inner ear complex (cochlea) is included. (D) Selected streamline clusters related to the MCA and the cochlea overlaid on anatomical template. (E) Selected streamline clusters related to the MCA and the cochlea.

Full size image

The utility of GlymphVIS is further validated by its success in reproducing known aspects of glymphatic transport. This is illustrated by the pathways and clusters derived from the MRIs at 1.2hr after contrast infusion into the CSF, shown respectively in Figs. 4 and 5. In particular, Fig. 4 demonstrates that pathways found by our methodology have accurately captured glymphatic peri-arterial transport along the MCA and in other areas such as the CSF reservoirs. Even more promising, are the trajectories shown by the streamline clusters in Fig. 5. This is the first time that specific contrast relevant streamlines have been captured moving towards the inner ear, and illustrates the promise of GlymphVIS and the new GR-OT flow analysis pipeline.

4 Conclusions and Future Work

In this paper, we considered a modification of the Benamou-Brenier formulation of OT in which both the continuity and energy cost functionals were modified. This was done to take into account noise as well as possible diffusion in the glymphatic flows for “normal” rat brains. In the future, we also intend to consider cases in which there may be some pathologies, in particular, rat brain models in which there is evidence of AD and vascular dementias. The concept and hypothesis to be tested would be to see if using these mathematical techniques, one could quantitatively differentiate between normal and aberrant CSF flow inside as well as outside the brain, which specifically relate to evolving neuropathology. Finally, one can consider the technique we have proposed as one of deformable registration. In contrast to other deformable methods such as LDDMM, we are not constrained by only considering diffeomorphic transformations. Moreover, in our setting, we have explicitly taken into account the advection-diffusion nature of the flow, and thus the underlying physics.