A multi-scale sub-voxel perfusion model to estimate diffusive capillary wall conductivity in multiple sclerosis lesions from perfusion MRI data

We propose a new mathematical model to infer capillary leakage coefficients from dynamic susceptibility contrast MRI data. To this end, we derive an embedded mixed-dimension flow and transport model for brain tissue perfusion on a sub-voxel scale. This model is used to obtain the contrast agent concentration distribution in a single MRI voxel during a perfusion MRI sequence. We further present a magnetic resonance signal model for the considered sequence including a model for local susceptibility effects. This allows modeling MR signal–time curves that can be compared to clinical MRI data. The proposed model can be used as a forward model in the inverse modeling problem of inferring model parameters such as the diffusive capillary wall conductivity. Acute multiple sclerosis lesions are associated with a breach in the integrity of the blood brain barrier. Applying the model to perfusion MR data of a patient with acute multiple sclerosis lesions, we conclude that diffusive capillary wall conductivity is a good indicator for characterizing activity of lesions, even if other patient-specific model parameters are not well-known. Author summary The use of advanced brain imaging techniques has supported in-vivo research targeted to the integrity of the blood-brain barrier. We propose a new type of post-processing for raw image data using contrast agent perfusion simulations on the data-poor capillary scale. Combining modern simulation techniques with the clinical image data allows us to determine patient-specific and pathologically relevant parameters such as the capillary wall conductivity. The presented simulation model is a step towards the quantification of contrast agent leakage in the brain, which is typical for acute multiple sclerosis lesions, but also occurs with other diseases affecting the blood-brain-barrier, such as cerebral gliomas.

Multiple sclerosis (MS ) is characterized by a cascade of inflammatory reactions that 2 result in the formation of acute demyelinating lesions (MS plaques). Acute lesions are 3 associated with an impaired blood-brain-barrier (BBB ) [1]. In healthy brain tissue, the 4 tight junctions between endothelial cells forming the blood vessel walls, are an efficient 5 barrier for most molecules in the brain capillaries. In active MS lesions tight junctions 6 have been found to be damaged or open [2]. Due to an auto-immune reaction, 7 immunological cells can pass the BBB and attack the myelin sheath covering the 8 electrical pulse conducting axons, leading to dysfunctions of the central nervous 9 system [3]. Magnetic resonance (MR) enhancement, using contrast agents such as 10 when assessing the integrity of the blood-brain barrier (BBB) [10,11]. In a typical 23 DSC-MRI study, contrast agent is administered intravenously (bolus injection) and 24 whole brain MR image sequences are recorded with a repetition time of about two 25 seconds over a few minutes [11]. Normal appearing white matter is distinguished from 26 inflammatory plaques by image contrast and differences in intensity-time curves. Using 27 adequate post-processing techniques, qualitative assessment of leakage coefficients 28 allows to identify contrast-enhancing lesions in an automated way [12]. Although today, 29 perfusion MRI is not considered a standard procedure in the neuro-imaging workup of 30 MS, it enables a classification of lesions according to parenchymal leakage of an MR 31 contrast agent due to differences in perfusion behavior [13]. Perfusion imaging, both 32 DSC and dynamic contrast enhanced (DCE ), may provide information about the 33 leakiness of the tissue under investigation. In this work, we investigate DSC-MRI. 34 However, the extension of the method to DCE-MRI is conceivable. 35 For the interpretation of images obtained in a DSC-MRI study, the gray scale image 36 the underlying reasons for a particular intensity-time curve of a voxel, by identifying 41 and analyzing the model parameters which are able to reproduce the MRI data. This 42 process is also known as solving the inverse problem. To this end, the model parameters 43 are tuned by using parameter estimation techniques. Forward models are typically 44 based on a two-compartment pharmacokinetic tracer model and are parameterized by a 45 small number of parameters [14][15][16]. Figure 2 visualizes a two-compartment model 46 conceptually, with compartments representing plasma and extra-vascular, extra-cellular 47 space. The plasma compartment is supplied by a flux determined by an arterial input 48 function (AIF ) [17]. The AIF can be estimated from voxels that are mostly constrained 49 within a larger afferent artery [18]. The plasma compartment exchanges mass with the 50 extra-vascular, extra-cellular space proportionally to its permeability-surface product.
magnetic resonance (NMR) exploited to acquire the MR image. There have been many 58 suggestions for improving the modeling of the latter process [20][21][22][23]. The authors 59 of [24,25] show that the local, sub-voxel tissue structure has a significant effect on the 60 NMR signal. However, all previous studies, including the recent study by [23], rely on 61 state-of-the-art two-compartment models for the perfusion process providing only 62 average concentrations in two tissue compartments within a voxel. 63 To overcome the limitations of two-compartment models, we present a perfusion 64 model on a sub-voxel scale, including the capillary network structure. Fully, 65 three-dimensionally resolved fluid-mechanical models of brain tissue perfusion imply 66 prohibitively complex and computationally expensive simulations due to the large 67 number of vessels, their non-trivial geometrical embedding, and the complex geometry 68 of the extra-vascular, extra-cellular space [26]. To reduce complexity, we use a represented by a homogenized three-dimensional continuum. The model reduction, 72 which is described in more detail in the following, leads to a coupled system of 73 one-dimensional partial differential equations for flow and transport in the vessels, and 74 three-dimensional partial differential equations for flow and transport in the 75 extra-vascular space. Related models have been used to study the proliferation of cancer 76 drugs [27][28][29], the transport of oxygen [30][31][32][33][34], and nano-particle transport for 77 hypothermia therapy [35]. A recent study [36] describes contrast agent perfusion based 78 on diffusive transport with a mixed-dimension model. The herein presented 79 fluid-mechanical model is similar to the drug proliferation model described in [27] and 80 introduced by [28]. It is derived here for the specific application of contrast agent 81 perfusion in brain tissue.

82
The fluid-mechanical model is coupled to an NMR signal model. simplified [37]. Furthermore, we employ a homogenized continuum model for blood, 106 using an apparent viscosity,μ B . Assuming a constant hematocrit of 45%, the apparent 107 viscosity can be described as a function of the effective vessel radius, R, by an empirical 108 relation [38], derived from experimental data, The effective vessel radius, R, is chosen with respect to the cross-section area,

110
A v = πR 2 . Blood density is assumed constant, ρ B = 1050 kg m −3 [39]. Under these 111 assumptions, the flow in the lumen of a capillary vessel can be described by fraction, x v , [37] 120 where M c is the molar mass of the contrast agent, ρ m,B the molar density of blood, and 121 D c B the binary diffusion coefficient of the contrast agent in blood. The exchange with 122 the extra-vascular compartment is modeled by the fluxq c , with units of kg s −1 m −1 .

123
The shape factor ω > 0 reflects the variation of axial velocity profiles in vessel 124 cross-sections [37], where χ(r), φ(r) are the dimensionless velocity profile and the dimensionless 126 concentration profile, respectively. As it has been observed that small nano particles are 127 likely to be distributed evenly [40], we choose ω = 1.

128
In the following, we consider the Gadolinium-based contrast agent Gadobutrol.
where µ P = 1.32 Pa s [43] denotes the blood plasma viscosity, T the temperature in K, 139 and k B the Boltzmann constant. phase, the interstitial fluid, through a porous medium is described by Darcy's law [44] 144 where ρ I , µ I are density and viscosity of the interstitial fluid, v t the filter velocity where φ denotes the porosity, the ratio of pore volume to total volume in a pathways. Here, we consider only transport by advection and diffusion, following [39]. across Γ, which is inversely proportional to wall permeability and wall thickness.

175
According to Starling's hypothesis [49,50], the transmural flux of a fluid is proportional 176 to the hydrostatic and colloid osmotic pressure gradient between capillary lumen and 177 interstitial space where L p is the filtration coefficient, with units of m Pa −1 s −1 , S = 2πR is the 179 circumference of the vessel, is the average hydrostatic pressure on the vessel wall, π v , π t , denote the osmotic 181 pressure in capillary lumen and interstitial space, respectively, and 0 ≤ σ ≤ 1 is the 182 osmotic reflection coefficient. The difference in osmotic pressure results from large 183 plasma proteins in the blood stream, such as albumin, and effectively draws fluid into 184 the vessels. For the in silico experiments, we assume the osmotic pressures to be 185 constant, with ∆π = π v − π t = 2633 Pa [51]. Furthermore, we choose σ = 1, 186 corresponding to the vessel wall modeled as a perfect selectively permeable membrane. 187 The contrast agent is assumed to be transported by advection with the plasma, as 188 well as by molecular diffusion. The reduction of the vessel wall to a surface leads to a 189 concentration jump across the vessel wall, which is inversely proportional to diffusive 190 wall conductivity and wall thickness. The transmural transport is described as [50] 191 where D ω is the diffusive wall conductivity, with units of m s −1 , is the average contrast agent mole fraction on the vessel wall, denotes the mole fraction in upwind direction, and 0 ≤ σ c ≤ 1 denotes the solvent-drag 194 reflection coefficient. As the considered contrast agent is a small molecule and the 195 endothelial tight junctions are damaged in lesion tissue, we set σ c = 0, neglecting 196 reflection. Determining D ω from MRI data is the major objective of this work.

197
The mass balance Eqs.
so that q m is a line source restricted by the Dirac delta function δ Λ to the center line of 206 December 20, 2018 10/45 a vessel. Analogously, we set q c = −q c δ Λ , for the source term in Eq. (7).

219
The injected fluid thus forms a sharp bolus. However, the bolus disperses significantly 220 before it reaches the brain capillaries. Therefore, the concentration inflow profile to the 221 capillary network has to be estimated from the parameters of the bolus injection. To this end, we use an ansatz from [22] 223 In summary, the complete coupled fluid mechanical model of tissue perfusion reads as: 238 times, T 1 , T * 2 , relaxes the magnetization into the initial state aligned with B 0 .

258
According to [22], the GRE-EPI voxel signal during a DSC-MRI perfusion sequence can 259 be modeled as where the repetition time T R , is the time between two RF pulses, and the effective echo 261 time, T E , is the time between RF pulse and signal readout. The base signal S 0 > 0 262 depends, i.a., on tissue proton density and the MR scanner hardware. In the following, 263 we look only at the normalized signal S n (t) = S(t)S −1 pre , where S pre is the signal before 264 the contrast agent bolus arrives in the tissue sample. The pre-contrast signal, S pre , contains all constant factors in Eq. (17), including S 0 . It follows from Eq. (17) that a 266 shortening of T * 2 results in a decrease of NMR signal strength, while a shortening of T 1 267 results in signal enhancement.

268
The following two sections introduce the models for the relaxation rates that depend on the spatial and temporal evolution of the contrast agent sub-voxel distribution of the contrast agent concentration. We follow [22], to develop a 276 model considering the spatial and temporal distribution of the contrast agent.

277
Transversal relaxation in tissues with locally heterogeneous microstructure 278 The transversal relaxation rate, R * 2 , depends on the complex local microstructure of the 279 tissue [24] and is altered by the presence of the contrast agent. We are only interested 280 in the signal change relative to the baseline, so we split the relaxation rate in a static 281 pre-contrast contribution and a time-dependent contribution depending on the contrast 282 agent concentration, The relaxation rate for a sub-voxel control volume can be described by contributions 284 of three compartments, the vascular compartment (B), the extra-cellular, extra-vascular 285 space (I), and the cellular compartment (S), weighted by their volume fractions, φ B , φ I , 286 According to [25], the rate in each compartment κ ∈ {B, I, S}, comprises 288 contributions on three spatial scales The rate R * 2,macro describes effects of static local inhomogeneities of the magnetic field 290 B 0 , which are time-independent. Since the static effects do not depend on the contrast 291 agent concentration, they are included in the pre-contrast relaxation rate, R * 2,pre . The 292 rate R * 2,κ,micro depends on the local chemical composition. The effects are independent 293 of the pulse sequence. Gadolinium-based contrast agent molecules increase the 294 relaxation rate, which can be described by a linear relationship [25], where r 2 is the molar relaxivity, and c κ the local contrast agent concentration in 296 compartment κ. The molar T 2 relaxivity, r 2 , of Gadobutrol at 3 T and 37 • C is 297 approximately 3.9 m 3 mol −1 s −1 [53]. Here, we assume that the contrast agent cannot 298 enter the cells, c S = 0, hence R * 2,S,micro = 0.

299
The term R * 2,meso stems from a meso-scale effect. The magnetic field perturbations 300 induced by the difference in magnetic susceptibility in the blood vessel and the 301 extra-vascular space, increase the relaxation rate of the extra-vascular space in 302 proximity of a blood vessel. The generated magnetic field perturbations are several 303 orders of magnitude smaller than B 0 . Furthermore, the influence decays rapidly with 304 distance to the vessel surface. Therefore, we consider each segment of the vessel network 305 to cause a perturbation independent of the other segments. The increase in R * 2 at a 306 given location in the tissue caused by mesoscopic magnetic field perturbation will then 307 be the superposition of all n segment perturbations where |c B −c I | is the difference of the concentrations averaged on the vessel surface. The 309 factor κ B ≥ 0 is an ad-hoc parameter, scaling the strength of these perturbations. The 310 proportionality factor ϕ i models the decay of the influence of the with distance from the 311 vessel wall. We set ϕ i = R 2 /r 2 , assuming a quadratic decay, where r is the distance to 312 the vessel center line and R the radius of the vessel segment. The susceptibility contrast 313 likewise increases the transversal relaxation rate, which we model by The same effect occurs at the cell surfaces, induced by the difference in magnetic 315 susceptibility between interstitial space and cells. Note that we consider cells not to be 316 invaded by contrast agent. We include this effect by adding a term to Eq. (22), and to the relaxation rate of the cell compartment, where κ T ≥ 0 is a second ad-hoc parameter, determining the strength of these 319 perturbations. Furthermore, we assume that there is no direct interface between the 320 cells and the vascular compartment.

321
Combining Eqs. (19), (21) and (23) to (25), we obtain a formulation for the transversal relaxation rate dependent on the concentration fields and the volume fractions of the three compartments: Longitudial relaxation with contrast agent administration 322 Similar to T * 2 , the contrast agent also shortens T 1 . However, the effects occur merely on 323 the micro-scale. Thus, we can model the relaxation rate R 1 = 1/T 1 of the tissue sample 324 where we implicitly assumed that contrast agent does not enter cells, c S (x, t) = 0. The 326 molar T 1 relaxivity, r 1 , of Gadobutrol at 3 T and 37 • C is approximately The relaxation rates, R * 2 and R 1 (Eqs. (26) and (27) block-matrix structure in residual form, where u is the respective discrete primary variable (fluid pressure in Eq. (15) where ILU 0 (A) denotes an incomplete LU-factorization of the matrix A using A's 356 sparsity pattern (zero fill-in) [54,Chapter 10]. 357 We assume that the influence of the sub-voxel contrast agent evolution during a 358 single image acquisition on the NMR signal is negligible, and thus, Eq. (17) is solved as 359 a post-processing step after each time step of the perfusion model. inherently patient-specific and cannot be directly measured, or parameters for which the 388 measurement data is not available for the given patient. These parameters are, a, b, t p , 389 κ B , κ T , T 1,pre , T * 2,pre , L p , D ω . Determining these parameters for a given signal-time 390 curve constitutes an inverse problem. In particular, we aim to determine D ω , which 391 may quantify contrast agent leakage, and thus, has direct clinical relevance.

392
In the following, we briefly discuss typical values or value ranges for these cf. Fig. 1, as well as the corresponding ||E opt || 2 , are given in Table 1.

431
A comparison of the simulated and measured NMR signals, Fig. 4, indicates that the 432 model can reproduce the measured curves well. Table 1 shows that the diffusive wall    Simulated normalized NMR signals compared with MRI data (see Fig. 1), using the best fit parameter estimates given in Table 1. Left -the result for the lesion sample (L), right -the result for the NAWM sample (N).  Table 1, this means that there is little to no contrast agent leakage for 450 sample N, while there is significant leakage for sample L. This is in accordance with the 451 present understanding of the pathology, which assumes leaky vessel walls in MS lesions. 452 However, the problem of finding best fit parameters is typically ill-conditioned, or neglected, the early time signal enhancement due to T 1 -shortening becomes even 486 stronger than the signal decrease due to T * 2 -shortening, as clearly seen in Fig. 7. This 487 illustrates that it is essential for the NMR signal model to include meso-scale effects.

488
The scaling parameter κ T for the meso-scale T * 2 -effects from the cell walls, only where p(θ|X) is the posterior distribution, i.e. the probability of θ given the observation 519 data X. p(X|θ) is the likelihood function, i.e. the probability of the X being from the 520 same population as the model prediction, given θ. p(θ) is the prior distribution 521 reflecting prior knowledge about the parameters θ, before knowing the observations. 522 p(X) is the marginal likelihood, a normalization constant, not depending on θ. Now, let 523 Y = M(θ) be the model prediction given the the parameters θ. We assume that we can 524 where is the combination of measurement error and unbiased model error and σ its 526 standard deviation. The likelihood that any model answer, Y , comes from the same 527 population as the measurement, X, is a Gaussian likelihood if the errors of all observations are assumed to be uncorrelated. The standard deviation, 529 σ, has to be estimated for the given MRI data and the proposed model. We assume 530 that our model represents the underlying physical processes accurately, so that model 531 and discretization error are negligible compared to the MRI data measurement error.

532
The measurement error is estimated from the MRI data obtained before the contrast 533 agent bolus reaches the tissue sample, where the measurement is assumed to fluctuate 534 around a constant baseline signal. To this end, we take 100 random signal samples from 535 the brain slice shown in Fig. 1, normalize the signal to the mean of the first 10 sample 536 data points, and compute the standard deviation of all such baseline data points across 537 all samples, yielding σ = 0.009. forward model runs within one step. The algorithm is briefly described in S2 Appendix. 548 We refer to the literature [63,64] for a comprehensive discussion.

549
In the following, Bayesian parameter inference is used to compute the probability  Table 2. The ensemble sampler is configured with k = 100 554 walkers. The parameter vector is θ = [a, b, t p , log 10 D ω , T 1,pre , κ B , κ T ] T , so that N = 7. 555 The parameter L p remains fixed to reduce the dimension of the parameter space. Its 556 influence on the NMR signal has been shown in the previous section to be significantly 557 weaker than the influence of D ω (see Fig. 6).

558
The sampler convergence is estimated using the integrated auto-correlation time, 559 τ f [63], is a finite chain of length M , e.g. the value of parameter a for each 561 sample in the Markov chain, and µ f its arithmetic mean. We use an estimate of the integrated auto-correlation, τ f,e , using the Python module acor [65,66]. We compute 563 this estimate for the chain of each parameter, θ i , and use the minimum and maximum 564 values, τ max = max 0≤i<N τ θi,e , τ min = min 0≤i<N τ θi,e . The sampler is run until the 565 sample size, j > 100 · τ max , and the change in the auto-correlation time estimate from 566 sample j − τ max to sample j is less than 1 %. The resulting histograms for each 567 parameter and their covariance with respect to the other parameters is visualized in 568 Fig. 8 for sample L and Fig. 9 for sample N (cf. Fig. 1). To eliminate artifacts from the 569 burn-in phase of the MCMC algorithm, the first 10 · τ max samples are discarded. To  Table 1 that were obtained previously with PEST.

573
To interpret the results, we recall the original question: What can we learn about the 574 model parameters, given the MRI data?

575
If the posterior distribution of a parameter is close to uniform, i.e. close to the prior 576 distribution (see Table 2), the data did not provide any additional information about 577 this parameter. This is the case for κ T and T 1,pre in Fig. 9, which is consistent with the 578 observation in Fig. 7 that the sensitivity of the NMR curve with respect to changes in 579 κ T or T 1,pre is low.

580
In contrast, if the posterior distribution differs significantly from the prior 581 distribution, the data provides significant information on this parameter. This is the 582 case for the parameters D ω and κ T in Fig. 8, and D ω and κ B in Fig. 9. Again, this is  Fig. 1, sample L). The histograms on the diagonal are the histograms for single parameters, the scatter plot in the matrix shows the covariance between the respective row and column parameters (plot generated with [67]  improved to include re-circulation and to be derived from AIF measurements.

614
The presented model considers processes in a sub-voxel tissue sample that is 615 surrounded by tissue with the same properties. However, contrast-enhancing lesions in 616 the brain typically span over several MRI voxels, see Fig. 1 understanding of the sub-voxel processes beyond the scope of this paper.
December 20, 2018 34/45 The relative error of the physical quantity u ∈ {p t , x t }, between reference and a solution 669 on a coarser grid, u ref , u, respectively, is defined as In time, we define the maximum relative error over all time steps t i , i ∈ {0, · · · , τ }, 671 where τ is the number of time steps, as 672 e u,∞ = max e i u .
Finally, we measure the difference of the signal-time curve, S, to the reference curve,

673
S ref , computed with the finest spatial and temporal discretization, in the following norm 674 The convergence rates for a given error e are computed from one refinement level n 675 to the next as 676 rate = ln e n+1 − ln e n ln ν n+1 max − ln ν n max , where ν max is the respective maximum discretization length. In space, ν max is defined as 677 the maximum edge length of all elements, h max . When refining, the vessel domain grid 678 is also refined by bisecting large elements until the maximum element length is smaller 679 than h max . In time, ν max is defined as the maximum time step size. The time step size, 680 ∆t, is chosen to be small around the time where the contrast agent front reaches the 681 domain, and increasingly larger as the process becomes slower, following the heuristic 682 ∆t = θ ln(t + 1.05), where θ > 0 is a factor controlling the time step size in the refinement study.

683
The reference solution is obtained with h max = 1 µm and θ = 0.125. The parameters 684 are chosen to be the optimal parameter set computed by an optimization algorithm 685 described in the following section Parameter estimation (see Table 1 Table 3. Errors and convergence rates in space for the pressure, p t , and the mole fraction of the contrast agent, x t , in the extra-vascular domain.   signal difference to the MRI data from an MS lesion shown in Fig. 1 (in red).
687 Table 3 show the errors and convergence rates of the extra-vascular fluid pressure, p t , 688 and the contrast agent mole fraction, x t . Fig. 11 shows the NMR signal curves and 689 errors with respect to the reference solution when refining in space and time.

690
It can be seen that all quantities converge to the reference solution. We obtain 691 convergence rates close to 1 for the pressure and the mole fraction of the contrast agent. 692 The signal curve converges with first order in time and a slightly higher order in space. 693 The higher convergence may be explained by the computation of the signal involving 694 the integration of the concentration over the entire domain. The relative error with 695 respect to the reference solution, is smaller than 1 % for a moderate spatial and temporal refinement. In conclusion, we consider a spatial resolution of h max = 8 µm, 697 and a temporal resolution θ = 1 as sufficient for the subsequent analysis. We justify this 698 with the assumption that the errors resulting from model parameter uncertainty, as well 699 as the errors in the measurement data, are larger than the discretization error. This is 700 also evident, when looking at the results of the parameter study and comparing the 701 variability with that of the signal-time curves shown in Fig. 11 for different spatial and 702 temporal discretizations. In order to verify that the discretization error is small also for 703 other parameter configurations, we ran the above analysis for various parameter where θ j is a walker position randomly drawn from the positions of the other set of 712 walkers, S 1−m , and ζ is a random variable drawn from a proposal distribution g(ζ), , where N = dim(θ) is the dimension of the parameter space. If the sample is not