Realistic numerical simulations of diffusion tensor cardiovascular magnetic resonance: The effects of perfusion and membrane permeability

Purpose To study the sensitivity of diffusion tensor cardiovascular magnetic resonance (DT‐CMR) to microvascular perfusion and changes in cell permeability. Methods Monte Carlo (MC) random walk simulations in the myocardium have been performed to simulate self‐diffusion of water molecules in histology‐based media with varying extracellular volume fraction (ECV) and permeable membranes. The effect of microvascular perfusion on simulations of the DT‐CMR signal has been incorporated by adding the contribution of particles traveling through an anisotropic capillary network to the diffusion signal. The simulations have been performed considering three pulse sequences with clinical gradient strengths: monopolar stimulated echo acquisition mode (STEAM), monopolar pulsed‐gradient spin echo (PGSE), and second‐order motion‐compensated spin echo (MCSE). Results Reducing ECV intensifies the diffusion restriction and incorporating membrane permeability reduces the anisotropy of the diffusion tensor. Widening the intercapillary velocity distribution results in increased measured diffusion along the cardiomyocytes long axis when the capillary networks are anisotropic. Perfusion amplifies the mean diffusivity for STEAM while the opposite is observed for short diffusion encoding time sequences (PGSE and MCSE). Conclusion The effect of perfusion on the measured diffusion tensor is reduced using an increased reference b‐value. Our results pave the way for characterization of the response of DT‐CMR to microstructural changes underlying cardiac pathology and highlight the higher sensitivity of STEAM to permeability and microvascular circulation due to its longer diffusion encoding time.


INTRODUCTION
Diffusion Tensor Imaging is an in-vivo non-invasive MR technique for characterizing the underlying tissue microstructure. This MR technique measures the self-diffusion of water molecules within a voxel providing information on the magnitude of diffusion, the degree of diffusion anisotropy and its orientation. 1,2 Diffusion tensor imaging applied to the heart presents several challenges such as the intrinsic deformation during the cardiac cycle or the respiratory associated movement. 3,4 However, recent developments 4,5 in diffusion tensor cardiovascular magnetic resonance (DT-CMR) have overcome these difficulties providing in-vivo insights into the microstructural abnormalities underlying various myocardial pathologies. 6,7 Typical DT-CMR sequences allow water molecules to diffuse tens of micrometers probing the local tissue microstructure. However, the MR signal is measured in a much larger voxel region (3 × 3 × 8 mm 3 ) effectively averaging the measured tissue microstructure. Furthermore, the link between a given change in microstructure and the corresponding change in DT-CMR parameters is not well established. The comprehension of this missing link between the MR image data and the microscopic structure lends itself to investigations using Monte Carlo (MC) random walk simulations [8][9][10] or finite element method (FEM) based simulators, 11,12 which have the power to elucidate the sensitivity of the DT-CMR signal to confounding factors or biological parameters present in the tissue microstructure.
The first cardiac microstructure numerical phantom 8,9 utilized MC simulations and simplified cardiomyocytes to cylindrical geometries with differing hexagonal cross sections and lengths. Recently, 10 we generated a more realistic numerical phantom based on manually segmented histology. In this virtual tissue, the membranes were considered impermeable and radially adjacent blocks were rotated in the local coordinate system to simulate the helix angle (HA) present in the myocardium. 13 By performing MC random walk simulations, we showed that the distribution of extracellular space (ECS) has a measurable impact on the DT-CMR parameters. We also showed that, in the case of permeable membranes, the distances diffused in long diffusion time sequences such as Stimulated Echo Acquisition Mode (STEAM) would be sufficient to span multiple cardiomyocytes. Simultaneously, experimental studies have shown that cardiomyocytes are not fully impermeable and exchange rates for intra-extracellular water exchange have been reported from 6 to 30 Hz. [14][15][16] Mathematical models 17,18 can be utilized in the MC random walk simulations to accurately represent the cell permeability by solving the interface boundary condition. 19 However, these models suffer from numerical limitations 20 when the membrane separates media with different diffusivites.
DT-CMR techniques are not only sensitive to molecular diffusion but also to the blood flow in capillaries. 21 Numerical simulations 22,23 have been performed in the brain where the capillary network is modeled as isotropic and several theories have been proposed for characterizing the microcirculation in the capillary network. 21 Due to the assumption that the capillaries are randomly oriented, most of these models do not translate well to cardiac and skeletal muscle tissue as the capillary network is anisotropic and, on average, parallel to the myocytes. 24,25 More recent numerical simulations [26][27][28] have created anisotropic capillary network models to simulate the effect of several perfusion parameters on the diffusion weighted MR signal.
The aim of this work is to study the sensitivity of DT-CMR parameters to changes in cell membrane permeability and microvascular perfusion. The simulations of the diffusion process use a new permeability model 20 that improves upon the treatment of the interface boundary condition 19 between regions of different diffusivity as compared to previous numerical simulations of diffusion tensor imaging. In a similar manner, the microvascular perfusion signal has been modeled assuming that a group of particles travel at constant velocities through a capillary network anisotropically oriented according to a Dimroth-Watson distribution. 29,30 The sensitivity of the tensor to perfusion has been studied for typical DT-CMR sequences by adding a weighted contribution of the perfusion signal to the simulated diffusion signal, providing a more complete MC random walk simulation than in previous work.

METHODS
Here we introduce the histology-based substrate, the treatment of the interface boundary condition substrate, the permeability model and the key elements of the MC random walk simulations. Next, we present the perfusion model and its incorporation into our DT-CMR simulations.

Histology-based substrate
The histology-based substrate is constructed by replicating a block of tissue of 495 × 392 × 127 m 3 . 10 This repeating block of tissue is generated by extruding 1182 contours obtained from manual segmentation of porcine histology (Masson staining). 31 The cardiomyocyte cross-sections were segmented from a region of interest in the mesocardium where the imaging plane was approximately Illustration of the two-dimensional histology-image and the region of interest (ROI) that was manually segmented. (B) Extrusion of the ROI in the long axis of the cardiomyocytes recreating a block of tissue of 500 × 400 × 120 m 3 . (C) Representation of the rotation applied to adjacent blocks in order to simulate the transmural variation in the helix angle (HA). Note that for illustration purposes, the three-dimensional model contains few block repetitions with fewer cardiomyocytes than the actual histology-image and the spacing between the blocks is not realistic. perpendicular to the cardiomyocyte long axis. As observed in Figure 1, each block of tissue is rotated around the myocardial radial direction (Y-direction in Figure 1) by applying a rotation of 10 • mm −1 to account for the transmural HA variation observed in the myocardium. 4,13,32 The different blocks are also shifted half a block every other row in the long axis of the cardiomyocytes (Z-direction in Figure 1) to eliminate long straight channels.

MC random walk
The Brownian motion or self-diffusion of water molecules can be modeled as a MC random walk of massless, noninteracting particles. Particles are uniformly seeded in the histology-based substrate inside a cuboid representing the dimensions of a typical DT-CMR voxel (2800 × 2800 × 8000 m 3 ) and a buffer zone to minimize biases due to particles entering/leaving the voxel during the simulation. 10 Note that a typical DT-CMR voxel has a larger slice thickness (Z-direction) than the in-plane resolution to maintain an acceptable signal-to-noise ratio during the short acquisition times. However, additional simulations have shown similar tensor results for an isotropic voxel (refer to  Table S1). At each time step t, the position of each particle ⃗ X(t) is updated using a random constant step-size vector ⃗ R, such that for each component R i , there is the same equal probability of moving one unit in either possible direction (R i = +1 or R i = −1). After one whole step, where ⃗ R is scaled by the local diffusivity √ 2D t which can either be intracellular, D ICS , or extracellular, D ECS . The intersections are resolved using additional substeps and the transit model described in Section 2.3. The cumulative phase for each particle is defined by where ⃗ G(t) is the time-dependent gradient vector. At the end of the random walk, the diffusion signal attenuation is described by the following equation where S diff ∕S 0 is the normalized diffusion signal, ⃗ e g the gradient direction, b a specific sequence b-value 33 and N p,diff the number of particles utilized in the diffusion MC random walk simulation. To validate our in-house MC random walk algorithm we compared to a finite element (FE) solution in a representative two-dimensional geometry (refer to Figures S1 and S2).

Permeability
In this work we assume uniform sarcolemma permeability sarco for all cardiomyocytes. A range of sarco values is considered, to study the effect of intra-extracellular water exchange on DT-CMR parameters, from impermeable sarco = 0 m ms −1 to highly permeable sarco = 0.05 m ms −1 . 14-16 Whenever a particle interacts with a cell membrane, a probability of transit p t is computed via a transit model that represents the permeability embedded in the membrane boundary condition. 19,20 This transit model, that we refer to as the hybrid model, 20 treats the membrane permeability and the step change in diffusion as successive barrier interactions. Doing so reduces the complexity of the transit probabilities to respectively, for transit from ECS to ICS and viceversa, where p b is the membrane transit probability and p d the probability of a particle transitioning between two media of different diffusivities. The membrane transit probability p b is derived by Fieremans et al. 17 and the probability of transitioning between two media p d of different diffusivity is computed through a step-balanced interface permeability. 34 These two probabilities are given by the following equations where x is the distance at the beginning of the step between the particle and the membrane surface projected onto the vector normal to the surface.

Extracellular volume fraction
The effect of changes in extracellular volume fraction (ECV) is studied by shrinking/growing the manually segmented cardiomyocytes. 10 The cardiomyocytes have been morphed obtaining a default healthy value of 24.69% and two extreme ECV values of 41.82% and 18.8% (see Figure S3).

Intercalated disks
Intercalated disks (ICD) are complex structures connecting the ends of neighboring cardiomyocytes and are composed of three major complexes: desmosomes, fascia adherens and gap junctions. 38 Desmosomes and fascia adherens reinforce the cardiomyocyte junction while gap junctions ensure the passage of electrical current between neighboring cardiomyocytes through intercellular channels. 39 To study the sensitivity of the measured diffusion tensor to possible water restriction in the ICDs, we compare simulations for ICDs having the same permeability as the rest of the sarcolemma, with simulations for lower ICD permeability of ICD = 0.005 m ms −140 on the end caps of each cardiomyocyte. As depicted in Figure 2, this low permeability region has been extended 2 m along the long axis of the cardiomyocyte to ensure that the entire end cap is within the region of low permeability irrespective of any geometry irregularity.

Perfusion model
The microvascular circulation during and between diffusion encoding gradient waveforms causes alterations to the diffusion MR signal magnitude. The contribution of microvascular perfusion to the MR signal is investigated independently from diffusion simulations by simulating a group of particles following specific paths determined by an anisotropic model. As depicted in Figure 2, each particle travels at constant velocity v through different multiple straight segments traversing k segments, each with a length l k and orientation ⃗ e x k . The multiple segments are derived from an anisotropic model that simulates the orientation of the capillary network providing a simplified representation of the capillary orientation without the need for a detailed network. Note that the mean capillary orientation is always aligned with the mean orientation of the cardiomyocytes, thereby following the same HA rotation described in Section 2.1. Similarly to random walk, the phase acquired by any particle is given by equation (2). For constant intercapillary velocity, the perfusion phase perfusion can be simplified by summing the accumulated phase of each traversed segment. Following this simplification, perfusion is normalized by the b-value, gradient strength G max and velocity v 27,41 and is then only dependent on the orientation of the segments ⃗ e x k and the gradient profile ⃗ G(t) normalized in amplitude and time h(t∕T), that is, ⃗ F I G U R E 2 (A) 3D illustration of the capillary network with representative expanded (spacing increased for visualization) cardiomyocytes, showing a Dimroth-Watson distribution (upper PDF) for the zenith angle z k and Weibull distribution (lower PDF) for the length L k of the capillary segments. The PDF has been calculated considering parameters that are representative of the capillary network in the myocardium. (K = 3.25 and = 60 m, = 40 m for the length distribution 24 ). (B) 3D illustration of the intercalated disk (ICD) structures between the neighboring cardiomyocytes. A more detailed illustration for the two different permeabilities ( sarco and ICD ) is depicted.
In the above s k−1 and s k are the normalized times when the particle starts and completes a certain capillary segment, m 0 = ∫ s 0 h(s ′ ) ds ′ the normalized 0-th gradient moment and a the normalized b-value a = √ b∕( G max T 3∕2 ). At the end of the simulation, the normalized perfusion signal S perfusion (b, ⃗ e g )∕S 0 is determined through Equation (3). The total normalized signal (S(b, ⃗ e g )∕S 0 ) is defined by the weighted contribution: where ⃗ e g is the gradient direction and f the perfusion fraction (i.e., the proportion of the detectable spins within the capillary network). The total signal loss is defined as 1 − S(b, ⃗ e g )∕S 0 and has been obtained using a healthy ECV of 24.7% for the diffusion simulations with sarco = ICD = 0.03 m ms −1 and a default perfusion fraction of f = 10% as echocardiographic studies 42 report a volume fraction of blood in myocardium of 12.6% ± 1.8% for healthy individuals at rest. Figure 3 illustrates the intercapillary Gaussian velocity distributions used in the simulations. Increased standard deviation v corresponds to greater velocity dispersion, while minimal v = 0.001 mm s −1 approximately models constant intercapillary velocity. A mean capillary velocity of v = 0.5 mm s −143 has been assumed while the

F I G U R E 3
Histogram showing six Gaussian intercapillary velocity distributions for N p,perf = 10 5 perfusing spins. The Gaussian distributions are all centered around a constant mean intercapillary velocity v = 0.5 mm s −1 .
intercapillary velocity distributions have been truncated at −0.1 and 1 mm s −1 to avoid extreme velocities.
The perfusing spins (N p,perf ) have been randomly seeded on one side of the voxel at plane z = 0 and the perfusion simulations have been performed with N p,perf = 10 5 which we determined sufficient to reliably compute the perfusion signal (see Figure S4). In the perfusion model, particle displacements are exclusively influenced by the anisotropic model and rotation within the voxel, eliminating the need for a buffer zone.
The effect of perfusion on the diffusion tensor D D has been studied by computing a new tensor D D+P fitting a mono-exponential model to the total normalized signal calculated in Equation (7).
Two different tensor fittings have been considered, with b ref values of 0 and 0.15 ms m −2 , respectively. Additional simulations (see Table S3) suggested that using six encoding directions was enough to reliably calculate tensor parameters. Consequently, six diffusion encoding directions 44

Capillary network
Diffusion studies traditionally 21 assumed an isotropic capillary segment distribution. However, the capillaries in skeletal muscle and in the myocardium are highly anisotropic with a preferred orientation aligned with the myocyte long axis. 24 Several experimental studies have quantified the anisotropy in skeletal 29,46 and cardiac muscle 30 capillaries and showed that a hemispherical Dimroth-Watson axial distribution 47,48 provides a good fit to capillary segment orientation in both muscle types. On this basis, a Dimroth-Watson distribution has been used here for the zenith angle z between each capillary segment and the cardiomyocyte long axis. The joint probability density of p( , z ) can be expressed as 48 and integrating with respect to the azimuth angle ( ) reduces the equation to: where K is a variable that determines the degree of anisotropy and may vary from isotropic K = 0 to highly anisotropic K ≫ 0. Note that due to the use of spherical coordinates the isotropic case K = 0 distributes more points around the equator z = ∕2 than closer to the cardiomyocytes axis z = 0. Equation (10)

The effects of ECV and variations in permeability
The diffusion tensor is studied for changes in the ECV and sarcolemma permeability ( sarco ). In Figures 4 and 5, the fractional anisotropy (FA), the mean diffusivity (MD) and each eigenvalue of the diffusion tensor are plotted for three different sequences (STEAM, PGSE, and MCSE). Increasing effects of restrictions on diffusion are observed for all three sequences when reducing the ECV from 41.8% to 18.8%. Conversely, increasing permeability reduces the anisotropy of the diffusion tensor, weakening FA and the 2 ∕ 3 ratio. Quantitatively, for the impermeable case sarco = 0 m ms −1 , this ECV reduction results in a decrease in MD of 40.9% (STEAM), 22.7% (PGSE) and 17.6% (MCSE). For ECV = 24.7%, FA decreases from 0.7 to 0.5 as sarco rises from 0 to 0.03 m ms −1 for STEAM. The large MD reduction as ECV decreases observed for STEAM compared with the other sequences is mostly a consequence of large reductions in 2 and 3 . Considering the same impermeable case, 2 only decreases 18.7% (PGSE) and 14.6% (MCSE) while 2 is reduced by 50% for STEAM. High 2 ∕ 3 ratios are observed for low ECV values when using STEAM and PGSE while this tendency is not present for MCSE as 2 ∕ 3 ratios are similar between the three ECV values.

F I G U R E 4
Fractional anisotropy, mean diffusivity, and 1 parameters of the diffusion tensor D D are illustrated considering three different extracellular volume fraction geometries (41.8%, 24.7%, 18.8%) and a range of realistic sarcolemma permeability values for cardiomyocytes. The graph shows mean values using six simulations with N p = 10 5 and N t = 10 4 walkers and timesteps respectively with a 95% confidence interval.

F I G U R E 5
2 , 3 and 2 / 3 are plotted for three different extracellular volume fraction geometries (41.8%, 24.7%, and 18.8%) and a range of realistic sarcolemma permeability values for cardiomyocytes. The plot illustrates mean values using six simulations with N p = 10 5 and N t = 10 4 walkers with the 95% confidence interval shown by the error bars. permeability-related increases in the diffusivity in the direction of the primary diffusion eigenvector 1 for PGSE and MCSE. However, for STEAM, we observe a larger 1 increase of ≈ 34% for all ECV cases as sarco rises from 0 to 0.01 m ms −1 . A more gradual increase in 2 , 3 when increasing permeability is noted for STEAM. Quantitatively, there is an increase of 156% for 2 (ECV = 24.7%) and 168% for 3 (ECV = 24.7%) when going from impermeable to the maximum permeability. These deviations are reduced to 29%, 25% and 8%, 7% for PGSE and MCSE, respectively. In the diffusion simulations, narrow uncertainty angle cones are obtained for all ECV and permeabilities values ranging from 0 • to 1 • for ⃗ E1 and 3 • to 6 • for ⃗ E2, ⃗ E3 suggesting that the eigenvectors are well-defined for all three sequences (see Table S4).

ICD effect
The effect of ICD permeability is investigated by introducing a lower permeability on the endcaps of the simulated cardiomyocytes. Figure 6 illustrates the effect of reduced ICD permeability ( ICD ) compared to the results obtained when considering a homogenous permeability across all membranes of the cardiomyocytes. The lower ICD permeability predominantly reduces the primary eigenvalue 1 (7% reduction for sarco = 0.05 m ms −1 ) which in turn reduces the MD for STEAM. This effect is only apparent when sarco ≠ 0 m ms −1 and only for the STEAM sequence. No major changes in FA, 2 and 3 related to ICD permeability are observed for any of the three sequences (see Figure S5). Figure 7 depicts the normalized phase Φ perfusion described in Equation (6) when the diffusion encoding gradient direction is aligned with each of the three diffusion tensor eigenvectors ( ⃗ E1, ⃗ E2, ⃗ E3). A narrow phase distribution is noted for STEAM and PGSE when the gradient direction is parallel to the main capillary direction ( ⃗ E1, which is also the mean cardiomyocyte orientation) while encoding along the two orthogonal directions ( ⃗ E2, ⃗ E3) results in wider symmetric distributions. Encoding along the main capillary direction with minimal inter-capillary velocity dispersion v = 0.001 mm s −1 results in the maximum total phase (i.e., all magnetization vectors point in a similar direction) and the highest normalized perfusion signal S perfusion (b, ⃗ e g )∕S 0 , while encoding along the orthogonal directions leads to signal loss reducing the normalized perfusion signal. For MCSE, the normalized phase distributions show a peak centered at 0 for all three gradient directions indicating relative insensitivity to mean

F I G U R E 6
Illustration of the effect of intercalated disks (ICD) in three different sequences illustrating fractional anisotropy, mean diffusivity and 1 . For this plot, an extracellular volume fraction = 24.7% and a sarco = 0.05 m ms −1 have been considered. The different plots show mean values using six diffusion simulations with N p,diff = 10 5 and N t = 10 4 walkers with a 95% confidence interval. Note that we have considered ICD = 0.005 m ms −1 for all sarco cases except for the impermeable case where ICD = sarco = 0 m ms −1 has been considered.

F I G U R E 7
Top figures: The relative particle frequency for the normalized phase Φ perfusion is plotted for three different gradient directions parallel to the eigenvectors of the diffusion tensor D. Bottom figures: The signal attenuation along ⃗ E1 is plotted for the three sequences considering several Gaussian intercapillary velocities varying its SD from v = 0.001mm s −1 to v = 0.15 mm s −1 .
capillary flow velocity changes in the intercapillary velocity dispersion while the narrow width of the same Φ perfusion distributions result in minimal perfusion signal loss in all directions. For STEAM and PGSE, increasing the intercapillary velocity dispersion v >> 0 broadens the phase perfusion distribution resulting in reduced normalized perfusion signals. As an example, using STEAM the normalized perfusion signal (S perfusion (b, ⃗ E1)∕S 0 ) in the direction of ⃗ E1 considering a b-value of 0.15ms m −2 is reduced from 0.951 to 0.234 when the velocity dispersion is increased from v = 0.001 mm s −1 to v = 0.15 mm s −1 . Figure 8 illustrates the effect of the capillary anisotropy K and the mean intercapillary velocity v on the normalized perfusion signal measured in the direction parallel to the cardiomyocyte long axis. Even with a spherical distribution, the isotropic regime K = 0 results in different perfusion signal values for both velocity distributions v . Increasing the anisotropy of the capillary network K > 0 increases the normalized perfusion signal and results in almost no signal loss when the particles travel at a constant velocity v = 0.001 mm s −1 while using a higher intercapillary velocity dispersion reduces the normalized perfusion signal. Note that these signal differences are greater for long diffusion time sequences. Quantitatively, the normalized perfusion signal at K = 3.25 is reduced from 0.91 to 0.053 and from 0.98 to 0.93, respectively, for STEAM and PGSE when increasing the inter-capillary velocity dispersion v from 0.001 to 0.15. STEAM and PGSE are less sensitive to changes in v than in s v and K and small deviations in perfusion signal are only observed when the intercapillary velocity distribution is near constant ( v = 0.001 mm s −1 ). Similarly to the results in Figure 7, the normalized signal from MCSE displays minimal changes in response to the changes in mean inter-capillary velocity v, the velocity dispersion v or the anisotropy K of the capillary network. . The MD, 1,2,3 , and FA are linearly reduced as the perfusion fraction increases for PGSE and MCSE while STEAM presents higher apparent diffusivity in all diffusion directions (see Figure S6). The eigenvalue reductions observed in both short diffusion encoding sequences are the same for both tensor calculations and both extreme inter-capillary velocity distributions v . Quantitatively, for both v and tensor calculations, MD decreases ≈ 23% for PGSE and ≈ 35.6% for MCSE when f rises from 0% to 20%. Despite these large changes in MD, the FA maintains a value close to the diffusion only result (f = 0%).

Perfusion and diffusion
Due to a longer diffusion time, the effect of perfusion when increasing f (%) is greater on STEAM than PGSE or MCSE amplifying all the eigenvalues of the tensor 1,2,3 . For STEAM, a distinction between the two extreme intercapillary velocity distributions v is observed in the tensor parameters assessed. A small velocity dispersion v = 0.001 mm s −1 only increases 2,3 , maintaining the same diffusivity along the cardiomyocytes long axis direction ( 1 ) and reducing the overall tensor FA. On the contrary, as observed in Figure 7, increasing v for STEAM yields maximum signal loss toward ⃗ E1 intensifying the overall measured eigenvalues in all three principal directions resulting in higher MD and lower FA values. Conversely to PGSE and MCSE, larger variations are observed between the tensor fittings with the two different ref-  Figure 10 investigates the effect of varying the anisotropy of the capillary network. Similarly to Figure 9, two extreme intercapillary velocity distributions v are

F I G U R E 9
Fractional anisotropy, mean diffusivity, and 1 parameters for several perfusion fraction values using a K = 3.25 (red vertical line in Figure 10 considered together with tensor fittings with two different reference b-values. Note that the differences observed for STEAM in the normalized perfusion signal when the capillary network is isotropic K = 0 in Figure 8 (top row) are highly reduced in 1 (see Figure 10) when adding the normalized diffusion signal. Tensor variations are mostly observed for STEAM and PGSE when 0 < K < 10, while differences in the tensor due to changes in v are only observed for STEAM. Specifically for this sequence, increasing v reduces the influence of K on 1 . Quantitatively, 1 decreases 21% for v = 0.001 mm s −1 and only 3% for v = 0.15 mm s −1 as K increases from 0 (isotropic network) to 50 (highly anisotropic network).

DISCUSSION
Simulations of the effects of cell membrane permeability and capillary perfusion on DT-CMR show a significant dependence on model parameters. Three sequences (STEAM, PGSE, and MCSE) were simulated using MC random walk 20 on a histology-based tissue microstructure 10 in which the manually segmented cardiomyocytes were morphed (shrunk and thickened) to represent a physiological range of ECVs (18.8% to 41.8%). While STEAM and MCSE are in vivo and ex vivo DT-CMR techniques, PGSE is only suitable for ex vivo but is included for completeness and scientific interest. The primary eigenvector of the diffusion tensor ( ⃗ E1) is aligned with the cardiomyocytes long axis 10 while the second and third eigenvector directions ( ⃗ E2, ⃗ E3) are contained in perpendicular plane to ⃗ E1, with ⃗ E2 lying in the sheetlet plane and ⃗ E3 defining the sheetlet structure as its normal. 31,49 Having 2 >> 3 and narrow cones of uncertainty suggests that there is a preferred diffusion direction in the perpendicular plane to ⃗ E1 due to a specific sheetlet geometry structure. The results indicate that STEAM and PGSE have more sensitivity to the sheetlet structure compared to MCSE as low ECV values are linked to higher 2 ∕ 3 ratios. For all permeability values, all sequences show a reduction in MD with decreasing ECV due to a lower diffusion coefficient in the intracellular space, D ICS < D ECS . PGSE and STEAM show higher FA values for low ECV fractions while MCSE demonstrates the opposite. Due to the short effective diffusion time of the MCSE sequence the ratio 2 ∕ 3 is the same for all the three extracellular volume geometries and the anisotropy of the tensor (FA) is mostly driven by 1 .

F I G U R E 10
Fractional anisotropy, mean diffusivity, and 1 parameters considering a constant perfusion fraction of f = 10% (red vertical line in Figure 9) for a range of different capillary networks ranging from an isotropic (K = 0) to highly anisotropic (K > 10) capillary distributions toward the long cardiomyocytes axis. Two extreme intercapillary velocities are depicted considering two distinct fitting tensors methods with b ref = 0 ms m −2 and b ref = 0.15 ms m −2 . Note that the red line depicts a realistic anisotropic value K = 3.25 for the capillary network observed in the myocardium and the gray area illustrates extreme anisotropic cases in a different horizontal axis scale.
The diffusion restriction can be reduced by including permeability in the diffusion simulations allowing the walkers to cross between compartments. The sensitivity of each sequence to permeability is related to the comparison between the diffusion time (Δ) and the exchange time. The exchange time ( exchange ) is directly proportional to ECV and inversely proportional to the permeability. 50 Small permeability-related deviations in DT-CMR parameters are observed for PGSE and MCSE. However, major deviations are observed for STEAM, particularly at low permeability values. A permeability of sarco = 0.01 m ms −2 is linked to a water exchange time of 96 ms for the highest ECV. Note that for STEAM (Δ ≈ (1000 ms)) there is ample time for the walkers to cross multiple membranes between ICS and ECS, in contrast to the short diffusion time of both PGSE and MCSE. This explains the sudden jump in 1 , MD, and FA when increasing permeability sarco from 0 m to 0.01 m ms −1 , whereas a more gradual transition is observed for PGSE and MCSE as Δ < exchange .
ICDs connect neighboring cardiomyocytes and enable the passage of ions between cells allowing cardiac action potential to spread along the cardiac muscle. 51 Abnormalities of impulse propagation, arrhythmic tendency and hypertrophied myocardium have been linked to significant changes in gap junctions. 52 To study the sensitivity of the diffusion tensor to the presence of the ICD, results are compared for constant sarcolemma permeability with results for lowered ICD permeability 12 ( ICD ). Among all the sequences evaluated only STEAM shows changes in the diffusion tensor in response to reduced ICD permeability, as the spins diffuse further for larger Δ and are more likely to collide with the cardiomyocyte end caps. Short diffusion time sequences limit the number of collisons with the ICD zones reducing any observable effect.
The perfusion simulations have been performed by recreating an anisotropic capillary network with the mean capillary orientation aligned with the cardiomyocytes long axis ( ⃗ E1). The perfusion simulations are performed separately from the diffusion simulations but considering the same HA rotation within the voxel. The simplicity of the perfusion model permits computing a phase for each particle that is normalized (Φ perfusion ) by the b-value, the constant capillary velocity and the overall sequence time. The signal at the end of the simulation is obtained by summing the contributions of each particle's magnetization vector which are modified by the phase distribution perfusion .
In order to simulate a distribution of velocities within the capillary network, a Gaussian intercapillary velocity distribution has been assumed with SDs between v = 0.001 mm s −1 and v = 0.15 mm s −1 . Note that if v = 0.001 mm s −1 the Gaussian bell is narrow, approximating a constant velocity which leads to similarly narrow distribution for ( perfusion ) and the normalized phase Φ perfusion . When the diffusion encoding gradient is parallel to ⃗ E1 and v = 0.001 mm s −1 , due to the similarity between perfusion and Φ perfusion , there is minimal perfusion signal loss as the phase dispersion within a voxel is low, that is, all the magnetization vectors point in a similar direction. This is in line with the results reported in a similar perfusion model. 27 However, the findings 27 of minimal signal loss due to perfusion in the direction of ⃗ E1 contradict experimental studies 26 and other numerical simulations 28 which suggest that the direction ⃗ E1 is related to the maximum signal loss.
The present results demonstrate the importance of considering variable velocities within the capillary network as this spreads and changes the normalized phase Φ perfusion distribution. The perfusion simulations show that incrementing the intercapillary velocity dispersion from v = 0.001 mm s −1 to v = 0.15 mm s −1 yields a faster perfusion signal reduction with b-value in the ⃗ E1 direction and shifts the maximum signal decrease toward this direction. These results are in line with a numerical study that reported faster signal reduction for laminar flows compared to plug flows for the same capillary networks. 26 As shown in Figure 7, large intercapillary velocity dispersions v = 0.15 mm s −1 reduce the effect of the anisotropy of capillary network or the mean magnitude of the intercapillary velocity to changes in perfusion signal. Conversely, a capillary network based on the original intravoxel incoherent motion method 21 (isotropic) is less susceptible to changes in the intercapillary velocity dispersion. Nevertheless, small variations in anisotropy (K ≈ 2) in the capillary network when using long diffusion time sequences present large deviations from the isotropic case. We argue therefore that it is necessary to consider both dispersion in orientation and in velocity. However, the present perfusion model is limited to constant intracapillary velocities and thus MCSE shows minimal sensitivity to perfusion.
The effect of perfusion on DT-CMR has been simulated by fitting a tensor D D+P to a weighted signal from both the diffusion and perfusion processes. The impact of perfusion on the tensor D D+P depends on the perfusion fraction f (%) and the capillary network geometry. The tensor D D+P has been computed considering a b-value of 0.6 ms m −2 and two different b ref values of 0 ms m −2 and 0.15 ms m −2 . Long diffusion time sequences (STEAM) allow water molecules to perfuse and disperse large distances, increasing the sensitivity of the perfusion signal to changes in v while less sensitivity is observed for short diffusion time sequences (PGSE and MCSE). As demonstrated in Figure 7, for STEAM when the gradient is parallel to ⃗ E1, increasing v results in a rapid decline in the normalized perfusion signal as b-value increases. When the b-value is fixed and v is large, the signal loss from perfusion in STEAM is greater than that of diffusion, lowering the total signal resulting in increased eigenvalues and MD. Conversely, the total signal for PGSE and MCSE increases because the signal loss associated with diffusion is greater than that of perfusion. This results in a counter-intuitive reduction in MD as the perfusion fraction (f ) increases as the proportion of signal originating from perfusing spins increases in comparison to diffusing spins.
For STEAM, an increased reference b-value b ref = 0.15 ms m −2 reduces the effect of perfusion leading toward similar FA values for both D D+P and D D tensors. Similarly an increase in the reference b-value in experimental studies 53 is thought to reduce the effect of perfusion. When using STEAM, increasing the inter-capillary velocity dispersion ( v ) not only accelerates the perfusion signal decay with b-value as described above, but reduces the contribution of the anisotropy of the capillary network. This can be observed in Figure 10 for changes in 1 where STEAM and v = 0.15 mm s −1 show moderate changes in the tensor (of the order of 3% for 1 ) between K = 0 and K ≫ 0 while larger deviations are observed for v = 0.001 mm s −1 .

CONCLUSION
The separate and combined effects of permeability and microvascular perfusion on DT-CMR measurements have been investigated with the extent of possible confounding contributions due to altered diffusion at the ICD. The models developed in this work provide new insights into the response of DT-CMR to microstructural changes underlying cardiac pathology, which will be of great relevance in interpreting changes observed in clinical DT-CMR studies.
In conclusion, ECV and permeability have measurable impacts on DT-CMR parameters. Increasing permeability reduces the exchange time between the water molecules located in the intracellular and ECS, augmenting the effective diffusion. The ratio of the water exchange and diffusion times determines the sensitivity of a certain sequence to changes in permeability. Of the sequences considered, results show that STEAM is the most sensitive to permeability variations and the only sequence with measurable variations in 1 in response to changes in the permeability of the ICD. These findings suggest that STEAM is a potential sequence for the early detection of certain permeability-related pathologies and abnormalities of impulse propagation.
With regard to microvascular perfusion, it is concluded that the intercapillary velocity dispersion ( v ) has a major impact on the perfusion signal when capillary networks are anisotropic. Long diffusion time sequences (STEAM) are more sensitive to changes in v as water molecules perfuse further during the diffusion time and over more tortuous paths. Our results show that STEAM yields the highest perfusion signal loss toward the primary diffusion eigenvector ⃗ E1 for the maximum intercapillary velocity dispersion considered. In contrast, PGSE and MCSE are less sensitive to perfusion leading to a more pronounced diffusion signal loss compared to perfusion. This results in an unexpected decrease in MD and all eigenvalues in D D+P as the perfusion fraction is increased for PGSE and MCSE sequences.
For all sequences but specifically for STEAM, it is shown that that using a tensor fitting with a b ref > 0 ms m −2 reduces the effect of perfusion when comparing the diffusion tensor D D parameters to D D+P . For clinical DT-CMR studies where changes in permeability may be of interest, STEAM should be considered the sequence of choice. However, the sensitivity of STEAM to microvascular perfusion must also be considered and these confounding effects can be at least partially mitigated by the use of a higher reference b-value.

SUPPORTING INFORMATION
Additional supporting information may be found in the online version of the article at the publisher's website. Figure S1. Comparison between Finite Element solution and Random Walk in a representative domain for sarco = 0.05 m ms −2 and T=500 ms. The FE solution was obtained considering Continuous Galerkin (CG) elements while assuming a small buffer region accounting for the membrane. Figure S2. Comparison between Finite Element solution and Random Walk in a representative domain for sarco = 0 m ms −2 and T=500 ms. The FE solution was obtained considering Continuous Galerkin (CG) elements while assuming a small buffer region accounting for the membrane. Figure S3. Illustration of three different extracellular volume fractions (ECV) obtained by morphing the manually segmented cardiomyocytes. Figure S4. Convergence of the perfusion signal obtained for a b-value of 0.6 ms ms −2 for different number of perfusion particles N p,perf . Figure S5. Illustration of the effect of ICD in three different sequences illustrating 2 and 3 . For this plot, an ECV = 24.7% and a sarco = 0.05 m ms −1 have been considered. The different plots show mean values using 6 diffusion simulations with N p,diff = 10 5 and N t = 10 4 walkers with a 95% confidence interval. Note that we have considered ICD = 0.005 m ms −2 for all sarco cases except for the impermeable case where ICD = sarco = 0 m ms −2 has been considered.  E1 are plotted for the three sequences considering several Gaussian inter-capillary velocities varying its SD from v = 0.001 mm s −1 to v = 0.15 mm s −1 . Table S1. Comparison of diffusion tensor parameters acquired from 6 diffusion simulations using an anisotropic and isotropic voxel. The diffusion simulations are performed considering STEAM, ECV = 24.7%, and sarco = 0.02 m ms −2 . The results show similar results using the anisotropic and isotropic voxels which can be attributed to the similarity in microstructure in the slice direction. Table S2. Sequence parameters for the three typical DT-CMR sequences considered in the simulations.