Animal handling
All studies were conducted according to a protocol approved by the Institutional Animal Care and Use Committee. Female wild-type (WT) C57BL/6 J mice (n = 7) were obtained from Charles River Laboratories (Sulzfeld, Germany) and were studied at an age of 12 weeks. Mice were anesthetized with 4% isoflurane in 2.0 Vol.% oxygen (2 L/min), applied by a nose cone and were positioned vertically (head up). For cardiac and respiratory monitoring during the measurement, a pressure- sensitive pneumatic balloon (Graseby Medical Limited, Watford, United Kingdom) was placed between the inner radio frequency (RF) resonator wall and the murine thorax. The pressure signal from the balloon was transformed into an electrical signal by a pressure transducer (24PCEFA6 D, Honeywell S&C, Golden Valley, Minnesota, USA) and was amplified and processed in real-time by a custom-built ECG unit [17]. Due to the small inner diameter of the gradient insert and the RF coil, core body temperature could be maintained at physiological 37 °C during CMR-measurements by adjusting the temperature of the gradient cooling unit.
Data acquisition
Measurements were performed with a 17.6 T vertical-bore small animal MR system (Bruker Avance 750 WB, Bruker BioSpin MRI GmbH, Rheinstetten, Germany, operated with Paravision 4.0) with a 1 T/m gradient system (diameter: 40 mm) and a custom-built single-channel transmit-receive electromagnetic (TEM) resonator (inner diameter: 24 mm). To localize the position of the aortic arch, balloon-triggered axial and longitudinal 2D-cine FLASH measurements were acquired. Subsequently, retrospective flow measurements were performed with a radial PC-FLASH-sequence (see Fig. 1B) in a 3D slab perpendicular to the aorta (image volume: 25 × 25 × 4 mm3, see Fig. 1A). Spatial encoding was performed with a 3D radial trajectory with an angular density optimized for the anisotropic field of view [18], which was calculated with an open source tool box [19]. For flow-encoding a balanced 4-point flow-encoding scheme [20] with an encoding velocity of vENC = 125 cm/s was used. Each flow-encoding step consists of a read-out with 1.6 x 105 radial projections (140 read-out points, TR = 3 ms) covering a 3D-sphere in k-space (Fig. 1C). A flip angle of 15° was used in order to achieve high blood-tissue contrast. To minimize artifacts caused by off-resonances and signal dephasing due to accelerated flow, the echo time was set to 1.1 ms. In order to guarantee such a short echo time, it was necessary to design the amplitudes of the dephase gradients in a way that the gradient echo occurs at tE = 0.1 × tacq, where tacq is the acquisition time (Fig. 1B). To further increase robustness of data acquisition, the measurement was segmented into 10 subsets, each consisting of 1.6 × 104 projections (each covering a full 3D sphere in k-space) and 4 flow-encoding steps (Fig. 1D), which were acquired one at a time.
This kind of segmentation allows repetition of corrupted data sets in case of disturbances, e.g. caused by instabilities of the heart rate. Acquisition time of one subset was 3.2 min, leading to a total measurement time of 32 min for a full 4D flow protocol.
Phantom measurements
The stability of the 4D flow-encoding sequence was tested in a phantom consisting of a flow pump (MPC-Z V1.10, ISMATEC, Cole-Partner GmbH, Wertheim, Germany) with constant flow and adjustable flow values (max flow: 50 ml/s) and a silicone tube (ø = 6 mm). Flow was measured at 10 different flow values (15.00 ml/s - 26.25 ml/s in equidistant steps) with the protocol described above using only one subset per measurement (scan time: 3.2 min per subset) and the same encoding velocity as the in vivo measurements vENC = 125 cm/s. Mean flow values and standard deviations were calculated over 10 slices. To prevent artifacts due to phase aliasing, phase unwrapping was applied when necessary [21]. For comparison, flow was also quantified by measuring volumetrically in liters. The reference measurement was repeated 8 times and mean values and standard deviations were determined.
Self-navigation
All signal processing was performed with MATLAB (The Mathworks, Inc., Natick, Massachusetts, USA). For retrospective self-navigation, the magnitude value of the center k-space signal (k = 0) was used. First, high-frequency disturbances were removed by using a matched filter for low-pass filtering [22]. The matched filter can be interpreted as a convolution of the noisy navigator signal with a conjugated time-reversed small portion of the signal [23]. After filtering a baseline subtraction [24] was used in order to eliminate low-frequency modulations caused by respiratory motion and by the transient to the steady state.
Trigger points and breath gating intervals were determined with variable thresholds (Fig. 2A-C). Using a linear assignment, each read-out was afterwards allocated to a value between 0 and 1, corresponding to a phase in the cardiac cycle (Fig. 2D). For removal of corrupted data points due to respiratory motion, the time average of the trigger point intervals (i.e., the mean cardiac period) was calculated for all 40 subsets, respectively. Only trigger point intervals lying in a ± 4 × TR interval (±12 ms) window around the temporal average were accepted for reconstruction. For respiratory gating data points during inspiration were assigned to a cardiac phase value of −1. For reconstruction, the read-outs were divided into 20 selection intervals, corresponding to 20 cardiac phase intervals. For each selection interval, the associated projections were combined and an image was reconstructed using a nonuniform fast Fourier transform (NUFFT) with an open source software toolbox [25, 26]. In this manner a set of four 3D-cines (one flow-compensated cine and 3 flow-encoded cines) with 20 frames and an isotropic spatial resolution of 100 μm, respectively, was reconstructed.
Off-resonance correction
At high magnetic field strengths, B0 offsets and field gradients cause deviations of the radial trajectory, which can lead to severe blurring artifacts in the reconstructed images. To remove these artifacts, an additional flow compensated radial 3D FLASH measurement with two different echo times was performed in the same field of view (FOV) (tE1 = 1.3 ms, tE2 = 2.3 ms total measurement time: 3.2 min, spatial resolution: isotropic 100 μm). A 3D off-resonance map was calculated from the phase differences between the two images [27]:
$$ \Delta f\left(x,y,z\right)=\frac{\phi_2-{\phi}_1}{2\pi \cdot \left({t}_{E1}-{t}_{E2}\right)}, $$
(1)
where tE1,2 and φ1,2 are the echo times and phases of the two images. Assuming only a global field offset and constant field gradients and neglecting local field inhomogeneities, the spatial dependent off-resonance frequencies can be approximated as:
$$ \Delta f\left(x,y,z\right)\approx \Delta {f}_0+\alpha \cdot x+\beta \cdot y+\delta \cdot z. $$
(2)
To remove blurring artifacts induced by a global frequency offset ∆f0, the frequency value from the center of the FOV (x = y = z = 0) was taken and used for a phase correction of the MR signal:
$$ {S}^{\prime }(t)=S(t)\cdot {e}^{i2\pi \cdot \Delta {f}_0t}, $$
(3)
where S(t) denotes the raw uncorrected MR signal. To also correct deviations of the k-space trajectory caused by global field gradients, the corrected trajectory \( {\overrightarrow{k}}^{\prime }(t) \) used for re-gridding was calculated with:
$$ {\overrightarrow{k}}^{\prime }(t)=\overrightarrow{k}(t)+\left(\begin{array}{c}\alpha \\ {}\beta \\ {}\delta \end{array}\right)t, $$
(4)
where \( {\overrightarrow{k}}^{\prime }(t) \) denotes the undistorted radial trajectory. α, β and δ were determined through linear fits of the off-resonance map.
Rigid motion correction
Due to the vertical setup of the MR scanner, a slight slipping and shifting of the mouse was sometimes observed during the measurement. Since this would lead into undesired motion artifacts such as blurring and phase subtraction errors, a rigid motion correction was applied prior to the cine reconstructions.
First, the 10 measurement subsets described above were used for the reconstruction of 40 time-averaged 3D-images (4 encoders times 10 measurement subsets). Using the first image I1 as reference, the shifts x, y and z were calculated for each subsequent image In in order to minimize the error between these images (Fig. 3A):
$$ \Psi =\underset{x,y,z}{argmin}{\left\Vert {I}_1-{T}_{x,y,z}\cdot {I}_n\right\Vert}_2^2\kern2em n=2,3,...\mathrm{40.} $$
(5)
Hereby Tx,y,z denotes the translation operator in respect to the image coordinates x,y and z, which needs to be applied for minimization of the error between the first and the n − th image. The algorithm yields shift values on a time base of 48 s (Fig. 3B). This information was used for a phase correction of the signal in k-space using the Fourier Shift Theorem [28] prior to the reconstruction.
Image processing and segmentation
Depending on the slice orientation and the alignment of the aortic arch within the image volume, the phase accumulations induced by flow encoding can cause slight deformations of the waveforms of the self-gating signal. Due to these distortions, the cines of the 3 flow-encoders are sometimes temporally shifted against the flow-compensated cine. In order to correct these temporal shifts, the time-dependent image intensities averaged over one slice were compared against each other. By using cross-correlation [29], the temporal shifts were determined for each encoder and the cines were synchronized.
For segmentation of the aortic arch, an adapted version of the previously described semi-automatic segmentation technique [30] was used. This technique assumes that the segmentation of all 4 cines should in principle lead to the same number of identified pixels. By evaluating a cost function, an optimum threshold value corresponding to a minimum deviation between the flow encoding measurements can be derived. Using this technique, each slice (in z-direction) of the 3D cine was segmented independently. Slices near the aortic root were excluded due to strong signal cancellations induced by accelerated flow. Subsequently, the three velocity components (vx, vy, vz) were calculated from the phase differences between the cines. Using the segmentation data, velocity was afterwards zeroed outside the aorta and filtered with a spatial median filter with a 3-connectivity neighborhood inside the lumen [31]. The spatial median filter removes outliers of velocity values due to segmentation errors near the lumen boundaries but leaves velocity data within smooth regions inside the vessel untouched [7].
Calculation of WSS and OSI
Assuming a Newtonian and incompressible fluid, the general form of the WSS →τ can be written as [6]:
$$ \overrightarrow{\tau}=2\eta \overset{\cdot }{\varepsilon}\cdot \hat{n}, $$
(6)
where η denotes the viscosity of blood, \( \hat{n} \) the inward unit normal of the lumen surface and \( \overset{\cdot }{\varepsilon } \) the deformation tensor:
$$ {\overset{\cdot }{\varepsilon}}_{ij}=\frac{1}{2}\left(\frac{\partial {v}_j}{\partial {x}_i}+\frac{\partial {v}_i}{\partial {x}_j}\right),\kern2em i,j=1,2,3. $$
(7)
Hereby xi,j denotes the spatial coordinates and vi,j the velocity components.
To calculate the WSS and radial stress, the PC and segmentation data were imported into Ensight (CEI systems, USA). The velocity derivatives and the surface normals were calculated directly from the 3D velocities and the isosurface of the lumen segmentation using a custom-built Python script. For the blood viscosity, a value of η = 0.04 Pas was assumed [9].
To separate the two components of the WSS and the radial stress, a centerline of the aortic arch was calculated, as described by [7] (Fig. 4A). Ring segments perpendicular to the centerline were afterwards generated at 14 different locations of the aorta (Fig. 4B). The isosurface of the segmented lumen imported to Ensight consists of a grid with approx. 5 · 103 nodes. For each node of the surface grid, the longitudinal (\( \hat{l} \): in parallel to centerline), radial (\( \hat{r} \): pointing towards the centerline) and circumferential \( \left(\hat{c}=\hat{l}\times \hat{r}\ \right) \) unit vector was calculated. The WSS and radial stress can afterwards be separated using:
$$ {\tau}_l=\overrightarrow{\tau}\cdot \hat{l},\kern2em {\tau}_c=\overrightarrow{\tau}\cdot \hat{c},\kern2em {\tau}_r=\overrightarrow{\tau}\cdot \hat{r}. $$
(8)
Mean and median values of the three components were calculated for each ring segment and cardiac phase. In addition, temporal averaged WSS values \( \left(\overline{\tau \to \left(r,\to \right)}\right) \) were derived using:
$$ \overline{\tau \to \left(r,\to \right)}\mid =\frac{1}{T_{RR}}{\int}_0^{T_{RR}}\overrightarrow{\tau}\left(\overrightarrow{r},t\right) dt, $$
(9)
where \( \overrightarrow{\tau}\left(\overrightarrow{r},t\right)=\left[{\tau}_l(t)\kern0.5em {\tau}_c(t)\kern0.5em {\tau}_r(t)\right] \)are the time-dependent WSS and radial stress components and TRR is the cardiac period. To also measure the temporal variability of the WSS waveforms and the degree of oscillatory flow, the OSI was calculated from the time dependent stress values using [5]:
$$ \mathrm{OSI}=\frac{1}{2}\left(1-\frac{\mid {\int}_0^{T_{RR}}\overrightarrow{\tau}\left(\overrightarrow{r},t\right) dt\mid }{\int_0^{T_{RR}}\mid \overrightarrow{\tau}\left(\overrightarrow{r},t\right)\mid dt}\right). $$
(10)
No changes in the direction of the stress over time result in a minimal OSI value (OSI = 0). In contrast, when strong periodic variations and sign changes occur, e.g. caused by recirculative flow during the diastolic cardic phase, the integral value nears itself to the limit \( {\int}_0^{T_{RR}}\overrightarrow{\tau}(t) dt\to 0 \) and the OSI approximates its maximal value (OSI = 0.5).