Simulation of Contrast Agent Transport in Arteries with Multilayer Arterial Wall: Impact of Arterial Transmural Transport on the Bolus Delay and Dispersion

One assumption of DSC-MRI is that the injected contrast agent is kept totally intravascular and the arterial wall is impermeable to contrast agent. The assumption is unreal for such small contrast agent as Gd-DTPA can leak into the arterial wall. To investigate whether the unreal assumption is valid for the estimation of the delay and dispersion of the contrast agent bolus, we simulated flow and Gd-DTPA transport in a model with multilayer arterial wall and analyzed the bolus delay and dispersion qualified by mean vascular transit time (MVTT) and the variance of the vascular transport function. Factors that may affect Gd-DTPA transport hence the delay and dispersion were further investigated, such as integrity of endothelium and disturbed flow. The results revealed that arterial transmural transport would slightly affect MVTT and moderately increase the variance. In addition, although the integrity of endothelium can significantly affect the accumulation of contrast agent in the arterial wall, it had small effects on the bolus delay and dispersion. However, the disturbed flow would significantly increase both MVTT and the variance. In conclusion, arterial transmural transport may have a small effect on the bolus delay and dispersion when compared to the flow pattern in the artery.


Introduction
Dynamic susceptibility contrast magnetic resonance imaging (DSC-MRI) has been shown to be a powerful technique to qualify cerebral blood flow and is playing an increasing role in diagnosis of acute ischemic stroke [1,2]. However, the accuracy of perfusion parameters obtained from DSC-MRI is challenged by the delay and dispersion of the contrast agent bolus caused by the estimated arterial input function (AIF) [3,4].
Cerebral blood flow quantification requires knowledge of the arterial input function that is the concentration of contrast agent in the feeding vessel to the tissue of interest. In theory the AIF should be measured at the feeding vessel close to the tissue of interest [5,6]. However, due to technical difficulties, it is commonly estimated from a major artery (e.g., the middle cerebral artery) for analysis of cerebral perfusion [5]. Delay and dispersion of the contrast agent may occur between the position of the AIF estimated and the tissue of interest, which would lead to underestimation of the blood flow [5].
Another implicit assumption of the technique of DSC-MRI related to the estimation of AIF is that no or a negligible amount contrast agent penetrates the arterial wall, due to the fact that, in general, the vascular wall is considered irrelevant for the contrast mechanism in DSC-MRI. The assumption is tenable for the intact blood brain barrier (BBB) [7]. However, for a damaged BBB which usually occurs in the stroke and normal arterial wall except in the brain, the assumption is unreal. The commonly used contrast agent gadolinium diethylenetriominepentax-acetic acid (Gd-DTPA) with hydrodynamic diameter of less than 2 nm can easily leak into the arterial wall, since it has been demonstrated that the hydrophilic solute with the diameter less than 7.0 nm can be transported through the intercellular junction of endothelial cell into arterial wall [8]. Whether 2 The Scientific World Journal the unreal assumption is valid for the estimation of the delay and dispersion of the contrast agent bolus has never been verified. It seems that the leakage of contrast agent into the arterial wall may affect contrast agent transport, which might cause distortion of the contrast concentration-time curve and hence the delay and dispersion of the contrast agent bolus.
To investigate the delay and dispersion effects of a bolus of contrast agent, two approaches were usually applied. One that has been commonly used was convolving the estimated AIF with a vascular transport function (VTF). The VTFs were quite differently assumed, ranging from the simple model of a single-exponential to the more sophisticated model of a feeding artery in series with small parallel vessels [3,9]. The disadvantage of the method is that accuracy of the different VTF models is hard to assess. The other was the simulation of the blood flow and contrast agent transport using computational fluid dynamics (CFD) to describe the delay and dispersion errors of the contrast agent bolus [10][11][12]. The approach has been used to investigate the effect of carotid artery stenosis on the cerebral blood flow quantification and to simulate the dispersion along a simplified coronary artery with different stenosis under steady and unsteady flow condition. Although these studies indicate that CFD simulations were useful to investigate the delay and dispersion of the contrast agent bolus, all the studies were based on the unverified assumption of impermeable arterial wall.
In the present study, to investigate whether the unreal basic assumption of DSC-MRI is valid for the estimation of the delay and dispersion of the contrast agent bolus, we formulated the arterial wall as a five-layer model and numerically simulated the flow and the transport of contrast agent in the model using CFD. This five-layer model included the endothelial glycocalyx layer (EGL), the endothelium, the intima, the internal elastic lamina (IEL), and the media, which were all treated as macroscopically homogeneous porous media. The effects of different factors that may affect contrast agent transport such as the integrity of endothelium and the disturbed flow after the stenosis of artery on the delay and dispersion of contrast agent were further analyzed.

Geometry of the Model.
As a basic geometrical model, the arterial segment concerned was simplified as a straight axisymmetric cylinder. The inner diameter of model ( ) and the wall thickness of the artery ( ) are chosen to be = 3.7 mm and = 0.34 mm [13,14]. The wall of the arterial segment was modeled as a five-layer structure and the thickness of each wall layer is shown in Figure 1(b) [15,16], where the ratio of intima and media thickness was chosen to be approximately 0.75 [17].
For the stenosed model, the variation of the vessel radius along the stenosis was described using a cosine function and the reduction in the cross-sectional area of the lumen was set as 75% [18].

Governing Equations
2.2.1. Lumen. The flow simulation in the lumen of the arterial segment was based on the incompressible Navier-Stokes and continuity equations: where u and represent, respectively, the fluid velocity vector and the pressure. and are the density and viscosity of blood ( = 1050 kg⋅m −3 , = 3.45 × 10 −3 kg⋅m −1 s −1 ). The mass transport of contrast agent Gd-DTPA in the flowing blood can be described by where is the concentration of Gd-DTPA, and is the diffusion coefficient of Gd-DTPA in blood. At 37 ∘ C, when assuming that the radius of Gd-DTPA ( ) is 1 nm [19], can be calculated as 6.5847 × 10 −11 m 2 s −1 from the following Stokes-Einstein equation [20]: where is Boltzmann constant, and is the absolute temperature.

Arterial Wall Layers.
The transmural flow across the arterial wall can be described by the Brinkman equation as follows [16]: where u and represent, respectively, the superficial velocity vector and the pressure based on the volume averaged method. and are the porosity and the hydraulic permeability of the wall layer concerned. For all layers, the viscosity of plasma was assumed to be 0.72 × 10 −3 kg m −1 s −1 [16]. The transport of Gd-DTPA across the arterial layers was modeled by the following equation [16]: where is the superficial concentration of Gd-DTPA; is the effective diffusivity of Gd-DTPA in the arterial layers; is the filtration reflection coefficient. because it is hard to measure the parameters, the transport parameters are scarce. In the present study, most of the parameters were obtained from theoretical models. The calculation of the parameters is presented in Appendix. The parameters used are summarized in Table 1. Figure 1  flow velocity profile was used for the pulsatile flow simulation (Figure 1(c)) [21]. The mean velocity was chosen to be 0.24 m s −1 [22].

Boundary Conditions. As shown in
BC-B the pressure at the outlet boundary of the artery lumen was set at 100 mm Hg.
BC-C at the media-adventitia interface, a constant pressure boundary condition with 30 mm Hg was employed.
BC-D symmetric conditions were set at the axis of symmetry.
BC-E no viscous flow was set on both the axial ends of the arterial wall.
The boundary conditions for the mass transport equations ( (2) and (6)) are as follows.
0 was set to 3 s, which meant that the bolus of contrast was injected after the periodicity of the flow field was completely developed.
BC-2: for other boundaries, zero concentration gradient was assumed.

Computation
Procedures. The numerical simulations were carried out using a validated finite element algorithm COMSOL Multiphysics. Three full cardiac cycles (3 s) simulation of the pulsatile flow with a time step of 4 ms were carried out to achieve a periodic flow independent of the initialization. Based on the initial velocity field, the mass transport equations were solved coupling with flow transport equations for 37 s with time step from 1 ms to 4 ms depending on the temporal variations of contrast concentration.

Quantification of Dispersion.
The effect of the dispersion of the contrast agent bolus on the AIF can be described by convolving the estimated AIF with a vascular transport function (VTF) [4,23]: where 0 is the bolus delay and is distance to the inlet. Using the variance of the VTF around its mean, the degree of dispersion can be quantified: where the mean vascular transit time (MVTT) is given by the ratio of the first to the zeroth moment of the VTF and can be used to quantify the bolus delay: AIF est is prescribed directly from the contrast agent concentration at the inlet ((7) ( , = 0)) and AIF is computed as the area weighted average of the concentration of contrast agent on several lumen cross-sections between the inlet and the outlet perpendicular to the axial vessel direction ( ( , )).
To obtain the mean delay (MVTT) and the variance ( 2 VTF ), it is not necessary to calculate VTF firstly, since they can be calculated directly with the zeroth, first, and second integral moment of the concentration of the contrast agent [4]: The moments of the concentration of the contrast agent are about 0.5 second and decreases 2%, when compared with that at the inlet ( Figure 2). The delay and dispersion are enlarged in the arterial wall, resulting in that the maximal concentration at 20 in endothelial layer is about 1.3 second later and 18% lower than that at the 10 ( Figure 3). In addition to the delay and dispersion in the axial direction, as demonstrated in Figure 2, the concentration of the contrast in the radial direction is much more decreased and spread. For instance, the maximal concentration of intima layer is approximately 20 times smaller than that of the endothelial layer and the time interval between the half maximal and maximal concentration at 20 at the intima layer is approximately 4 times longer than that at the endothelial layer ( Figure 3). The development of the bolus delay and dispersion in the unstenosed model are quantified by the mean vascular transit time (MVTT) and the variance of the VTF from (11). As illustrated in Figure 4(a), MVTT and the variance are increased almost linearly in relation to the distance to the inlet, which reflects consistently the concentrationtime curve in Figure 2. When compared to the simulation with contrast agent transport only in the lumen, the arterial transmural transport would lead to slight increase in the MVTT and small increase in the variance (Figures 4(a)-4(b)) resulting in only increase of 0.186 s 2 at 10 and 0.288 s 2 at 20 . The quantitative difference in variance in Figure 4(c) further demonstrates that although the increase in the variance would increase gradually in relation to the distance to the inlet, the increase ratio would decrease sharply resulting in approximately 25% at 20 .

Effects of Endothelium on Contrast Agent Transport and the Bolus Dispersion.
When the endothelium is damaged in diseased condition, the transport of the contrast agent may be affected and hence the bolus dispersion. To investigate the role of endothelium in the transport of the contrast agent, the endothelium was assumed to be totally damaged and the four transport parameters of the damaged endothelium and EGL were simplified to be similar to that of the intima. In addition, the transport parameters of other layers were consistent with that of the model with intact endothelium.
Figures 5(a)-5(d) demonstrate the concentration-time curve of the damaged endothelial model in the four arterial layers. As evident from these figures, due to the damage of the endothelium, the contrast agent is sharply accumulated in the arterial wall.
As illustrated in Figures 5(e)-5(f), the difference in the MVTT between the damaged endothelial model and intact one is slight and the difference in the variance between the two models decreases along the axial direction resulting in almost no difference at the outlet. These results indicate that the damage to the integrity of endothelium would have a small effect on the bolus delay and dispersion.

Effects of Disturbed Flow on Contrast Agent Transport and the Bolus Dispersion.
Until now, the numerical simulations have been carried out only for a simplified straight axisymmetric blood vessel. However, the geometry of the physiological blood vessel is much more complex with such characteristics as branching, twisting, taper, and curvature, which would lead to far more than parabolic flow profile in the simplified model but very complicated flow patterns [24,25]. One common flow pattern developed in the arterial system is the disturbed flow, which is inclined to occur in locations such as aneurysm, stenosed artery, the inner wall of curved segments, and the outer walls of arterial bifurcations [26,27]. To further investigate the effects of arterial transmural transport on the bolus delay and dispersion under disturbed flow, contrast agent transport in the stenosed model was calculated. Figure 6(a) shows the concentration distribution in the stenosed region at different moments of the contrast injection. It can be observed that the concentration in the disturbed flow region after the stenosis is quite uneven and significantly later than that of the regions of high flow.
Figures 6(b)-6(e) demonstrate the concentration-time curve in the four arterial layers at the disturbed flow region (10 ). The concentration in the disturbed flow model decreases and spreads much more significantly than that in the straight model at the same axial location.
In the stenosed model the influence of the disturbed flow on the bolus delay and dispersion was observed immediately behind the stenosis (Figure 7). After the initial significant increase in the MVTT and variance behind the stenosis a reduction in the MVTT and variance follows in the disturbed flow zone. The arterial transmural transport would also affect the bolus delay and dispersion by decreasing the MVTT and increasing the variance. However, the effect of the arterial transmural transport on the MVTT and the variance is much less than that of the disturbed flow caused by the stenosis.

Discussion
One unreal assumption of DSC-MRI is that the arteries are impermeable to contrast agent. In the present study, we numerically coupled contrast agent transport in the arterial wall with that in arterial lumen to investigate whether the unreal assumption is valid for the estimation of the delay and dispersion of the contrast agent bolus. The results obtained reveal that the arterial transmural transport would slightly affect the bolus delay qualified by MVTT and moderately increase the bolus delay qualified by the variance. The MVTT and the variance induced by the arterial transmural transport are much less than that by the disturbed flow after the stenosis. In addition, although the integrity of endothelium can significantly affect the accumulation of contrast agent in the arterial wall, it has small effects on the bolus delay and dispersion. The bolus dispersion would lead to a systematic blood flow underestimation. Graafen et al. simulated the dispersion in coronary arteries using a computational fluid dynamics approach and demonstrated that the variance between 1.0 and 2.5 s 2 would lead to the myocardial blood flow underestimation between about 6% and 10% [11]. The present study demonstrated that the variance in both the stenosed and the unstenosed model is less than 2.5 s 2 and the variance caused by the arterial transmural transport is less than 0.5 s 2 , which indicated that the myocardial blood flow error induced by the arterial transmural transport may be well below 10%.
As the transport of contrast agent is mainly governed by the convective transport in the lumen and the diffusion transport in lumen and the arterial wall, the local flow pattern, the flow rate, and the arterial transmural transport would affect the degree of bolus dispersion. Among the three factors, the arterial transmural transport may be the least important one, since the present simulation indicated that the disturbed flow produced much more dispersion of the bolus than the arterial transmural transport. However, the contribution of local flow pattern and the flow rate to the dispersion is comparable. Our simulations demonstrated the disturbed flow induced by the stenosis would significantly increase the bolus dispersion. In contrast, Graafen et al. indicated that stenosis in the coronary arterial model leads to a reduction of dispersion [11]. Calamante et al. found that for the two patients with similar degree carotid stenosis, only one of them had larger dispersion, while the other was found to produce a dispersion of the bolus similar to that found in the normal subjects [10]. These studies and our results indicated that, to achieve an accurate estimation of the dispersion, the simulations with patient-specific geometrical and boundary layer model should be performed and the effects of the arterial transmural transport may be neglected in practice.
This study revealed that local flow pattern caused by the stenosis would significantly affect the local contrast agent transport and hence the local bolus dispersion, which is constant with mass transport simulations based on the patient-specific model. For instance, Liu et al. demonstrated that the disturbed flow would significantly hinder the mass transport in the human aorta and computational studies indicated that the oxygen transport was significantly low in the outer wall of carotid artery where disturbed flow developed [28,29]. Due to the fact that flow patterns are quite complex in the arterial segments such as the aneurysm, stenosed artery, the inner wall of curved segments, and the outer walls of arterial bifurcations [26,27], the present results combining previous studies indicated that the arterial input function may be estimated beyond these regions to avoid large bolus dispersion [12]. In this pilot study, the accuracy of the results may be reduced by the simplifications of the transport parameters, the blood rheological properties, and the geometries.
Due to the scarcity of the transport parameters, most of the parameters were obtained from theoretical models. The estimated parameters, especially the diffusivity of the contrast agent, would affect the simulation results. The diffusivity used in the present study (6.5847 × 10 −11 m 2 s −1 ) is lower than the estimated diffusion coefficient of contrast agent in the previous studies (5.5 × 10 −10 m 2 s −1 and 1.5 × 10 −10 m 2 s −1 ) [11,12]. The lower diffusion coefficient in the lumen would lead to an increase of the dispersion due to a smaller exchange of contrast agent particles perpendicular to the flow direction. However, the lower diffusion coefficient in the arterial wall The Scientific World Journal reduces the arterial transmural transport and hence decreases the dispersion in the lumen, which is shown in the present result that when compared to the damaged endothelial model the normal one with relatively low diffusivity leads to a small dispersion.
The blood in the present study was simplified as Newtonian fluid and the non-Newtonian rheological properties such as shear thinning nature of the blood were neglected. It is reported that the shear thinning non-Newtonian nature of blood could slightly reduce oxygen flux (similar micromolecule as contrast agent) in most regions of the arteries, and this effect became much more apparent in areas with disturbed flow [29]. Therefore, it may be valid to assume the blood Newtonian fluid to estimate the bolus dispersion in most regions of the model However, the assumption of Newtonian fluid may lead to an underestimation of dispersion in the disturbed flow region.
Another limitation of the present study is that all simulations were performed on simplified geometries. The geometry of the physiological blood vessel is very complex with such characteristics as branching, twisting, taper, and curvature, which would lead to complicated flow patterns. Therefore, it is necessary to use the realistic geometries to simulate the exact dispersion of the contrast agent bolus. As the main aim of the present study is not to estimate the exact dispersion errors in an individual patient, the results obtained from the typical parabolic flow and disturbed flow can still shed some light on the effects of arterial transmural transport of the contrast agent on the bolus delay and dispersion.

Conclusion
The arterial transmural transport would slightly affect the bolus delay and small increase the bolus dispersion, which may have small effects on the estimation of the blood flow estimation. In comparison with disturbed flow induced by the presence of stenosis, the arterial transmural transport plays a less important role in the bolus delay. Although the assumption of impermeable arterial wall of DSC-MRI is unreal, it may still be a good simplification for the estimation of the delay and dispersion of the contrast agent bolus.

A. Parameters of the Model
A.1. Parameters of the EGL. The EGL layer of the arterial wall was assumed to be composed of symmetrical EGL fibers and leaky junctions of cells resulting from either dying or being in mitosis. According to the results by Liu et al. [15], the porosity and hydraulic permeability ( egl ) of EGL fibers were obtained as 0.6735 and 6.0383 × 10 −18 m 2 .
The total filtration reflection coefficient of the EGL layer is determined by (A.1) The reflection coefficient of the fibers ( egl ) is given by [30]  where is the ratio of the contrast agent Gd-DTPA radius ( = 1 nm) to the outer radius of the fluid annulus = 3 1/4 (2 + Δ)/ √ 2 ( = / ) and = / . is the fiber radius (6 nm). Δ is the space between fibers (8 nm).
In (A.1), the reflection coefficient ( ), the hydraulic permeability ( ), and the effective diffusivity ( ) of the leaky junctions are determined using the pore theory as follows [20]: where / is the portion of leaky junctions on surface area and is the hydraulic permeability of one leaky junction ( = 2 /3). If we assumed that the leaky cell is located at the center of each periodic circular unit of radius [31], then / can be determined by (2 cell )(2 )/( 2 ). Since the fraction of the leaky junctions, , which is defined as the ratio of the area of leaky cells to the area of all cells ( = 2 cell / 2 ), can be determined from experimental results, (A.4) can be given by The total effective diffusivity of the EGL is calculated from where is the effective diffusivity of leaky junction: = . (A.7) The effective diffusivity of a single leaky junction is determined by [20] = free (1 − ) ( ) , where free is the Gd-DTPA diffusivity in the free plasma, which is 3.1537 × 10 −11 m 2 s −1 at 37 ∘ C calculated from (3) with viscosity of . The effective diffusivity ( egl ) for the EGL is obtained from the stochastic theory [32]: After taking cell as 15 m and setting as 0.0005 [33], egl and egl are determined as 0.0555 and 1.0128 × 10 −10 m 2 s −1 , respectively.
A.2. Parameters of the Endothelium. Filtration flow can move through the endothelium via both normal junctions and leaky junctions of the endothelial cells. Therefore, the hydraulic permeability ( end ) and the reflection coefficient where end is the thickness of the endothelium and = 1.576 × 10 −9 m s −1 mm Hg −1 according to the experimental results [34].
The reflection coefficient of the normal junction is defined as is the Gd-DTPA diffusivity in the normal junctions, which can be determined as 3.02 × 10 −13 m 2 s −1 according to the radius of Gd-DTPA [35]. From (A.10) and (A.7), the hydraulic permeability, the reflection coefficient, and the effective diffusivity of the endothelium can be calculated as 1.7383 × 10 −20 m 2 , 0.1212, and 3.0276 × 10 −13 m 2 s −1 , respectively.
A.3. Parameters of the Intima. Using the results from Dabagh et al., the porosity and hydraulic permeability ( int ) of the intima are assumed to be 0.8025 and 0.42 × 10 −16 m 2 [36]. Applying the same method by Liu et al. [15], the effective diffusivity and the reflection of the intima can be obtained as 1.2219 × 10 −10 m 2 s −1 and 0.2514, respectively.

A.4. Parameters of the IEL.
The IEL is assumed to have a constant thickness with fenestral pores. It seems that the fenestral pores were filled with the same fiber matrix as intima [37,38]. Thus we took the same properties of the intima for the IEL to calculate the parameters and assumed that the fenestral pores were approximated as cylinder pores with a radius ( fen ) of 1.5 × 10 −7 m [33]. The hydraulic permeability ( IEL ), the reflection coefficient ( IEL ), and the effective diffusivity of IEL ( IEL ) are then determined as follows: where is the permeability of the porous medium and is the fractional void volume of smooth muscle cells in the media. equals 6.09 × 10 −19 m 2 after taking and as 1.43 × 10 −18 m 2 and 0.4.
The effective diffusivity of the media is obtained from the following equation: