Computational assessment of hemodynamics-based diagnostic tools using a database of virtual subjects: Application to three case studies

Many physiological indexes and algorithms based on pulse wave analysis have been suggested in order to better assess cardiovascular function. Because these tools are often computed from in-vivo hemodynamic measurements, their validation is time-consuming, challenging, and biased by measurement errors. Recently, a new methodology has been suggested to assess theoretically these computed tools: a database of virtual subjects generated using numerical 1D-0D modeling of arterial hemodynamics. The generated set of simulations encloses a wide selection of healthy cases that could be encountered in a clinical study. We applied this new methodology to three different case studies that demonstrate the potential of our new tool, and illustrated each of them with a clinically relevant example: (i) we assessed the accuracy of indexes estimating pulse wave velocity; (ii) we validated and refined an algorithm that computes central blood pressure; and (iii) we investigated theoretical mechanisms behind the augmentation index. Our database of virtual subjects is a new tool to assist the clinician: it provides insight into the physical mechanisms underlying the correlations observed in clinical practice.


Introduction
Cardiovascular disease is a major cause of morbidity and mortality worldwide (Gersh et al., 2010). Early diagnosis in clinical routine allows the clinician to prevent development of the disease and to offer adequate treatment to the patient. This can be achieved through derivation of indexes describing the function of the heart and large vessels, such as the systolic/diastolic and pulse pressures, ankle-brachial index, and augmentation index (AIx) of central blood pressure. In recent years, emphasis has also been placed on the importance of arterial stiffness in the development of cardiovascular disease (Laurent et al., 2006). Indexes describing the latter include the pulse wave velocity (PWV), arterial distensibility, and stiffness index. Diagnosis can also be achieved by applying specific algorithms to clinical data. For example, the transfer function is a commonly used algorithm that allows the estimation of central pressure from a peripheral pressure measurement (Chen et al., 1997).
Deriving and validating these diagnostic tools on cohorts of patients is a challenging task. In addition to being expensive and time-consuming, in-vivo measurements are subjected to experimental errors, they are not always available at all sites of interest, and fluctuate with the condition of the patient. Most importantly, the estimated index can rarely be compared to its "true" value, since the latter is unknown.
In this study, we propose to use a new methodology based on numerical modeling to assess the accuracy of these diagnostic indexes and algorithms. Computational modeling has proven to be an efficient tool to represent arterial blood flow. In particular, onedimensional (1D) models offer a good balance between accuracy and ease of computation. Several studies have shown their ability to predict arterial pressure and flow waveforms under healthy (Reymond et al., 2011;Mynard and Smolich, 2015) and pathological conditions (Steele et al., 2003;Willemet et al., 2013). Most of 1D modeling studies focus on patient-specific applications, trying to reproduce a single set of parameters and a unique hemodynamic state (Reymond et al., 2011;Huberts et al., 2013;Willemet et al., 2013;Alastruey et al., 2016). This work uses the "populationspecific" approach presented in Willemet et al. (2015). Using 1D modeling and variations in cardiac and arterial parameters within healthy ranges, we generated a database of 3325 healthy virtual subjects presenting a diversity of hemodynamic, structural and anatomical and structural properties, as well as pressure, flow and area waves at the larger arteries; therefore allowing the computation of the exact value of the diagnostic tool. The main assumptions of our 1D formulation are: (i) thin elastic incompressible arterial wall, (ii) 55-artery network, (iii) outflow from the left ventricle prescribed as inflow boundary condition at the aortic root, and (iv) RCR Windkessel models as outlet boundary conditions. We refer the reader to Willemet et al. (2015) for further details on the methodology followed to generate the database.
In this work, we used this new database approach to study the following three clinical problems: assessment of arterial stiffness through PWV (Section 2), validation of an algorithm for estimating central pressure (Section 3), and investigation of theoretical mechanisms behind AIx (Section 4). These examples illustrate the potential of our database of virtual subjects to assess computationally the efficiency of indexes and diagnostic tools based on pulse wave propagation, as highlighted at the end of each section. We then discuss limitations (Section 5), summarize our main results (Section 6) and provide some perspectives for future applications (Section 7).

Assess the accuracy of clinical indexes
Case study: Estimation of pulse wave velocity We investigated the accuracy of three different methodologies proposed to estimate PWV in clinical practice: (i) the foot-to-foot, (ii) loop, and (iii) sum-of-squares methods. We compared them with the theoretical PWV based on the Moens-Korteweg relation (Willemet et al., 2015).

Foot-to-foot PWV
Foot-to-foot PWV (PWV ff ) is the recommended method in clinical practice to estimate arterial stiffness (Laurent et al., 2006). By computing the time delay between the feet of pressure or flow signals at two arterial sites, PWV ff provides an average stiffness along an arterial segment. While most clinical devices compute PWV ff using pressure signals acquired with tonometry or mechanotransducers (Laurent et al., 2006), new methods compute PWV ff using flow waveforms acquired during magnetic resonance imaging (Ibrahim et al., 2010;Wentland et al., 2014). We aimed to assess the accuracy of these measured PWV ff compared to the theoretical value within the wide variation of cases offered by the database.
2.1.1. Method. Foot-to-foot PWV is computed as ΔL=Δt, with ΔL the distance travelled by the pulse wave, calculated as the difference between lengths of wave propagation from the heart; and Δt the transit time between the feet of waveforms, computed using the intersection algorithm described in Gaddum et al. (2013). In this study, we limited our analysis to the PWV ff along the aorta, i.e. from the aortic root to the iliac bifurcation. The corresponding theoretical PWV was computed using Eq. (2) from Willemet et al. (2015).
2.1.2. Result. Along the aorta, both pressure-and flow-based indexes reproduced the theoretical values relatively well, introducing a slight bias (underestimation) of À 12% on average in stiffer arteries (Fig. 1). While flow-based PWV ff tended to strictly under-estimate the theoretical PWV, pressure-based PWV ff both under-and over-estimated the theoretical PWV, with a standard deviation twice larger than the flow-based one.
2.1.3. Discussion. As detailed in our initial analysis (Willemet et al., 2015), this bias in PWV ff is caused by reflected waves present during late diastole; these compromise the foot-to-foot algorithm assumption of a reflection-free period in early systole. These reflected expansion waves are due to negative reflection coefficients at the aorto-iliac bifurcation, which, according to Greenwald et al., 1990, are more likely to occur in older individuals.
Thanks to the virtual database, we could confirm that the aortic PWV can be estimated with good accuracy using the foot-to-foot method, either with pressure or flow signals.

PU-, QA-and lnDU-loop PWV
The loop methods allow the estimation of the localized PWV from simultaneous measurements of two hemodynamic variables, by assuming that there are no reflected waves during early systole. Recent studies have shown that loop PWV indexes deviate from the theoretical PWV value quite significantly (Alastruey, 2011;Borlotti et al., 2014;Segers et al., 2014). We tested these indexes on all virtual subjects. Fig. 1. Comparison of theoretical and foot-to-foot PWV, PWV ff , along the aorta. Both foot-to-foot PWV (computed with pressure (left) or flow (right) waves) are reasonable estimates of the theoretical PWV along the aorta, i.e. from the aortic root to the iliac bifurcation. The larger the stiffness of the elastic arteries (c el ), the larger the dispersion or deviation of the foot-to-foot PWV estimate. The dashed line indicates identity. The left plot presents the same results as in Fig. 6 from Willemet et al. (2015), and is reproduced here for ease of comparison.

Method.
Loop PWV indexes were computed over the linear part of the loop during early systole (Khir et al., 2001;Rabben et al., 2004;Feng and Khir, 2010): with ρ the blood density, P the blood pressure, U the blood velocity, Q the blood flow, A the luminal cross-sectional area and lnD the natural logarithm of the arterial diameter. As suggested by Borlotti et al. (2014), in order to quantify deviations of the loop PWV from the theoretical PWV, we calculated the reflection coefficient at an arterial bifurcation: with Y ¼ A=ðρPWVÞ the characteristic admittance of the parent vessel (Y p ) and daughter vessels (Y d1 ,Y d2 ).

Result.
At the ascending aorta, the PU-loop PWV underestimated the theoretical PWV by 7% while the lnDU-loop method over-estimated it by 10% on average (Fig. 2a,top). PWV loop estimates tended to invert along the aortic length towards the iliac bifurcation ( Fig. 2b and c, top). In peripheral arteries, opposite deviations of the loop methods were observed: PWV PU overestimated the theoretical PWV, while PWV lnDU under-estimated it ( Fig. 2d and e, top). In the carotid artery, the PU and lnDU-loop errors increased up to þ59% and À 36% respectively, while in the iliac artery, these errors were about 78%. In all arteries, PWV QA was very close to PWV lnDU with a slight over-estimation of less than 0.4 m/s. We observed positive reflection coefficients R t at the bifurcation downstream to the left common carotid artery, while negative R t were observed for bifurcations along the aorta, leg and arm arteries. The aorto-iliac bifurcation presented reflection coefficients between À 0.3 and þ0.3.

Discussion.
In Alastruey (2011), the efficiency of local PWV methods was analyzed for a single virtual subject. Similar results were observed, though our results highlight the significant deviation of loop methods along the aorta for a wide variety of cases.
As suggested previously (Alastruey, 2011;Segers et al., 2014), using both lnDU (or equivalently QA) and PU-loop results can provide a more accurate PWV estimate. For example, taking the mean between lnDU and PU-loop PWV reduced the relative error to less than 4% for all arteries, at the exception of the carotid and aortoÀiliac bifurcation where the error decreased to around þ10%. Nevertheless, if pressure and area waveforms are available, one might simply compute the Bramwell-Hill PWV (Nichols and O'Rourke, 2005). Borlotti et al. (2014) showed that loop methods were affected by positive and negative reflection sites: close to positive reflection sites, the PU-loop PWV over-estimated the exact wave speed while the lnDU-loop PWV under-estimated it, and vice-versa for negative reflection sites. In the arterial network, bifurcations (characterized by R t ) are responsible for part of the reflections observed. The above theory is verified for the carotid artery, with positive R t : PWV PU over-estimated the theoretical PWV while PWV lnDU underestimated it (Fig. 2d, top). Results along the aorta, with negative R t , were also in agreement with Borlotti's theory ( Fig. 2a and b, top). The large variation of R t at the aorto-iliac bifurcation might explain the inversion of the deviation of the loop methods towards the distal aorta (Fig. 2c, top). However, some peripheral arteries (e.g. iliac and brachial arteries) did not verify the above theory if arterial bifurcations alone are considered: R t are negative, while the PU-loop PWV over-estimated (and the lnDU-loop PWV underestimated) the theoretical wave speed (Fig. 2e, top). Additional sources of reflections, such as anatomical tapering and secondary re-reflections at bifurcations are not taken into account by R t and might play an important role. The effect of positive reflections due to tapering needs to be verified in further studies.
Our results suggest that the loop methods should be used with care at arterial locations close to reflection sites.

Sum-of-squares PWV (PWV Σ 2 )
The sum-of-squares method was introduced to estimate locally the arterial stiffness by minimizing the net wave energy over a cardiac cycle, in short vessels where foot-to-foot methods cannot be applied (Davies et al., 2006). Initially introduced for the coronary arteries, its use in other arteries has been controversial since PWV Σ 2 is computed over a period when waves are not unidirectional, and the resulting PWV is influenced by reflected waves (Alastruey, 2011;Kolyva et al., 2008).
2.3.1. Method. The sum-of-squares method computes PWV using simultaneous pressure P and velocity U measurements as with dP and dU the changes in pressure and velocity, respectively, across the wave fronts, and the sum extending during one cardiac cycle. Fig. 2, bottom compares the sum-of-squares PWV to the theoretical PWV at five arterial locations: the ascending aorta, thoracic aorta, aorto-iliac bifurcation, carotid and iliac arteries. At the peripheral arteries, PWV Σ 2 over-estimated the theoretical PWV by up to þ140% (carotid) and þ34% (iliac) on average ( Fig. 2d  and e, bottom). Along the aortic arch, there was also an overestimation of the theoretical PWV by 50% (ascending aorta, Fig. 2a, bottom) and 20% (descending aorta, results not shown here) on average. The agreement at the thoracic aorta was better on average with an under-estimation of 2%, even though individual models deviated from the theoretical value by up to 60% (Fig. 2b, bottom).

Result.
Similarly to the loop PWV indexes, on average, a large deviation (þ 58%) was observed at the aorto-iliac bifurcation (Fig. 2c, bottom). The over-estimation of PWV Σ 2 was also observed at the other arteries of the database.
2.3.3. Discussion. A detailed analysis based on the reflection coefficient (not displayed here) confirmed that the bias of PWV Σ 2 was related to the amount of reflected waves in the arterial system. In view of these results, the sum-of-squares method should be used only with great care in clinical practice, as this method does not always estimate the PWV accurately and is strongly influenced by wave reflections.

Summary
The virtual database enables us to: -Test the accuracy of an index based on pulse wave propagation since its "true" theoretical value can be computed exactly.

Validate your algorithm without experimental error
Case study: Non-invasive estimation of central blood pressure Central blood pressure provides major information on the state of the heart and cardiovascular system. While there exists a number of non-invasive methods for estimating central systolic pressure (Chen et al., 1997;Guilcher et al., 2011), they are of limited accuracy and difficult to apply during cardiac magnetic resonance imaging. A new algorithm, which estimates the whole pressure waveform in the ascending aorta from non-invasive measurements, has recently been suggested (Vennin et al., 2015). This algorithm divides the target pressure waveform into an early systolic upstroke determined by the water hammer equation, a late systolic portion described by a secondorder polynomial constrained by conditions of continuity and conservation of mean arterial pressure, and a diastolic decay equal to that measured in peripheral arteries. The algorithm uses the local aortic flow and stiffness, the mean and diastolic blood pressure, and the diastolic pressure decay to generate the central pressure waveform (Fig. 3). Thanks to waves being available at central and peripheral locations, our virtual database allows us to assess the efficiency of the new algorithm at reconstructing the aortic pressure waveform.

Method
In the initial proof-of-concept study (Vennin et al., 2015), the aortic stiffness was evaluated using the sum-of-squares PWV at the ascending aorta (Eq. (3)). Because PWV Σ 2 requires simultaneous measurements acquired invasively, in this study we evaluated alternative non-invasive methods to derive aortic stiffness: pressure-based carotid-femoral PWV (cfPWV ff ) and brachial-ankle PWV (baPWV ff ) acquired with tonometry, flow-based ascending to descending aorta PWV (adPWV ff ) and descending to thoracic aorta PWV (dtPWV ff ), lnDU-and QA-loop PWV at the ascending aorta, Fig. 3. Schematic of the methodology of a recently proposed algorithm for estimating central blood pressure. Inputs include the brachial pressure waveform (P), the aortic velocity waveform (U) and the aortic stiffness, estimated by PWV. The output of the algorithm (the estimated aortic pressure waveform) is compared to the exact aortic pressure simulated using the 1D model for each virtual patient. The error between the systolic shoulders P 1 of the estimated and reference central pressures is then computed. Fig. 4. Distribution of the error in systolic pressure shoulder for different PWV methods. The error ϵ ¼ P 1 algo À P 1 ref between the reconstructed (P 1 algo ) and the reference (P 1 ref ) systolic pressure shoulders is computed using 7 different PWV methods. Each box indicates the 25th percentile, median, and 75th percentile; whiskers extend to minimum and maximum data points. Outliers are plotted in grey. cfPWV ff : pressure-based carotid-femoral foot-to-foot PWV; baPWV ff : pressure-based brachial-ankle foot-to-foot PWV; adPWV ff : flow-based ascending to descending aorta foot-to-foot PWV; dtPWV ff : flow-based descending to thoracic aorta foot-to-foot PWV; lnDU-loop, QA-loop and sum-of-squares (Σ 2 ) PWV were evaluated at the ascending aorta.
with distension and flow waveforms acquired during MRI (Wentland et al., 2014). In the central pressure algorithm of Vennin et al., (2015), the PWV mainly affects the systolic increase and specifically the estimation of the systolic shoulder P 1 (Fig. 3). We therefore assessed the accuracy of P 1 obtained from each PWV index by computing the error: where P 1 algo is the systolic shoulder pressure estimated by the algorithm and P 1 ref is the exact reference systolic shoulder pressure.

Discussion
In view of the deviation of the sum-of-squares method at the ascending aorta (over-estimation of 50% of the theoretical value, see Fig. 2a, bottom), this analysis highlights that the proposed algorithm introduces a bias in the estimation of the systolic shoulder, which is balanced when using the sum-of-squares method. Note that the database presents only type-A waveforms, with a systolic shoulder occurring before the systolic peak (Murgo et al., 1980). Furthermore, that shoulder is moderately marked. Therefore, the conclusions drawn here should only apply to that category of cases.
Given the wide collection of data available without experimental errors, the database has been a useful tool to help understand and improve quickly the methodology of a clinically relevant algorithm.

Summary
The virtual database enables us to: -Test an algorithm based on pulse wave propagation when accurate in-vivo data is not available. -Avoid inaccuracies from experimental error.

Understand theoretical mechanisms of wave propagation
Case study: Investigating the augmentation index as a measure of wave reflections The AIx is used clinically as a measure of wave reflection in the arterial network (Nichols and O'Rourke, 2005). It is assumed that the systolic shoulder visible on the pressure wave results from the arrival of the reflected wave interacting with the forward flow from the heart. Recently however, serious limitations about the validity of that concept have been suggested (Hughes et al., 2013), and it has been proposed that AIx may be indicative of the compliance of elastic arteries (Davies et al., 2010) and ventricular function (Fok et al., 2014). We assessed the relation between AIx and reflections using the database.

Method
AIx was extracted from the pressure wave shape and was computed as the ratio of the augmented pressure (the pressure difference between the systolic peak and shoulder) to the pulse pressure. The shoulder point was detected by the zero-crossing of the fourth derivative of the signal, using the algorithm detailed in Segers et al. (2005;2007), with a Savitzky-Golay differentiation filter of order 6 over 35 points. Reflections from the arterial network were assessed by decomposing the pressure wave into its forward (P f ) and backward (P b ) components, and by computing the ratio of backward to forward pressure peaks (P b /P f ), as well as the timing of the backward peak. We also compared the AIx to the average reflection coefficient at all bifurcations of the network (R t ), with each R t computed using Eq. (2).

Result
Results are presented at the ascending aorta and carotid artery (Fig. 5). At both arteries, AIx showed a weak to moderate linear correlation with the ratio of amplitudes P b /P f (Pearson's correlation coefficient r ¼ 0:47 (aorta) and r ¼ À0:29 (carotid), Fig. 5a). The timing of the shoulder point also only correlated weakly with the timing of the peak of the backward wave (r ¼ À0:34 (aorta) and r ¼ 0:27 (carotid), Fig. 5b). However, at the carotid artery, we observed an important correlation (r ¼ 0:82) between AIx and R t , with AIx increasing as R t tends to zero (Fig. 5c). This correlation was not as marked in the ascending aorta (r ¼ 0:28).

Discussion
According to our results, AIx is not related to reflections when those are described through the decomposition into forward and backward pressure components. However, some correlation can be observed if one describes reflections by geometrical and structural properties of the arterial network. The difference between observed and theoretical reflections highlights the limited applicability to clinical data of the theory behind AIx. Indeed, as we are considering results from numerical modeling -without any experimental error-, one would expect much stronger correlations between AIx and reflection indexes, if AIx were a robust measure of wave reflection in the arterial network. As stated previously, our current database contains only type-A waveforms, with a moderate systolic shoulder occurring before the systolic peak. As highlighted by Hughes et al. (2013), different associations might be observed if type-C waveforms are studied.
While this analysis does not go into the details of providing a clear explanation of the origin of the augmentation index, it highlights that arterial wave reflections, as assessed in clinical literature, are not the only mechanism influencing the systolic shoulder. Interpreting the augmentation index as a measure of wave reflection only should therefore be made with caution.

Summary
The virtual database enables us to: -Understand physical mechanisms of pulse wave propagation. -Challenge established concepts, and highlight their limitations.

Limitations
The results presented in this study depend on the characteristics of the 1D model of arterial hemodynamics used to generate the database.
The model produced pressure waveforms which are representative of a healthy adult subject older than 30 years of age, rather than a young adult. These are type-A waveforms with an early systolic shoulder (Murgo et al., 1980). In future works, we plan on using an improved description of the wall elasticity, together with a more detailed model of the heart (Mynard and Smolich, 2015) in order to address that limitation and generate type-C waveforms as well.
Our model contains the 55 larger systemic arteries; it includes neither the coronary circulation or hand arteries, nor a detailed cerebral arterial network. Extending the arterial network to include more peripheral vessels is worth considering for the following two reasons. (i) Non-invasive measurements can be acquired at peripheral sites: e.g. volume pulse using photoplethysmography at the finger tip or on forehead arteries. (ii) It may improve the accuracy of signals at our current terminal vessels; e.g. the carotid arteries at the entrance of the complex network of cerebral arteries. Another improvement to the current model relates to the length of the arterial network. Arterial lengths could be varied uniformly to consider the influence of subject height.
Despite these limitations, the database offers a wide range of hemodynamic waveforms under physiological conditions. We believe that it is representative of a subset of the "real population"; this still needs to be demonstrated in a comparison with a cohort of human-derived data in future works.

Summary
We have presented a new methodology to assess theoretically the efficiency of diagnostic tools (indexes and algorithms) based on pulse wave propagation. It consists of a database of virtual arterial waveforms generated by a wide range of physiological, structural and geometrical properties. We have applied this new approach to three different tools used in the clinic to study hypertension.
Thanks to the availability of the "true" theoretical value of arterial stiffness for each virtual subject, we have studied the accuracy of pulse wave velocity estimates. According to our results, PWV estimates are affected by wave reflections. Nevertheless, the central foot-to-foot PWV seems to be the most accurate -and easier to use-in clinical practice. By running a new algorithm for central pressure estimation on all virtual subjects without any experimental error, we identified some recurrent bias in the proposed algorithm due to PWV. Finally, by applying wave separation analysis on central waveforms, we challenged the established concept that the augmentation index is an estimate of arterial wave reflections.
These case studies have highlighted the potential of the database approach to assess theoretically the efficiency of indexes and diagnostic tools based on pulse wave propagation. Furthermore, the database has allowed us to identify important mechanisms of wave propagation relevant to pulse wave velocity and augmentation index.

Perspectives
We believe that the database of virtual subjects can become a new useful tool for the clinician, providing insights into mechanisms important for the design of large cohort studies, and advice on the most relevant in-vivo data acquisition protocol. Indeed, if an index or algorithm fails to be validated within the theoretical framework of the virtual database, it would then be even less likely to perform well in a clinical setting.
Our database is also a very useful tool for biomedical students as it offers hemodynamic waveforms at various locations of the arterial tree on which one can apply post-processing algorithms and theoretically test hemodynamic concepts.
This study has focused on three particular clinical case studies, but future works could assess other clinical indexes (e.g. stiffness or ankle-brachial indexes), and evaluate new diagnostic tools such as patient-specific transfer functions (Hametner et al., 2015).
Future extensions of the database could include diseased conditions such as peripheral arterial disease or heart pathologies.
The complete database is available for download on the project website: www.haemod.uk/virtual-database.html, while the numerical code used to run simulations and generate virtual cases can be found at www.haemod.uk/nektar.html.