MRI model-based non-invasive differential diagnosis in pulmonary hypertension

Pulmonary hypertension(PH) is a disorder characterised by increased mean pulmonary arterial pressure. Currently, the diagnosis of PH relies upon measurements taken during invasive right heart catheterisation (RHC). This paper describes a process to derive diagnostic parameters using only non-invasive methods based upon MRI imaging alone. Simultaneous measurements of main pulmonary artery (MPA) anatomy and flow are interpreted by 0D and 1D mathematical models, in order to infer the physiological status of the pulmonary circulation. Results are reported for 35 subjects, 27 of whomwere patients clinically investigated for PH and eight of whomwere healthy volunteers. The patients were divided into 3 sub-groups according to the severity of the disease state, one of which represented a negative diagnosis (NoPH), depending on the results of the clinical investigation, which included RHC and complementary MR imaging. Diagnostic indices are derived from two independent mathematical models, one based on the 1D wave equation and one based on an RCR Windkessel model. Using the first model it is shown that there is an increase in the ratio of the power in the reflected wave to that in the incident wave (Wpb/Wptotal) according to the classification of the disease state. Similarly, the second model shows an increase in the distal resistance with the disease status. The results of this pilot study demonstrate that there are statistically significant differences in the parameters derived from the proposed models depending on disease status, and thus suggest the potential for development of a non-invasive, image-based diagnostic test for pulmonary hypertension. & 2014 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/3.0/).


Introduction
The pulmonary circulatory system is a low-pressure, high-flow vascular bed (Crapo, 2004) with multiple branches, running from the right ventricle to the left atrium. Pulmonary hypertension (PH) is a disease of the pulmonary circulation represented by a complex haemodynamic and pathophysiological state (Galie et al., 2011), which is clinically defined by an increased mean pulmonary arterial pressure (mPAP), Z25 mmHg at rest, measured by right heart catheterisation (RHC), and has many causes .
General characteristics of PH include vascular remodelling, intimal proliferations of the vessels, vasoconstriction, and increased pulmonary arterial stiffness. Each of these phenomena contributes to an elevation of pulmonary vascular resistance (PVR) and a reduction of vessel compliance (C) that leads to the augmentation of right ventricular afterload and ultimately right ventricular failure and death (Saouti et al., 2010). The overall system resistance (R total ) and compliance (C) change with the disease state (Westerhof et al., 2009) and a 0D-distributed model of impedance could be used to estimate their values (Grant and Paradowski, 1987). In addition to the changes in the Windkessel parameters, it has also been demonstrated (Laskey et al., 1993), (Hollander et al., 2001), (Huez et al., 2004), (Dwyer et al., 2012) that the pathophysiological modifications introduced by PH are detectable in the proximal pulmonary vessels from the relative magnitudes of the waves reflected back from the increased distal resistance of the diseased peripheral vasculature.
The wave reflections can be quantified from synchronised flow (Q(t)) and pressure (p(t)) measurements at the same anatomical location. In previous publications (Furuno et al., 1991), (Hollander et al., 2001), (Dwyer et al., 2012) the two measurements were made independently, and there are spatio-temporal alignment issues that might introduce error into the calculation process. Furthermore the reported pressure measurements are invasive. This paper proposes a non-invasive method, based on synchronised magnetic resonance imaging (MRI) measurements of blood flow and wall motion in the pulmonary artery (PA), to characterise the pulmonary circulation in healthy volunteers and stratified PH.
The first model is a 3-element (R c CR d ) Windkessel ( Fig. 1.a). R c is commonly chosen to be the characteristic impedance, Z c of the upstream vessel, defined as the pressure-flow ratio in the absence of wave reflections (Nichols 2005), but in the current model it is a free parameter that simply represents the resistance that is proximal to the capacitor. R d represents the distal resistance of the pulmonary circulation and C is the compliance of the pulmonary vasculature. The sum of the two resistors equals the ratio of mean pressure to mean flow (Westerhof et al., 2009), which can be interpreted as the clinically measured PVR (Lankhaar et al., 2006). This simple model should capture the increase of right ventricular (RV) afterload determined by the PVR increase and compliance decrease (Saouti et al., 2010), and so might be used to differentiate between healthy and diseased systems (Reuben 1971), (Saouti et al., 2010). It has already been shown (Lankhaar et al., 2006) that characterisation of the parameters in this configuration, using invasive measurements, can separate patients without PH from those with chronic thromboembolic pulmonary hypertension (CTEPH) and idiopathic pulmonary arterial hypertension (IPAH).
The second model is based on the assumption that the pressure and flow are measured at a cross-section of a straight elastic tube ( Fig. 1.b). This system is characterised by forward-and backwardtravelling waves transmitted from the boundaries (Shi et al., 2011), (van de Vosse andStergiopulos, 2011) and the hypothesis is that PH might be differentiated by the relative magnitudes of the waves.

Patients
Images were acquired from 35 subjects, 8 healthy volunteers and 27 PH patients referred to the Sheffield Pulmonary Vascular Disease Unit. Each patient underwent RHC, based on which they were divided into 3 groups according to the value of the measured PVR, expressed in Wood Units (WU): PVR o 4WU (n¼ 9), PVR Z4WU (n¼ 14) and NoPH (n¼ 4). N.B. WU is the clinical unit of measurement of PVR, which is equivalent to mmHg s/ml. The latter group consisted of those patients with suspected PH but for whom the diagnosis was not confirmed by the measured pressure, mPAP o 25 mmHg. This study was undertaken following approval from the local research ethics committee.

MRI acquisition
MRI is able to provide spatio-temporally registered measurements of flow and area at the same physical location. The flow waveforms were derived from velocity sensitised phase contrast (PC) MR imaging. PCMR produces two images, representing the magnitude and phase of the signal, at each time point. The magnitude image is used for anatomical orientation whilst the phase image encodes the velocity of the fluid in every pixel (Donald et al., 2003). In principle the cross-sectional anatomy of the vessel can be measured from the magnitude image, but in practice the contrast is poor during diastole when there is minimal flow. To overcome this problem a balanced steady state free precession (bSSFP) cine sequence, with high vessel to blood delineation (Cukur and Nishimura, 2009), was used to extract the dynamic radius changes during the cardiac cycle. Typical MRI images of the main PA during peak systole and late diastole from PCMR and bSSFP sequences are presented in Fig. 2.
PCMR and bSSFP sequences acquired in the same plane were used to acquire images of the main pulmonary artery (MPA), orthogonal to the pulmonary trunk and approximately 2 cm distal to the pulmonary valve. In each sequence, using retrospective acquisition and ECG gating, 40 time points were captured, during breath-hold, on a GE HD Â 1.5T scanner with the an 8 channel cardiac receiver array coil. The PC and bSSFP images were spatially and temporarily synchronised, using the same pixel size (256 Â 128 matrix dimensions, 480 Â 288 mm 2 FOV) and same number of cardiac images (40). For the PC sequence, 150 cm/s was used for the velocity encoding, 5.85 ms repetition time (TR), 2.87 ms echo time (TE) and 10% arrhythmia rejection. For the bSSFP acquisition, the TR was 3.73 ms and TE was 1.62 ms.

Image post-processing
The images were post-processed using custom algorithms developed using Matlab (R2011b, The MathWorks Inc.). For practical utility, and indeed for consistency, it is important to minimise the need for manual input into the process of extraction of area and flow data. Furthermore the artery moves throughout the cardiac cycle, and it is imperative that this movement is tracked so that the appropriate voxel velocity values within the vessel are summed at each time step to derive the flow. This tracking is achieved by non-rigid registration using the Sheffield Image Registration Toolkit (ShIRT) (Barber et al., 2007), (Lamata et al., 2011). The operator manually identifies the region of the vessel on one image (labelled the fixed image), and the registration tracks automatically the movement and deformation of this region over the cardiac cycle.
The resulting binary images from the segmentation served both to measure the cross-sectional area of the vessel (A(t)) at each time point and as a moving mask for the phase images to identify the appropriate pixels over which to sum velocity data to obtain flow (Q(t)).

Radius as a pressure surrogate
The proposed models require as input the p(t) and Q(t) waveforms. Area is measured non-invasively, and area (or equivalent radius) change is used as a surrogate for pressure. (Greenfield and Griggs, 1963) reported the relationship between the pressure and diameter of the human MPA. The pressure and MPA diameter were measured simultaneously in 11 patients undergoing open heart surgery, 3 with PH. An R 2 value of 0.98 was found between the radius and pressure changes over time. Similarities in the dynamic pressure and radius changes in a large group of children that included PH cases was also shown by (Jarmakaniet al., 1971).
A pressure-area relationship is given by the equilibrium condition at the wall of the vessel (Milisic and Quarteroni, 2004), where the vessel is considered to be formed from a series of rings, tethered axially and in static equilibrium at each time step, with the dynamic and viscoelastic effects negligible.
pðtÞffi Eh where: E is Young's modulus, h is vessel wall thickness, is Poisson's Ratio, R 0 is diastolic radius, and δR is the change in lumen radius at each instant in time.
Fig . 3 illustrates the diastolic pressures, (P d ) and Peterson's elasticity modulus (E p ), (Eq. 2) (Liao et al., 1999), as functions of the area change, using the data set published by Greenfield and Griggs (Greenfield and Griggs, 1963), which has been extended to include our own invasive catheter measurements. The first data set contains valuable invasive measurements in patients investigated for other purpose than PH, which combined with our own data set offers an E p and P d in our process of developing a non-invasive PH tool. A power law, with R 2 ¼0.77, best described the E p -area change relationship, whereas the relationship between P d and area change was best described by an exponential law, with R 2 ¼ 0.40. The E p values of our data set on the graph from Fig. 3.a were computed using the same formula as Greenfield and Griggs in their original paper. Eq.3 (derived from Eqs.1 and 2), together with the relationships given by the curve fitting in Fig. 3 were used to compute p(t) from the measured radius changes in each subject in our study.

0D parameter optimisation
Usually, the challenge in using electrical analogue circuits in bio-medical hemodynamic applications is the determination of the appropriate numerical values of the electrical components. The more complex the circuit is, the more difficult is the parameter optimisation and the more measurements are generally required. For the simple RCR model of Fig. 1a, the governing equations are: where Q 0 ðtÞ and P 0 ðtÞ are the measured flows and pressures at the domain inlet and P c is the pressure drop on the capacitor. The analytical solution of the R c CR d circuit can be written in the frequency domain (in complex representation) as: where ω-angular frequency.
The numerical values of the resistors and capacitors are determined by an optimisation process, minimising the residual between measured and analytical Fig. 3. The Peterson's elasticity modulus,(E p ) (a) and diastolic pressure, P d (b) as a function of area change in a mixed data set: 11 cases using the values reported by Greenfield and Griggs, 1963 (squares) and 27 subjects from our patient cohort (triangles). The relationship between the E p and area change -plotted on a log-log axes is best described by a power law (R 2 ¼ 0.77). The relationship between the P d and area change -plotted on a log-lin axes is best described by an exponential law, with R 2 ¼0.4. For each subject of our study, a P d and an E p value were derived from the two above described relationships and used in the non-invasive estimation of pressure waveform (Eq. 3).

Fig. 2.
Typical Phase Contrast (PC) magnitude and balanced steady state free precession (bSSFP) images of the MPA, during the peak systole and late diastole. The bSSFP sequence has been chosen to be used for the area segmentation during the entire cardiac cycle, due to better blood-vessel delineation than the PC sequence, noticed especially during the diastole. solutions of P 0 (t) for the measured Q 0 (t). The cost function chosen for the optimisation process represents a linear combination of the root mean square error over the cardiac cycle and the error in the peak systolic and minimum diastolic pressures. This is achieved using components from Matlab's Optimisation Toolbox. The results of gradient-based optimisation algorithms can be very dependent on the starting values supplied. Our approach was to run a genetic algorithm optimisation process over a wide interval, representing a range of starting values, and to use the solution from this process as a starting point for a (faster and more efficient) gradient-based method to find the optimum solution in the vicinity of this point. In practice, because the total resistance must equate to the ratio between mean pressure and mean flow (Eq. 7), and this constraint is imposed with the implementation, this is a two-term optimisation process. Although other optimisation approaches are possible, the method adopted in this paper is designed to be robust, requiring no user intervention, and fast.
The hypothesis is that the computed values of R d , R c and C might be used as a differentiation criteria in the characterisation of the proportion of distal and proximal pulmonary vascular resistance in PH.

Wave reflection quantification
Pulse wave propagation in an elastic tube can be described in the frequency domain (Nichols, 2005), or in the time domain using wave intensity analysis (WI) (Parker, 2009). It is noted (Hughes and Parker, 2009) that the frequency domain decomposition method adopted in this manuscript will produce identical results, in terms of the contributions of forward and reflected waves, to an equivalent time domain decomposition based on wave intensity, such as that described by (Hollander et al., 2001;Dwyer et al., 2012).
For application to the current problem, Eqs. (8) and (9) are re-cast into cylindrical co-ordinates and then integrated over a circular cross-section to provide a relationship between pressure (which can be shown to vary minimally in the radial direction and thus can be assumed to be approximately constant over a cross-section) and flow (the integral over the cross-section of the axial velocity). An assumption is made about the shape of the velocity profile over the cross-section to perform the integration. Eq. (10) shows the reduced form when the velocity profile is taken to be pseudo-static (Poiseuille distribution). It is noted that in practice, for the large vessels such as the pulmonary arteries studied in this paper, the viscous terms are in fact small, and dropping them yields the well-known inviscid one-dimensional wave equation.
where c-wave velocity. Our hypothesis in the application of this model is that the healthy system is likely to feature lower reflections (more efficient power transfer) than the diseased system. (Hollander et al., 2001), (Dwyer et al., 2012) have used the same idea to characterise the disease from invasive measurements, expressing the magnitude of the reflections using the ratio between the backward and forward wave intensities. The intensity of each component of the wave can be determined from Eq. 11.
where, W pw -the power of the n-th harmonic, P w -the amplitude of the pressure wave at the n-th harmonic and Z c -the characteristic impedance of the vessel (Eq. 12).
The total power of the forward and the backward wave (W pf , W pb ) was computed for the 1D model as a sum of the power of the individual harmonics after the decomposition of the pressure wave into forward and backward components and the ratio between the total power of the backward wave and total power of the total wave was used as a differentiation criteria.
3. Results and discussions

0D model
The results from the lumped parameter model (Fig. 4) showed that, on average, the Rd increases with the disease status showing statistically significant difference between all the analysed groups, with the exception of the healthy volunteers versus the NoPH group. The R d average values (0.24 mmHg s/ml -NoPH 0.57 mmHg s/ml-PH(PVR o4WU), 1.2 mmHg s/ml -PH(PVR Z4WU)) were comparable with the one reported by (Lankhaar et al., 2006).
Evaluating the contribution of R d to total resistance we found that this contributes 85 75.9% to R total in the healthy volunteers group, 9374.7% in the NoPH group, 9375.6% for the mild PH group and up to 95 76.3% for the severe PH group. The metric defined as the ratio of resistances showed statistically significant difference between the healthy volunteers group and each of the other 3 groups (Fig. 5 a).
The total resistance was found on average to increase as well with the disease state, but no strong correlation was found between it and the measured PVR values (Fig. 5.b).
The R c parameter showed no statistically significant difference between the groups, and they are on average two times smaller than the ones reported by (Lankhaar et al., 2006).
The compliance, C decreased with the disease state, from 4.19 ml/mmHg in the healthy volunteers group, to 0.78 ml/mmHg in the severe PH group, values comparable with the ones reported by (Reuben 1971), (Saouti et al., 2010), (Lankhaar et al., 2006), (Westerhof et al., 2009).

1D model
Using Fourier analysis, the pressure and flow waves were decomposed into forward and backward components and the wave power was computed for each harmonic. Fig. 6 shows typical examples of pressure wave decomposition in a healthy volunteer and a PH patient. There is a clear increase in the magnitude of the reflected wave (dashed-red) in disease compared with the healthy condition.
The visual results are confirmed by the data analysis (Fig. 7). The 1D model showed on average that the backward wave contributes just 10 70.036% to the total wave power in the Healthy Volunteers group, more than double in the NoPH and mild PH groups, 2470.067% and 27 70.13% respectively, and over 41 70.1% in the severe PH group. Significant difference was found between all the analysed groups, with the exception of between the NoPH and the mild PH groups.
The ratio of the power of the backward wave to the total wave correlates linearly (po 0.001) with the measured PVR, with an R 2 value of 0.34. (Laskey et al., 1993), (Huez et al., 2004), using invasive measurements reported increased values of the reflection coefficient in PH patients when compared with normal subjects. (Castelain et al., 2001), (Nakayama et al., 2001), used the pressure waveforms alone to differentiate the patients with CTPH and IPAH. Although their results showed statistically significant differences between the two groups, the method is based on recognition of features of the waveform that can be missing or very subtle in some cases.
Generally both of our proposed models perform well in differentiation of the groups based on the means differences. The wave power metric was not able to show significant difference between the NoPH and the PH group with PVR o4WU group, although the R d parameter can bring complementary information.
Another non-invasive criterion, based on relative area change (RAC) in the MPA between systole and diastole, reported by (Swift et al., 2012) performs better for the evaluation of the NoPH versus mild PH (Fig. 8). However, the proposed 1D model based method was found to be superior to RAC for this small group of patients with PH, differentiating patients with PVR o4WU from those with PVR Z4WU, which highlights the potential added sensitivity of this technique in the assessment of disease severity.

Radius as a pressure surrogate
The primary purpose of this paper is to show that there are statistically significant differences in computed parameters for different groups classified by the presence and degree of PH. These parameters are derived entirely from a non-invasive imaging process, without catheterisation. The process requires the transformation of measured cross-sectional areas to pressures. It is natural to enquire whether this transformation can be used directly to estimate pressure in the individual. Fig. 9 illustrates the correspondence between measured and computed mean, systolic and diastolic pressures, and it can be seen that, although there are consistent trends, the R 2 values indicate that the correlation is relatively poor. This is not surprising because the PA can continue to enlarge independent of the pressure, and the pulmonary arterial dimension at a particular stage in the disease course does not reliably estimate an instantaneous measurement of pulmonary arterial pressure (Boerrigter et al., 2010). Nevertheless the derived parameters have been shown to have significant differences between patient groups. We suggest that this is because the parameters depend not only on the pressure magnitudes, but also on the variations of the waveforms over the cardiac  5. (a) The ratio of the distal resistance to the total resistance in healthy volunteers and stratified PH. On average, the distal resistance R d contributes 85% to the total resistance in the HV group, showing statically difference when compared with the other three defined groups (b) The clinically measured pulmonary vascular resistance (PVR) and total resistance (R total ) of the 0D model linearly correlates with R 2 ¼0.19, p ¼ 0.021. cycle. Fig. 6 illustrates clearly the different waveforms of a healthy volunteer compared with a patient with PH. Our conclusion is that the methodology might be used to develop a diagnostic test for pulmonary hypertension, but cannot be used reliably to derive mean pressure directly.

Conclusions
Our findings suggest that the ratio of reflected-to-total wave power and the electrical analogue parameter, R d from our proposed circuit models derived from dynamic MRI measurements of flow and pulsatility in the PA could be used as a non-invasive measure for the differentiation of the PH patients from NoPH and healthy volunteers. These metrics seem to have the potential to distinguish also between different degrees of disease severity. Further work shall seek to evaluate these metrics in a large prospective cohort of patients drawn from the ASPIRE PH registry (Hurdman et al., 2013) and improve on the methods for noninvasive estimation of pulmonary arterial pressure. The underpinning 0D and 1D models are very well documented in the literature. In our work we have developed an effective optimisation protocol that minimises operator dependence in the numerical characterisation of the model components from measured data. The primary novelties in this work are in linking together the processing steps rather than in the theory of the underpinning models and, most importantly, in the demonstration in this clinical pilot study that PH can be differentiated between healthy  (Swift et al., 2012). Although RAC measurement shows statistically significant differences of the means in between the healthy volunteers, NoPH and the two PH groups, no significance is shown in between the diseased groups, as the 0D and 1D models display.
volunteers and stratified PH using simple computational models of physiology entirely from non-invasive data.  Fig. 9. Correlation between the invasively measured mean (mPAP), systolic and diastolic pulmonary arterial pressures and the non-invasively estimated values.