Fisher-Rao Regularized Transport Analysis of the Glymphatic System and Waste Drainage

TL;DR

This paper introduces a Fisher-Rao regularized Lagrangian framework for analyzing the glymphatic system's waste drainage, enabling visualization and differentiation of flow patterns under various conditions, with implications for understanding Alzheimer's disease.

Contribution

It presents a novel Lagrangian approach using Fisher-Rao regularization for analyzing time-varying flows in the glymphatic system, offering interpretability and visualization advantages over existing methods.

Findings

Successfully captures known glymphatic flows.

Distinguishes flow patterns under different anesthetics.

Provides insights into altered waste drainage states.

Abstract

In this work, a unified representation of all the time-varying dynamics is accomplished with a Lagrangian framework for analyzing Fisher-Rao regularized dynamical optimal mass transport (OMT) derived flows. While formally equivalent to the Eulerian based Schr\"odinger bridge OMT regularization scheme, the Fisher-Rao approach allows a simple and interpretable methodology for studying the flows of interest in the present work. The advantage of the proposed Lagrangian technique is that the time-varying particle trajectories and attributes are displayed in a single visualization. This provides a natural capability to identify and distinguish flows under different conditions. The Lagrangian analysis applied to the glymphatic system (brain waste removal pathway associated with Alzheimer's Disease) successfully captures known flows and distinguishes between flow patterns under two different anesthetics, providing deeper insights into altered states of waste drainage.

Equations

$$\underset{\mu ,\nu }{\inf }{\int }_0^1{\int }_{\mathcal{R}}$$

$$\mathrm{subject}\mathrm{to}$$

$$\hat{\nu }(t,x)=\nu (t,x)-ϵ\nabla \log \mu (t,x).$$

$$\nabla \cdot (\mu \hat{\nu })=\nabla \cdot \left[(\nu -ϵ\nabla \log \mu )\mu \right]=\nabla \cdot (\mu \nu )-ϵ\Delta \mu ,$$

$${\partial }_t\mu +\nabla \cdot (\mu \hat{\nu })=0.$$

$$X(0,x)=x,{\partial }_{t}X={\hat{\nu }}_{min}(t,X(t,x)),$$

$$\underset{\mu ,\hat{\nu }}{\inf }{\int }_0^1{\int }_{\mathcal{R}}$$

$$\mathrm{subject}\mathrm{to}:$$

$$H({\mu }_1\Vert {\mu }_0):={\int }_{\mathcal{R}}{\mu }_1(x)\log {\mu }_1(x)-{\mu }_0(x)\log {\mu }_0(x)dx.$$

$$\hat{\nu }(t,x)+ϵ\nabla \log \mu (t,x)=\nu (x,t)$$

$${\int }_0^1{\int }_{\mathcal{R}}\mu \Vert \hat{\nu }+ϵ\nabla \log \mu {\Vert }^{2}dxdt={\int }_0^1{\int }_{\mathcal{R}}\mu \Vert \hat{\nu }{\Vert }^2+2ϵ\mu ⟨\hat{\nu },\nabla \log \mu ⟩+{ϵ}^2\mu \Vert \nabla \log \mu {\Vert }^{2}dxdt,$$

$$J:={\int }_0^1{\int }_{\mathcal{R}}⟨\mu \hat{\nu },\nabla \log \mu ⟩dxdt=H({\mu }_1\Vert {\mu }_0)$$

$$J=-{\int }_0^1{\int }_{\mathcal{R}}\nabla \cdot (\mu \hat{\nu })\log \mu dxdt=-{\int }_0^1{\int }_{\mathcal{R}}\mu {\partial }_t(\log \mu )dxdt+H({\mu }_1\Vert {\mu }_0).$$

$${\int }_0^1{\int }_{\mathcal{R}}\mu {\partial }_t(\log \mu )dxdt={\int }_0^1{\int }_{\mathcal{R}}\mu \frac{{\partial }_t\mu }{\mu }dxdt={\int }_{\mathcal{R}}({\mu }_1(x)-{\mu }_0(x))dx=0$$

Full text

1 Introduction

The optimal mass transport problem (OMT) entails minimizing the transportation cost associated with redistributing one probability distribution to match another. The static OMT formulation has been applied to medical image processing and registration [9, 11, 12]. The dynamical version of OMT [1] opened up new possibilities for numerical solutions and extensions using tools from the field of fluid dynamics [3, 21]. In particular, the authors consider a regularization of the dynamical OMT formulation [1] for modeling apparent fluid flow in dynamic contrast enhanced MRI (DCE-MRI) images where a diffusion term has been added to the Euler equation [7, 8]. Contrast transport derived using earlier models without a diffusion term were vulnerable to noise in the data [23], making the inclusion of a diffusion term advantageous from both physical and numerical standpoints. Analysis of the resulting velocity fields is typically carried out in an Eulerian framework, i.e., flow properties are considered at specific locations for each time point as compared to the Lagrangian framework, which tracks specific particles as they move over time; see Figure 1.

The Lagrangian approach gives a useful alternative to the Eulerian scheme and has recently been presented in conjunction with the aforementioned regularized version of OMT [16]. The pipeline taken from this work has untapped practical applications for distinguishing elusive physiologically relevant flow patterns and so we adopt it in the present work to study the glymphatic system (GS) and brain waste drainage. In [16], the authors employ a regularized version of OMT by adding a diffusion term to the continuity equation in the Benamou-Brenier formulation of OMT [1]. It turns out that there are two equivalent ways to derive the Lagrangian coordinates of the optimal trajectories. The first is via a transformation of the advection-diffusion equation in the regularized problem, and the second by adding a Fisher-Rao information type term to the kinetic energy in the original dynamic formulation [4]. In the present work, we adopt the second method, i.e., Fisher-Rao regularization. We should note that employing Lagrangian coordinates derives improved visualizations of GS function in a time efficient display and also reveals new disparate dynamic waste drainage features in two different states of arousal that were previously difficult to discern.

Waste products are removed from the brain through the GS, a peri-vascular transit passage for cerebrospinal fluid (CSF) which facilitates mixing of CSF with interstitial fluid [14]. Impaired waste clearance due to glymphatic dysfunction has been associated with vascular dementia, Alzheimer’s Disease (AD) [22] and sleep deprivation [24]. The GS accelerates waste clearance from the brain during slow wave sleep and with certain hypnotics when compared to other anesthetics and wakefulness [28]. A recent study showed that unconsciousness induced with dexmedetomidine and low dose isoflurane (DEXM+ISO) that mimics natural sleep enhances GS function to a greater extent when compared to deep anesthesia induced with isoflurane (ISO) [2]. These results suggest that hypnotic drugs that promote ‘natural’ sleep might be superior to deeper states of sleep/anesthesia regarding preservation of sleep like GS function. Novel ways to improve or maintain GS function for general brain health is urgently needed and dependent on developing robust analysis tools to quantify its function.

The main contributions of this work are as follows:

We directly use the Fisher-Rao regularization of the energy functional in visualizing the flows, highlighting the framework’s natural treatment of small diffusion values needed for studying the GS. 2. 2.

We successfully differentiate between GS flow patterns under two different anesthetics and provide additional new insights, in particular, to the relationship between GS flow and solute clearance via pathways outside of the brain.

2 Background

2.1 Regularized OMT

As alluded to above, the problem of optimal mass transport (OMT) is concerned with moving mass from one site to another so that minimal ‘work’ is expended, and mass is preserved, see [26, 27] for a complete set of references. For our purposes with respect to medical images, we are interested in the dynamical OMT formulation due to Benamou and Brenier [1]. They proposed an alternative numerical solution by introducing time and solving a partial differential equation constrained space-time minimization problem. This formulation has the geometric interpretation of finding geodesics [20] between the given densities $\mu_{0}$ and $\mu_{1}$ in the space of probability densities and the trajectory of the transport is explicitly factored into the cost.

Explicit representation of the density’s evolution suggests exciting capabilities for improved image registration techniques to account for dynamical aspects of physiological processes with additional aptitude for analyzing interesting time-varying phenomena. This is particularly powerful for medical applications where a complete physical model is often impractical to implement. We consider the following modified version of Benamou and Brenier’s OMT formulation where a diffusion term is added in the continuity equation. The following regularized OMT problem over the normalized time interval $t\in[0,1]$ is used to motivate the Lagrangian framework for flow representation [16] that will be detailed in the next section:

[TABLE]

where $\mu=\mu(t,x)\in\mathbb{R}$ is the density interpolant between the given densities $\mu_{0}$ and $\mu_{1}$, assumed to have the same total mass over a bounded, connected subspace $\mathcal{R}$ of $\mathbb{R}^{d}$, $\nu=\nu(t,x)\in\mathbb{R}^{d}$ is the velocity and $\epsilon$ is the diffusivity. It is interesting to note that this may also be regarded as a reformulation of the Schrödinger bridge problem [4].

3 Lagrangian coordinates

Following [16], the optimal trajectory in Lagrangian coordinates for the regularized case is defined as follows. Let $\mu_{min}$ and $\nu_{min}$ denote the minimum arguments of the action (1) subject to (3). Define the augmented velocity

[TABLE]

Noticing that

[TABLE]

leads to the following conservation form of the constraint (3):

[TABLE]

Analogous to [1], the Lagrangian coordinates of the flow $X=X(t,x)$ corresponding to the minimizing velocity $\nu_{min}$ is the solution of the differential equation

[TABLE]

where according to (4), $\hat{\nu}_{min}=\nu_{min}-\epsilon\nabla\log\mu_{min}.$

3.1 Fisher-Rao regularization

It is very important to note that the regularized OMT problem is equivalent to a Fisher-Rao information theoretic regularization of OMT, which is very closely connected to the Sinkhorn approach [6], and is in fact a dynamic formulation of the Schrödinger bridge problem [4]. For small $\epsilon$ in the advection-diffusion equation above, Sinkhorn will become unstable, which is why one needs to employ a different approach in the medical imaging realm. We will therefore derive an equivalent Lagrangian formulation of regularized OMT via the Fisher-Rao regularized functional. We should note that there are a number of works demonstrating the well-known equivalence of Fisher-Rao and OMT joint interpolation functional and regularized OMT including [4, 18, 5, 19]. We briefly sketch the proof from these works to highlight some of the key steps of going from the regularized OMT problem to the Fisher-Rao regularized one. We refer the reader to the latter works for more details and insights into the methodology.

Following the notation of the previous section and references [4, 18, 5, 19], we claim that the problem defined by (1–3) is equivalent to

[TABLE]

for time-varying densities $\mu=\mu(t,x)\in\mathbb{R}$ and velocities $\nu=\nu(t,x)\in\mathbb{R}^{d}$ where

[TABLE]

Clearly, $H(\mu_{1}\|\mu_{0})$ is a constant. Next, noting that

[TABLE]

and

[TABLE]

we just need to show that

[TABLE]

to verify our claim. However, using integration by parts and the constraints (9) we see that

[TABLE]

Finally,

[TABLE]

since the total masses of $\mu_{1}$ and $\mu_{0}$ are equal.

Letting $\hat{\nu}_{min}$ denote the optimal vector field for (8-9), the Lagrangian coordinate solution is then exactly that given by (7) above. As previously noted, there are very fast algorithms for solving such entropic regularized versions of OMT for sufficiently large $\epsilon$ [6]. Since we need to consider small $\epsilon$ values, we directly solve (8-9) based on the numerical scheme given in [25] where the initial and final densities $\mu_{0},\mu_{1}$ are taken from the given DCE-MRIs.

3.2 Pathlines

In order to represent the characteristics of the flow over all time in one comprehensive figure, we use the Lagrangian framework to construct what are commonly known as pathlines. A pathline $X(t,x)$ with initial position $X(0,x)=x$ is given by (7) and traces the trajectory of an individual particle through-out the time interval. As detailed in [16], additional information such as $\%$-signal from baseline is used to determine regions where flow is more likely to occur. Throughout this region of interest, $p$ seeding points, denoted $x_{l}$, are selected uniformly. Pathlines are then computed by integrating the augmented velocity (4) with initial positions $x_{l}$. The speed $\|\nu(t,X(t,x_{l})\|$ of each particle is simultaneously computed along each pathline, referred to as the speed pathline and denoted $X_{\nu}=X_{\nu}(t,x_{l})$, by (tri)linearly interpolating the derived velocity field.

Next, in order to extract pathlines that are representative of the flow behavior in specific anatomical regions, we subsequently cluster the pathlines by proximity using the QuickBundles algorithm [10]. All computations were implemented in MATLAB and the Lagrangian analysis took approximately 1 minute to run for each of the 8 datasets (5 DEXM+ISO, 3 ISO) of resolution $128\times 128\times 128$.

4 Discussion

While streamlines and pathlines are interchangeable for steady flows (i.e. time independent), the same cannot be said for unsteady flows (i.e. time dependent). A novel aspect of the dynamical OMT formulation [1] is its explicit description of the density’s time-varying evolution, suggesting the need for Lagrangian analysis of this behavior. The Lagrangian approach provides an elegant framework to observe various attributes such as speed and density along the pathlines. Having a single representation for the history of a particle and its attributes across time lends itself to a natural means for differentiating between flows under different conditions. We show the capability of this framework as it pertains to the GS and waste drainage.

GS transport can be observed with DCE-MRI in combination with administration of para-magnetic contrast agents (e.g. DOTAREM) into the CSF [17]. However, supplementary analysis is required to extract and distinguish characteristics of the GS transport patterns and flow. Current techniques for quantifying GS transport include assessment of brain parenchymal solute uptake or clearance [14], kinetic analysis [17], and k-means cluster analysis [13, 15]. Kinetic analysis strategies provide a static ‘snapshot’ of GS transport over 2-3 hours and provide global influx and efflux parameters however, regional information at the voxel level cannot be accurately derived due to heterogeneous transport kinetics across brain regions.

Here, we present a straightforward demonstration of the equivalence between the regularized dynamical OMT problem and the Fisher-Rao regularized variant in deriving the Lagrange coordinates associated with the flow. The Fisher-Rao perspective offers an accessible alternative methodology and deeper insight into the regularized OMT formulation as previously applied [16]. We apply the Lagrangian pipeline to measure GS function based on DCE-MRIs acquired in rats while under two different states of unconsciousness - ‘light’ sleep/hypnosis with DEXM+ISO versus ‘deep’ sleep/anesthesia with ISO based on data from [2]. Specifically, we tested the Lagrangian OMT framework’s ability to differentiate GS transport and solute drainage via pathways in the brain proper as well as along cranial efflux routes between these two different states of arousal.

Clustered pathlines derived from the Lagrangian framework are shown in Figure 2. It is evident, based on the pathlines patterns observed in the rats under the two anesthetics, that GS transport function is different. With DEXM+ISO anesthesia, GS transport is accelerated (compared to ISO) as noted by the high density of pathlines extending into the parenchyma globally and along the middle cerebral artery (MCA). In contrast, the pathlines in rats anesthetized with ISO do not align with the MCA to the same extent and parenchymal uptake is also reduced in other areas including the midbrain and hippocampus. However, paradoxically, GS transport appears vigorous in the area of the olfactory bulb including solute efflux into the nasal conchae.

Figure 3 captures pathline speed (integrated over 130 minutes) through the GS in rats anesthetized with either DEXM+ISO or ISO. The speed pathlines associated with the whole brain are shown as a 3D volume rendered color-coded map. Higher and lower magnitudes of speed are represented by red and blue colors, respectively. More speed pathlines are apparent in the frontal part of the brain in the DEXM+ISO anesthetized rats compared to ISO. Paradoxically, in the ISO anesthetized state, pathline speed appears to be higher in the nasal conchae (arrows). In Figure 4, we have captured the faster (red) and slower (green) moving speed lines, given in arbitrary units (a.u.) in an ISO anesthetized rat (Figures 4A, B) compared to a DEXM+ISO anesthetized rat (Figures 4C, D). It is obvious that the proportion of faster moving particles in the clustered pathlines is higher in ISO than DEXM+ISO anesthetized rats, as confirmed via quantitative analysis (Figure 4E). These differences suggest that solute efflux to lymphatic vessels in the nasal conchae may be acting as an efflux ‘valve’ in the setting of impaired fluid flow in ISO anesthesia.

5 Conclusion and Future Work

We presented a Lagrangian representation of Fisher-Rao regularized OMT to represent time-varying dynamics and various attributes in a single visualization. The resultant framework was then used to capture known GS transport flow patterns and provide additional insights in distinguishing GS function under different anesthetic conditions in addition to clearance of solutes via the olfactory nerves and into the nasal conchae. In the future, we will use the computed particle attributes as features in machine learning algorithms to differentiate and classify various transport related flow patterns for diagnostic purposes.

Acknowledgements

This study was supported by AFOSR grants (FA9550-17-1-0435, FA9550-20-1-0029), a grant from National Institutes of Health (R01-AG048769, RF1-AG053991), MSK Cancer Center Support Grant/Core Grant (P30 CA008748), and a grant from Breast Cancer Research Foundation (BCRF-17-193).