Unbalanced regularized optimal mass transport with applications to fluid flows in the brain

As a generalization of the optimal mass transport (OMT) approach of Benamou and Brenier’s, the regularized optimal mass transport (rOMT) formulates a transport problem from an initial mass configuration to another with the optimality defined by the total kinetic energy, but subject to an advection-diffusion constraint equation. Both rOMT and the Benamou and Brenier’s formulation require the total initial and final masses to be equal; mass is preserved during the entire transport process. However, for many applications, e.g., in dynamic image tracking, this constraint is rarely if ever satisfied. Therefore, we propose to employ an unbalanced version of rOMT to remove this constraint together with a detailed numerical solution procedure and applications to analyzing fluid flows in the brain.


Introduction
The optimal mass transport (OMT) problem is concerned with finding a transport mapping from an initial mass distribution to a final one, with optimality defined relative to a given cost function [1][2][3][4] .A reformulation of the OMT problem in a fluid dynamical framework was proposed by Benamou and Brenier 5 using the L 2 distance as the basis for the cost function, which may be written as follows in the special case of interest to us in the present work.
Given an initial mass distribution function ρ 0 (x) ⩾ 0 and a final one ρ 1 (x) ⩾ 0 both defined on a bounded region Ω ⊆ R 3 and with the same total mass (i.e., Ω ρ 0 (x) dx = Ω ρ 1 (x) dx), the dynamic OMT problem by Benamou and Brenier 5 aims to solve min where a temporal dimension t ∈ [0, T ] is added to the transport process.In the above formulation, ρ(t, x) is the dynamic density function, and v(t, x) is the dynamic velocity field defining the fluid flows from ρ 0 to ρ 1 .Equation (1b) is called the continuity equation in fluid dynamics, and characterizes the advective transport of a conserved quantity in bulk flows.The cost function (1a) is the total kinetic energy of the transport process.The square root of the achieved minimum of (1a), if it exists, is called the L 2 -Wasserstein distance between ρ 0 and ρ 1 .
As an extension of model (1a)-(1c), a regularized version of OMT (rOMT) has been developed and applied in various places in which one includes a diffusion motion into the transport; see relavent work [6][7][8] and the many references therein.More precisely, a diffusion term is added into the continuity equation (1b) to make it where the constant σ > 0 is the diffusion coefficient.Equation ( 2) is thus an advection-diffusion equation in fluid dynamics.
The rOMT model has been proven useful for a number of important tracking problems in computational fluid dynamics, such as in quantifying and visualizing the movement of solutes in the brain on dynamic contrast enhanced magnetic resonance imaging (DCE-MRI) [8][9][10][11][12] .

arXiv:2301.11228v3 [math.NA] 20 Sep 2023
As is well-known, both the OMT and rOMT models must satisfy the total mass conservation constraint, namely Ω ρ 0 (x) dx = Ω ρ 1 (x) dx.In this case, we usually call the problem as balanced to refer to the conservation of total mass.Indeed as formulated in (2), neither advection nor diffusion will change the total mass locally or globally in a given region Ω.However, for applications in dynamic imagery in which either OMT or rOMT is employed as an optical flow tracking method, this is almost never the case.For example, in DCE-MRI data where a gadolinium-based tracer is injected and delivered into the body, an early climbing period of the total image intensity is usually observed, since it takes time for the tracer to reach and fill the region of interest.Under these circumstances, if we assume that the intensity of image signal is proportional to the density and mass in the aforementioned two models, the total mass conservation law is no longer satisfied, and thus the OMT and rOMT models cannot be directly applied.Notably, analyzing the initial accumulating stage of tracers may help uncover interesting and physiologically relevant transport patterns in the brain, and therefore a new model is necessary.
In the present work, we propose an unbalanced regularized OMT (urOMT) model for applications to the DCE-MRI fluid flow data, where an independent variable and its indicator function are added as an "invisible" sink or source of mass.The urOMT problem is formulated as follows: + α χ(t, x)r(t, x) 2 ρ(t, x) dx dt (3a) where r(t, x) is the relative source variable, χ(t, x) is the given indicator function of r(t, x) which takes values either 0 or 1 to constrain r to a certain spatial and temporal location, and α > 0 is the weighting parameter of the source term in the cost function.This model takes inputs ρ 0 (x), ρ 1 (x) and χ(t, x), and solves for the optimal ρ(t, x), v(t, x) and r(t, x).The added second term in the cost function (3a) is called the Fisher-Rao term which arises from the Fisher-Rao metric in information geometry.Our model therefore can be viewed as the interpolation between the L 2 -Wasserstein and the Fisher-Rao metrics 13,14 .
The partial differential equation (3b) indicates that there are three types of physical phenomena taking place in the dynamic system, advection (∇ • (ρv)), diffusion (σ ∆ρ) and mass creation/destruction (χρr).We should note that we prefer to use the relative source r which controls the rate of mass gain (r > 0) and loss (r < 0), rather than the source s = ρr, as the unbalanced variable in our urOMT formulation, since the resulting expressions are cleaner and easier to interpret.Even though r plays a role as a sink of mass when r < 0 and a role as a source of mass when r > 0, in our work we refer to r as the relative "source" to broadly refer to its ability to generate unbalanced mass instantaneously.One can imagine that r is the source of both positive and negative mass.The main point however, is that we no longer require the total mass conservation condition for the input images ρ 0 and ρ 1 .
Our work presented here should be considered as an extension of the rOMT algorithm 10 since both share the similar numerical structure and setup.We are specifically interested in applying the urOMT method into DCE-MRI studies to quantify the fluid flows in the rat brain which is our motivation of introducing an unbalanced term to analyze unbalanced data.Although there has been a good amount of work in unbalanced OMT as mentioned above, to the best of our knowledge, this is the first work to incorporate the unbalanced regularized OMT for the quantification of dynamic fluid flows using DCE-MRI imaging.
Briefly summarizing the present paper, in the Section "Numerical Method", we give the detailed numerical method for solving the urOMT problem (3a)-(3c), and in the Section "Results", we test our urOMT model and show results of applications to studying dynamic fluid flows in both synthetic data and DCE-MRI rat brain data.In the Section "Discussion", we further discuss and summarize the urOMT method.Some relevant efforts and potential future work are also provided.Lastly in the Section "Conclusion", we conclude our present work.

Numerical Method
In this section, we elaborate on the numerical method developed for the urOMT problem (3a)-(3c) on 3D images, especially on DCE-MRI-based images in which the signal intensity levels are reflecting the concentration of the Gadolinium based tracer, and are therefore proportional to the mass/density in our model.Note that this method can be easily adapted to images in any dimension with minor changes.The numerical method used in this work is largely inspired by and based on the previous work 10 .

Model
Given a pair of 3D images, ρ img 0 (x) and ρ img 1 (x), and an indicator χ(x,t), in avoidance of the over-matching of the image noise, we consider posing a free end-point condition, so that we can remove the end-point constraint ρ(T, x) = ρ img 1 (x).Another fitting term is consequently added into the cost function.The numerical model we solve is therefore written as: where β > 0 is the weighting parameter for the fitting term in the cost function.The dynamic density function ρ(t, x) can be explicitly derived starting at ρ img 0 and following equation (4b) with a velocity field v and a relative source r with its indicator χ, so ρ is removed from the optimized variables.Basically, with this setup ρ becomes a state variable.We call this model (4a)-(4c) the unbalanced regularized OMT with free end-point (free-urOMT).

Discretization
Since 3D images are typically defined on cubical domains, we divide the cubical space Ω into a cell-centered grid size n 1 × n 2 × n 3 with uniform spacing ∆x, ∆y and ∆z in x, y and z-direction, respectively.Let n = n 1 n 2 n 3 be the total number of voxels.The time interval [0, T ] is partitioned into m equal sub-intervals with length ∆t = T m .Then we have m + 1 discrete time steps We use a bold font to denote a flattened vector discretized from its corresponding continuous function onto the cell-centered grid defined above.Therefore, the given initial and final images ρ img 0 (x) and , each denoting the mass distribution at t i , and where ρ 0 denotes the given initial image.The velocity field v(t, x), relative source r(t, x) and its indicator χ(t, x) may also be discretized into v i v i v i , r i r i r i and χ i χ i χ i for i = 0, • • • , m − 1, each denoting the velocity field, relative source, and the indicator transforming ρ i ρ i ρ i to ρ i+1 ρ i+1 ρ i+1 , respectively.We further denote . So far, we have defined the discretized variables on the space and time grids.The cost function equation (4a) may therefore be approximated by where Here ⊗ denotes the Kronecker tensor product, and ⊙ denotes the Hadamard product.Further, [•|•] represents the block matrix, and ∥ • ∥ means taking the L 2 norm of a vector.I k is the k-dimensional identity matrix for k ∈ N + .

Solving the Partial Differential Equation
Next, we deal with the partial differential equation (4b).We place ghost points outside of the boundary and employ the Neumann boundary condition such that the derivative across the boundary is always 0. The Laplacian operator ∆ may be approximated with a matrix Q under the aforementioned numerical grid and boundary condition.Similar to the technique employed in the previous work 10 , we use the operator-splitting method to numerically solve the equation but in this work we divide the equation into three steps.To be precise, at each time step from t i to t i+1 for i = 0, • • • , m − 1, we divide the whole process (4b) into first, mass gain/loss: ∂ ρ ∂t = χρr; second, advection: ∂ ρ ∂t + ∇ • (ρv) = 0; and third, diffusion: ∂ ρ ∂t = σ ∆ρ, in total three steps, and then integrate them together.
For the first mass gain/loss step, given an initial condition ρ(t i , x) = ρ i (x), the equation may be discretized into Numerical Pipeline of urOMT.From t i to t i+1 for i = 0, • • • , m − 1, the interpolated image ρ i ρ i ρ i is firstly added with mass by applying matrix R(r i r i r i ), and is secondly advected via the velocity field v i v i v i by applying the averaging matrix S(v i v i v i ) and is lastly diffused by applying matrix L −1 .
where 1 n 1 n 1 n is a vector of length n consisting of 1's.The second advection step with an initial condition ρ(t i , x) = ρ src i (x) can be discretized and solved with where S(v i v i v i ) is the averaging matrix with respect to v i v i v i using the particle-in-cell method which redistributes the transported mass to its nearest neighbors by a certain ratio.The ( j, k) entry of S(v i v i v i ) is the ratio of mass allocated from the old location k to the new location j.Multiplying S(v i v i v i ) by a vector which represents a mass distribution, we can derive a new mass distribution transported by the velocity v i v i v i under the pre-defined numerical grid.The third diffusion step with an initial condition ρ(t i , x) = ρ adv i (x) employs the Euler Backwards scheme in the following manner: Combining all three steps (7b), ( 8) and (9b), we have If we denote where diag(•) is the operator turning a vector into a diagonal matrix.Equation ( 10) may be re-written as and R(r i r i r i ) are n-dimensional matrix with respect to v i v i v i and r i r i r i , respectively.See Fig 1 for the pipeline of the numerical dynamics.

Computing the Gradient and the Hessian
Notice that in the discrete problem (13a)-(13c), we can equivalently re-write the constraint (13b) as for i = 1, • • • , m, which explicitly gives the expression of mass distributions at all time steps from a given initial mass distribution ρ 0 ρ 0 ρ 0 , a velocity field v v v and a relative source r r r.One can prove that Then according to equations ( 14) and (13c), r r r) can both be explicitly written in linear to v v v and r r r.Therefore, the optimization problem (13a)-(13c) can indeed be viewed as an unconstrained minimization problem.Observe that in the cost function, Γ 1 is linear to r r r and its component v v v ⊙v v v is quadratic to v v v; Γ 2 is linear to v v v and its component r r r ⊙r r r is quadratic to r r r; Γ 3 is quadratic to both v v v and r r r.It is then natural to employ the Gauss-Newton method to optimize on the problem, which involves computing the gradient and the Hessian matrix of Γ with respect to variables v v v and r r r.
Next, we focus on calculating the gradient of Γ: and the Hessian matrix of Γ: , where Equation ( 14) shows that and is thus independent of v j v j v j and r j r j r j for j ⩾ k.Defining v j = 0 and J k r j r j r j = 0 always hold for j ⩾ k.If we further denote then J v v v and J r r r are lower-triangular block matrices of the form where denotes the row block of J v v v , and J k r r r for that of J r r r for k = 1, • • • , m.With the notations defined above, then for the gradients we have and where we use for simplicity of notation.For the Hessian matrix H, we use a function handle which is a MATLAB data type in anticipation of solving the linear system Hx x x = −g in the next step.To be specific, we compute vector Hx x x where x r ∈ R mn rather than matrix H.To satisfy the symmetry of the Hessian matrix and to omit some complex second-order terms, we have the approximations The motivation of using a function handle Hx x x instead of the matrix H itself is to avoid numerical multiplication of two big matrices which could be time-consuming.With the function handle, for example, in (22c) we can instead multiply a matrix by a vector twice to arrive at the desired computation.
As for the formulation of J k v j v j v j and J k r j r j r j , by the structure of (13b) we have that for where is dependent only on ρ j ρ j ρ j ,v j v j v j because S(v j v j v j ) linear to v j v j v j , and With the explicit and recursive expressions in equations ( 23) and ( 25), one can compute and their transpose multiplied with a vector in an iterative manner.J T v v v and J T r r r multiplied with a vector can also be computed recursively due to their lower-triangularity.

Algorithm
With the analytic formulation of the gradient g and the Hessian handle Hx x x given above, we can then utilize the Gauss-Newton method to find the optimal solution.See Algorithm 1 for the pseudo-code.
If we have more than two successive given images, ρ img χ img q−2 between adjacent images where q > 2 and q ∈ N + , we can run the algorithm iteratively between each pair of adjacent images to derive prolonged velocity fields and relative sources.In other words, the process can be graphed as where the urOMT algorithm is run for q − 1 times and therefore q − 1 successive outputs are returned.If one would like the prolonged velocity fields to be smoother in the temporal dimension, one can put the last interpolated image ρ m ρ m ρ m of the previous loop into the next loop as the initial image to avoid constantly introducing new data noise into the system 10 .Compute S(v j v j v j ), R(r j r j r j ) and B(ρ j ρ j ρ j ,r j r j r j ) for j = 0, • • • , m − 1;

Algorithm 1 Gauss-Newton Method
6: Compute gradient g and the Hessian function handle Hx x x according to equations ( 20)-( 22); 7: Solve linear system Hx x x = −g for x x x; Update [v v v;r r r] = [v v v;r r r] + lx x x; 13: end for 14: return [v v v;r r r];

Results
In this section, we give some examples illustrating the application of the urOMT methodology.Obviously, the model admits both the Eulerian and Lagrangian perspectives for post-processing.The Eulerian formulation focuses on the fluid flows at fixed locations over time, while the Lagrangian formulation enables one to follow the trajectory of a given particle over time, and therefore to analyze the features along the given trajectory.

6/16
Specifically, suppose we run the urOMT algorithm on given density images ρ img χ img q−2 with discretization described before, the algorithm will run q − 1 successive loops and return the solutions for k = 1, • • • , q − 1 where the subscript k stands for the k-th loop.By taking the Therefore, r * s * k,i are called the Eulerian relative source maps and Eulerian speed maps, respectively, for k = 1, • • • , q − 1 and j = 0, • • • , m − 1, and they allow us to observe the fluid flows at fixed coordinates at different discrete time steps.For ease of visualization, we define the time-averaged Eulerian speed map and time-averaged Eulerian relative source map between respectively, where 0 In contrast, the post-processing framework developed and applied in the previous work [8][9][10]12 that follows Lagrangian coordinates present data as binary trajectories of the fluid flows, which we refer to as the pathlines in our work. Byconnecting the starting and terminal points of the pathlines, we have the velocity flux vectors, which are also called the displacement fields in physics.These vectors may be used to visualize the direction and the distance travelled in a compact and interpretable manner.If we endow the pathlines with more information, for example, speed (the L 2 norm of the velocity field) and Péclet (Pe) number (i.e., the ratio of the rate of advection to diffusion), we can derive what we call speed-lines and Péclet-lines, respectively.More details of this Lagrangian method may be found in Koundal et al. 8 and Chen et al. 10 Now that we have briefly described the two post-processing methods, the Eulerian and the Lagrangian perspectives, we next utilize two datasets to exhibit the results of the urOMT analysis.The first dataset is a synthetic geometric dataset derived from Gaussian spheres as a simple demonstration of the urOMT method.The second is the DCE-MRI data from a rat brain where we elucidate the application of urOMT in in vivo datasets for quantifying and visualizing the fluid flows.

Tests on Gaussian Spheres
We first start with creating five successive 3D images from Gaussian spheres denoted as in order to test the urOMT algorithm.These spheres were created to first gain mass and later lose mass in the center region of the spheres over time, and in addition the spheres are spatially transported forward and are also exhibiting active diffusion over time (Figure 2a).Specifically, ρ i was created from a 3D Gaussian function then advection is naturally included in the whole process to reflect the forward translation of the center of the Gaussian spheres from top-left to bottom-right (Figure 2a).We define the center region of these Gaussian spheres as the region within a radius of 1.5 of the center [0.8i, 0.8i, 0.8i] which are denoted as χ i (x, y, z) for i = 0, • • • , 4. Then we apply χ i to G i to imitate mass gain and loss in the center region which are given by (1 + a i χ i (x, y, z))G i (x, y, z) where [a 0 , a 1 , a 2 , a 3 , a 4 ] = [0, 0.1, 0.2, 0.1, 0].We then discretize these functions by taking a uniform spatial length and fitting them into the same numerical grid of size 50 × 50 × 50.Diffusion was further added to ρ i by applying a MATLAB inbuilt 3D Gaussian filter imgaussfilt3 to ρ i with standard deviation = (i + 1) √ 0.2 for i = 1, • • • , 4. One can imagine the five successive images as visualizing a "wormhole", which is moving forward (advection) with mass diffusing into its surrounding area.At the same time, at the center region of the wormhole mass is gained (ρ img 0 to ρ img 2 ) or lost (ρ img 2 to ρ img 4 ) from or to another interconnected space.

8/16
In our first test on the Gaussian sphere data, we fed all five images, ρ img 0 , ρ img 1 , • • • , ρ img 4 , into the urOMT algorithm.From ρ img i to ρ img i+1 , we utilize the indicators as the linear translation of χ i to χ i+1 (the region within a radius of 1.5 centered at [0.8(i + j m ), 0.8(i + j m ), 0.8(i + j m )] for j = 0, . . ., m, i = 0, . . ., 3) to only allow mass gain and loss to occur in the center regions.The parameters used in this experiment are listed in Table 1.Computations were run with MATLAB 2018b on the departmental High Performance Computing cluster at Memorial Sloan Kettering Cancer Center with Red Hat Enterprise Linux 7.5 operating system using 3 CPUs and 128GB of memory which took about 3 hours and 20 minutes.
From the Eulerian perspective, we show the time-averaged Eulerian speed maps and relative source maps between every pair of input images (Figure 2b,c) and between ρ img 0 and ρ img 4 (Figure 2d,e).The algorithm recognized the core regions of the spheres as having higher speed.Moreover, it successfully captured the mass gain and loss patterns during the entire process with the red color indicating the initial mass gain from ρ img 0 to ρ img 2 and the blue color indicating the later mass loss from ρ img 2 to ρ img 4 .The Eulerian relative source map is restricted in the center region of the spheres because of the indicator functions.With Lagrangian post-processing, we derived the pathlines, velocity flux vectors, speed-lines and Péclet-lines under Lagrangian coordinates (Figure 2f-i).The pathlines exhibit what trajectories of particles would look like over time if they were placed at given initial points in the system at t = 0.The funnel-shape of the pathlines is a reflection of the accumulated effect of diffusion which gradually disperses the mass.The velocity flux vectors are also provided to show the direction and distance of the whole transport process.The speed-lines, i.e., pathlines endowed with speed, indicate that higher speed occurred mainly in the core of the spheres which is in agreement with the Eulerian speed map.The Péclet-lines show that initially the transport was dominated by advection and later by diffusion.
To further understand how indicators affect the results and how mass is transferred in the urOMT system, we performed a second test with only two input images ρ img 0 , ρ img 1 and indicators as a constant 1 to allow mass gain/loss everywhere.The parameters used in this test are listed in Table 1.The computational runtime was about 30 minutes with the same machine and configuration as the first test.Eulerian relative source map (Figure 3a) and Eulerian speed map (Figure 3b) were returned from the urOMT algorithm.We observe that the top-left of the Eulerian relative source map was negative and colored in blue where the input sphere leaves and the bottom-right was positive colored in red where the sphere arrives.In other words, without any constraint on the relative source term, the Eulerian relative source map can provide global information on the transport status in terms of mass loss and arrival.2a, Eulerian outputs, shown in 3D rendering, were returned from the urOMT analysis.c: The transport in the system can be separated into two channels where the advection and diffusion take place in the R 3 space and the relative source pushes or draws mass between the two channels.
In this case, mass is freely transferred via two "channels", the real R 3 space and the imaginary source layer (Figure 3c).Advection and diffusion occur only in the R 3 space.The source layer can generate infinite amount of mass and push it to the corresponding location in the R 3 space; it can also draw a certain amount of the mass from the R 3 space to the corresponding 9/16 location in the source layer.The activities in both channels happen simultaneously.Suppose in R 3 , one would like to move mass ρ 0 > 0 at location x 0 to mass ρ 1 > 0 at location x 1 , i.e., . We denote x ′ 0 and x ′ 1 as the corresponding locations of x 0 and x 1 in the source payer, respectively; the amount of mass transported in R 3 from x 0 to x 1 as ρ trans 0→1 > 0; the amount of mass drawn from x 0 to x ′ 0 as ρ src 0 > 0; and the amount of mass pushed from x ′ 1 to x 1 as ρ src 1 > 0. Therefore, the following equations hold The weighting parameter α in the urOMT formulation therefore balances the split of ρ 0 into ρ trans 0→1 and ρ src 0 whose effect will be demonstrated in the next section.
In this test, to move ρ img 0 to ρ img 1 , part of the mass was transported in R 3 and part in the source layer, and these two phenomena facilitate each other.We know a priori from the creation of the data that the sphere should be transported forward and mass gain only occurs in the center region.However, allowing a global relative source, the urOMT algorithm finds an easier way to transform ρ img 0 into ρ img 1 by making use of the source layer to transfer mass.Most importantly, removing the indicator constraint the relative source globally compensates the transport in R 3 and is indicative of the leaving and arrival of mass.This can be very useful in the context of applications to some real dataset.As a matter of fact, a priori indicator for the source layer is sometimes unavailable due to the complexity of the system.

Application to Rat Brain MRI
A very useful application of the urOMT model is to quantify the transport properties of brain fluids with DCE-MRI data.It still remains debated in the scientific community how fluid and solutes are transported in brain parenchyma and how the "dirty" fluid is drained out of the brain to maintain homeostasis 25,26 .The relative source in the urOMT model may be helpful in revealing fluid and solute clearance patterns in the brain.In this experiment, our urOMT method was applied to 3D DCE-MRI dataset derived from a 3-month-old healthy rat brain.The tracer, gadoteric acid, was injected into the cerebrospinal fluid (CSF) of the rat after the rat was anesthetized.The DCE-MRI data series of the rat brain was collected every 5 minutes and lasted for 140 minutes, ending up with 29 images in total.The MRI signal images were then processed to derive the %-signal change from the baseline to approximate the concentration images of tracers.More information of the DCE-MRI data may be found in Chen et al. 9 In this numerical experiment, we assume that the intensity of the DCE-MRI images is proportional to the mass in the urOMT model.We used every other image within the brain region as input images in order to save running time, which resulted in a total of 15 images: ρ img 0 , ρ img 1 , • • • , ρ img 14 for the urOMT algorithm (Figure 4a).In the previous work 9 , the rOMT model was applied right before the peak of the total signal intensity of the input images, i.e., ρ img 3 , • • • , ρ img 14 in current notation.In the present experiment, due to the introduction of the relative source term in the model we were able to include earlier frames when the total intensity was still rapidly increasing (Figure 5a, the red dashed curve).Since the mechanism of the fluid transport in brains is complex and still under intense investigation, there is no information provided a priori for the relative source term.So we set its indicator function χ to be equal to 1 everywhere with the assumption that the entire brain system is "leaky" to allow the tracer to enter and exit the system through unknown ways.In order to derive smooth prolonged dynamics, we used the last interpolated image from the previous numerical loop as the initial image in the next loop.The parameters used in this experiment are listed in Table 1.The computation was also run with MATLAB 2018b on our departmental High Performance Computing cluster at Memorial Sloan Kettering Cancer Center with Red Hat Enterprise Linux 7.5 operating system using 3 CPUs and 128GB of memory which took about 8 hours and 50 minutes.
We show the Eulerian relative source maps and speed maps between every other input frame (Figure 4b-c).The Eulerian relative source maps indicate that the tracer first entered the CSF surrounding the brain causing the intensive mass gain noted in the first frames, and then subsequently moved into deeper brain tissue regions resulting in mass loss in CSF and mass gain in the brain tissue.The Eulerian speed maps indicate that the initial speed of tracer when entering the CSF was very high, and slowed when the tracer penetrated deeper, probably due to diffusion dominating the transport in the tissue.From the Lagrangian perspective, we show the pathlines of tracers starting at t = 0 as well as those lines endowed with speed and Péclet number (Figure 4d,f-g).The pathlines signifies the pathways/trajectories of the tracer entering CSF and brain tissue.The speed-lines show that higher speed was mainly in the CSF and along the perivascular space of the large vessels, and further that the transport was identified as advection-dominated according to the high values in the Péclet-lines in those regions.The velocity flux vectors were also derived and demonstrated that the tracer entered the brain tissue in a symmetrical pattern about the the midline of the skull base (Figure 4e).
Recall that in the urOMT model (3a)-(3c), the weighting parameter α > 0 penalizes the source term in the cost function.In theory, as α → +∞, r gets suppressed and this model approximates the rOMT model where unbalanced mass gain and loss are not allowed.In the test, we used above α = 10000.We ran further tests with different values including α = 1000, 3000, 6000, 20000 and 50000 to demonstrate the effect of α on the results.
From the urOMT algorithm, we can compute the final interpolation ρ i,m for i = 1, • • • , q − 1.To test the numercial performance, we define the normalized mean squared error (NMSE) between each pair of final interpolation ρ i,m and the ground truth ρ img i as We also define the percent change in total mass (PCTM) between ρ i,m and ρ img i as where sum(•) denotes taking the sum of all entries in a vector.Both NMSE and PCTM measure the accuracy of the urOMT algorithm and the lower they are, the more accurate the urOMT results are.However, NMSE measures the closeness of the ground truth and the interpolations in a local manner, while PCTM in a global manner.We emphasize that urOMT is a data-driven method; in other words, the transport between two images is a black box and urOMT is simply giving the most likely evolving path under a pre-defined least cost assumption.Therefore, it is difficult to define the ground truth between two input images for us to compare with.By fixing the rest of the parameters but using various α values (1000, 3000, 6000, 10000, 20000 and 50000), we ran our novel urOMT algorithm and the post-processing procedure on the same rat brain data.In Figure 5a, we plot the curves of the total image intensity of both input images and interpolations with different α's.From the curves, the smaller the α is, the closer the curve of interpolations is to the curve of input images.This is in agreement with the theoretical observation that the smaller the α is, the less the source term is penalized, which means more mass gain/loss is allowed in the system in order to adjust the total mass.From Figure 5b-c, we found that a smaller α may help derive more accurate numerical simulations, because for both NMSE and PCTM the trend is that the larger the α is, the higher the NMSE and PCTM values are.This makes sense because when the system is highly unbalanced, and forcefully using a high α produces a nearly balanced environment which contradicts with the data setting.To further examine the effect of α on the quantitative and visual results, we plot the Eulerian relative source maps and Eulerian speed maps from ρ img 6 to ρ img 8 and speed-lines with different α's in Figure 5d.They show that when the source term is greatly penalized with a high α, the speed is consequently elevated.For example in Figure 5d with α = 50000 where the model most approximates rOMT, there is very high speed transport color coded in red at the base of the  brain, which could be artificial and potentially over-estimated because in such a system mass has to move more quickly in order to match the final input image.In contrast, in a system with low α where instantaneous mass gain/loss is promoted, mass does not need to transport in the same degree to match the final input image, since it can instead pull in (or push out) mass from (or into) the "invisible sink" (the source layer in Figure 3c) via the relative source r.For example in Figure 5d with α = 1000, the rapid movement of tracer at the base of CSF almost disappeared because the change of the system is mostly accounted by the relative source.In general, the relative source r and the velocity field v compensate each other and their effects in the system are actively controlled by the weighting parameter α. system by adding an independent relative source term into the formulation, while the rOMT model requires the total mass to be conserved.As discussed before, both theoretically and numerically the rOMT model can be approximated by the urOMT model when the parameter α goes to infinity.In other words, urOMT "incorporates" rOMT, and one can use urOMT with a large enough α (thereby assuring that the mass conservation condition of input images to be met) in place of the previous version of rOMT.As such the new urOMT introduced in the present work is a more powerful and flexible model for analysis of fluid transport.
The urOMT method may be particularly useful for studying the cross-talk between the glymphatic system and meningeal lymphatics.See relevant work 25,27,28 for more details about how the two systems interplay.With the additional information of the relative source in our urOMT model which reveals mass gain and loss, we are now able to observe and quantify the solute and fluid entering and exiting the two systems.For example, in Figure 4b, at first we see mainly mass gain (colored in red) indicating that the tracer is flowing into the CSF, and then slowly we see blue color along the skull base, indicating that the tracer is either redistributing into the tissue bed or exiting via the draining lymphatic vessels.Therefore, our urOMT method has the potential to probe the clearance pattern in more depth at the level of the CSF and tissue compartments.
In the test on the rat brain DCE-MRI, we posed no spatial constraint on the relative source r and used an indicator χ all 1's given that so far there is no agreed upon answer yet on the underlying transport mechanisms and the exact drainage pathways from the brain 25,26 .Indeed from the DCE-MRI data (Figure 4a), we did not observe a specific anatomical efflux route of the tracer out of the brain but the total intensity curve did decline in later frames (Figure 5a, the red dashed curve).In this case, we assume that the system is leaky and the tracer can be transferred via unknown tunnels in brain (a correspondence of the source layer in Figure 3c) to form the given images.Similarly in the widely used Toft's model 29,30 for pharmacokinetic analysis of DCE-MRI studies of tumors where the signal curves are also unbalanced, a leakage between the blood vessels and the tissue is assumed to occur everywhere in the images.The popular parameter K trans returned from the Toft's model is used to quantify the local leakage of gadolinium-based tracers from the blood to the tissue 29,30 .
In the urOMT algorithm, the difference between two input images is mainly captured by either the mass gain/loss rate r in the source layer or the velocity field v in R 3 , and α is the parameter balancing the two.Indeed as demonstrated by other work 17 , with a decreasing α parameter the transport (characterized by the velocity field v) is being compensated by mass gain and loss (characterized by the relative source r).With the rat brain MRI dataset, we demonstrated that a low α value gives higher numerical accuracy, but produces decreased speed by allowing the relative source term to play a greater role in the dynamics.Thus, one needs to be aware of the trade-off between the accuracy and the strength of fluid flows by choosing an appropriate parameter α when applying urOMT.One approach would be to make use of the indicator function χ to restrict the behaviors of the relative source r within a certain region if the entering and exiting information of the fluid is known beforehand.Another approach is to use a time-varying α given that the total intensity curve of the input image is usually not linearly increasing or decreasing.Indeed, we plan to explore this possibility in some future work.
Other than the parameter α, there are also additional parameters worthy of examining, such as the weighting parameter β for the fitting term in the cost function (4a) and the constant diffusion coefficient σ .Some preliminary efforts have been made in using a non-linear diffusion term in the partial differential equation (2), given the non-constant diffusion phenomena in brain fluid flows 31 .Future work includes adding a non-linear and spatially dependent diffusion term in the current urOMT formulation.
The urOMT model was largely motivated by the changing pattern of the total signal curve from DCE-MRI experiments.Specifically in rat brains, the temporal signal intensity curve typically peaks at approximately 1 hour after tracer injection, and later on, either keeps decreasing (Figure 5a, red dashed curve) or reaches a plateau (this depends on the total amount of tracer administered), signifying that the total mass may be highly unbalanced over time in the system.Given that the DCE-MRI protocol has been widely used in cancer imaging in clinics [32][33][34] , the urOMT method also has the potential to be applied to tumor DCE-MRI data in human experiments to investigate the tumor vasculature, and to help pave the way for new medical treatments.

Conclusions
The urOMT methodology incorporates both advection and diffusion motions into the transport process, as well as allowing for mass gain and loss in the dynamic images by introducing a relative source variable.For special cases, it may also constrain the relative source to a given region or time interval.As an extension of the rOMT model [8][9][10] , the urOMT model removes the total mass conservation constraint, while keeping the attractive advection-diffusion framework, making it applicable to modeling the fluid flows in the brain under DCE-MRI protocol and many other real-world modeling problems in computational fluid dynamics.

Figure 2 .to ρ img 2 and was lost from ρ img 2 to ρ img 4 in
Figure 2. The First Test on 3D Gaussian Spheres.a: Five successive images, shown in 3D rendering, were created from Gaussian spheres as inputs into the urOMT algorithm.In addition to advection (from top-left to bottom-right) and diffusion included in the transport process, mass was gained from ρ img 0 to ρ img 2 and was lost from ρ img 2 to ρ img 4 in the center region.b-e: Under Eulerian coordinates, the time-averaged speed maps and relative source maps visualized in 3D indicate the speed and mass gain/loss distribution in the domain, respectively.b and c are derived between every pair of input images; d and e are derived between ρ img 0 and ρ img 4 .f-i: Under Lagrangian coordinates, the binary trajectories of the transport are recorded by pathlines color-coded with start and end points.Connecting the start and end points of pathlines, we derive the velocity flux vectors illustrating the direction and distance of the overall movement.By endowing the pathlines with speed and Péclet (Pe) number, we derive the speed-lines and Péclet-lines, respectively.

Figure 3 .
Figure 3.The Second Test on 3D Gaussian Spheres.a-b: Given input images ρ img 0 and ρ img 1 visualized in Figure2a, Eulerian outputs, shown in 3D rendering, were returned from the urOMT analysis.c: The transport in the system can be separated into two channels where the advection and diffusion take place in the R 3 space and the relative source pushes or draws mass between the two channels.

Figure 4 .
Figure 4. Application to 3D Rat Brain MRI.a: Rat brain MRIs, shown in 3D rendering from the right-lateral view plane, were fed successively into the urOMT algorithm.b-c: As returned outputs, the Eulerian time-averaged relative source maps and speed maps between every other image were plotted, indicating the rate of mass gain/loss and speed distribution over time, respectively.d-g: Under Lagrangian coordinates, pathlines, color-coded with start and end points, show the trajectories of the tracers in brain.The speed-lines show the speed values along pathlines, and similarly Péclet-lines indicate whether the transport is advection or diffusion-dominated along pathlines.The velocity flux vectors show the direction and distance of the transport.All lines and vectors are shown from both the right-lateral and bottom views, and are overlaid on the anatomical data in gray.

Figure 5 .
Figure 5. Examination of the Effect of the Parameter α. a: The comparison of the total image intensity curve of the input images and interpolations with different α's.b: The normalized mean squared error curve over the numerical steps with different α's.c: The percent change in total mass curve over the numerical steps with different α's.d: The comparison of the Eulerian relative source maps (first row) and Eulerian speed maps (second row) from ρ img

Table 1 .
Parameters used in the urOMT algorithm.