Optimization of topological complexity for one-dimensional arterial blood flow models.

As computational models of the cardiovascular system are applied in modern personalized medicine, maximizing certainty of model input becomes crucial. A model with a high number of arterial segments results in a more realistic description of the system, but also requires a high number of parameters with associated uncertainties. In this paper, we present a method to optimize/reduce the number of arterial segments included in one-dimensional blood flow models, while preserving key features of flow and pressure waveforms. We quantify the preservation of key flow features for the optimal network with respect to the baseline networks (a 96-artery and a patient-specific coronary network) by various metrics and quantities like average relative error, pulse pressure and augmentation pressure. Furthermore, various physiological and pathological states are considered. For the aortic root and larger systemic artery pressure waveforms a network with minimal description of lower and upper limb arteries and no cerebral arteries, sufficiently captures important features such as pressure augmentation and pulse pressure. Discrepancies in carotid and middle cerebral artery flow waveforms that are introduced by describing the arterial system in a minimalistic manner are small compared with errors related to uncertainties in blood flow measurements obtained by ultrasound.

Introduction

Computational models of the cardiovascular system are commonly separated into three-dimensional (3D), one-dimensional (1D) and lumped models (0D). One of the first attempts to model pressure and flow waveforms was through the classical 0D Windkessel (WK) model [1]. A noteworthy extension to this was presented in [2] where a resistance element representing the characteristic impedance was added, and many variations and extensions have been proposed [3]. The most important drawback of the family of 0D models is inherent in the assumption of infinite wave velocity and that spatially distributed parameters are modelled as single point parameters.

Through the years distributed models with various degrees of detail have been suggested. In [4,5], the systemic circulation was modelled as two asymmetric parallel branches, one supplying the head and upper limbs, and one supplying the rest of the body. In [6], a model consisting of the 33 largest systemic arteries was tested using an in vitro experiment. In [7], the arterial network was expanded to include 55 arterial 1D segments. In [8], a complete description of the systemic arterial tree containing the largest arteries of the head and upper and lower body was validated using in vivo measurements. The study also includes a detailed overview of 1D models up until 2009, highlighting their variation in detail and complexity. More recently, in [9], a model accounting for pulse wave propagation in all regions of the circulation including approximately 400 arteries and 350 veins was presented. Yet others have modelled the arterial system in a very high level of detail including more than 2000 arteries [10,11].

We have come a long way in creating realistic and detailed descriptions of the entire arterial tree and circulatory system. However, given the near endless number of small arteries and capillaries in the human body, the network has to be truncated at a certain level. Since reliable measurements of flow or pressure at all terminal sites are practically impossible to obtain, outflow boundary conditions are commonly set through simpler models representing the peripheral circulation. Indeed the above-mentioned family of 0D WK models have been the preferred choice for providing boundary conditions at terminal branches.

There is little consensus in the scientific community on the level of detail of the computational domain. Furthermore, few studies have focused on the errors and limitations associated with truncating the arterial network at given sites. In [8], they state that a detailed description of the cerebral circulation is required in order to attain accurate and physiological flow predictions in the common carotid artery. In [12], they found that the arterial tree could be truncated after the first generation of bifurcations without significantly altering pressure and flow waveforms, if matched three-element WK outflow models were used. In [13], they presented a method for lumping 1D arterial segments into three-element WK models and applied their method on a network of 55 arteries (excluding the circle of Willis).

Here, we present a sound mathematical framework that enables us to find the necessary arteries to include for a given clinical application. The framework involves finding the model with the fewest number of arteries that is still able to produce pressure and flow waveforms below a certain error threshold compared with a corresponding detailed (baseline) model (figure 1). This approach reduces the number of uncertain input parameters, while still assuring that the simplifications do not limit the model predictions. We illustrate the framework for different clinically relevant quantities of interest: central aortic and larger systemic artery pressure waveforms, common carotid and middle cerebral artery flow waveforms and coronary pressure waveforms. We note that our framework is intended to be used in an early stage as a tool for model selection that aims at minimizing total uncertainty.

Figure 1.

Two baseline models were used in this work: a model containing 96 arterial segments in which parameters and topology were adapted from [9] (a), and a patient-specific coronary network (b). The arrow indicates the location of invasive pressure measurements, and the section coloured in red is a significant stenosis.

Material and methods

Framework for balancing topological complexity with model error

Here, we present a framework for reducing the number of vessel segments still assuring wanted features of pressure and/or flow to be within acceptable agreement with the corresponding full model:

Define a baseline model.

Locate the quantity of interest appropriate for the problem (e.g. aortic pressure and/or carotid flow).

Define a threshold for pressure and/or flow (e.g. RMS-error, pulse or mean pressure).

Create reduced models by applying the methods described in §2.4.1 or §2.4.2, and solve the 1D networks.

Find the network with the fewest number of arteries subject to the constraint of the threshold.

Arterial baseline models

We applied our new methodology on two arterial models, both illustrated in figure 1.

Systemic arterial network

The first baseline model considered includes 96 of the largest systemic arteries, in which parameters and geometry were adapted from Mynard et al. [9]. They compared model-derived pressure/flow waveforms with published in vivo waveforms from healthy adults, validating the model’s capability of providing realistic waveforms throughout the arterial tree.

Coronary network

The second baseline model considered in this work was based on a series of invasive and non-invasive measurements of a patient (sex: female, age: 58, height: 162 cm, weight: 78 kg) with positive findings of stable coronary artery disease after clinical inspection and coronary computed tomography angiography (CCTA) examination. The data were collected as part of an ongoing clinical trial at St Olavs Hospital, Trondheim, Norway [14]. Cardiac output (CO) was measured by transthoracic Doppler echocardiography using a GE Vivid E95 scanner (GE Vingmed Ultrasound, Horten, Norway). The patient was further referred to invasive angiography, and a Verrata Plus (Philips Volcano, San Diego, USA) pressure wire was used to obtain pressure tracings at the coronary ostium and distal of an epicardial stenosis. Proximal, Pp and distal, Pd pressure tracings are shown in figure 8. The last 30% of the cardiac cycle is highlighted and was used to compute the instantaneous wave-free ratio (iFR), which is a drug-free index of the significance of the stenosis [15]. Measurement of fractional flow reserve (FFR) [15], obtained during drug-induced hyperaemia (maximum coronary flow) was also available. The coronary geometry was segmented using the open-source software ITK-SNAP [16], the surface was then meshed using the open-source library Vascular Modeling ToolKit [17]. 1D domains were extracted from the 3D volume mesh by computing equivalent axisymetric cross-sectional areas along centrelines. Stenotic regions were automatically detected using a Gaussian filter-based approach [18].

Figure 8.

Proximal Pp and distal Pd pressure waveforms obtained during invasive angiography, together with a computed left ventricle pressure waveform in grey. The last 30% of the cardiac cycle is highlighted in red and was used to compute the instantaneous wave-free ratio, iFR, a drug-free index of the significance of the stenosis.

Numerical formulation

One-dimensional flow solver

The solutions of pressure and flow waveforms presented here were obtained using the 1D flow solver STARFiSh [19]. The hyperbolic partial differential equations for blood flow in compliant vessels are written in terms of pressure and flow variables (P, Q): and and solved using the explicit MacCormack scheme [20]. Here, t is the time, x is the axial coordinate, f is the frictional term and is given by −2(ζ + 2)μπU, where ρ is the density (1060 kg m−3), μ is the viscosity of blood (3.5 mPa s), A is the cross-sectional area and U is the cross-sectional averaged velocity. The following velocity profile was prescribed: where r(x, t) is the lumen radius, ξ is the radial coordinate and ζ = 9 is the polynomial order. At arterial connections compatibility of propagating characteristic variables were enforced [7] in addition to conservation of mass and a coupling equation for the pressure, i.e.: and where N is the number of vessels in the connection, and ΔP is an additional pressure loss which was set equal to zero for normal connections. At arterial stenoses, the flow regime is 3D and the 1D assumptions no longer hold. Stenotic regions were thus removed and treated as junctions with N = 2, however, now with an additional experimental-based pressure loss term given by Liang et al. [21]: where the viscous, Kvisc and expansion, Kexp coefficients were calculated based on geometrical features, as described in [21].

The pressure–area relation assumes thin-walled elastic vessels and can be derived from Laplace’s Law: where Pdia is the diastolic pressure with corresponding cross-sectional area Ad, E is the elastic modulus and h is the thickness of the vessel wall. The stiffness parameters E h are related to the pulse wave velocity c and have been obtained using the relation [22]: where rd is the radius at diastolic pressure, and the values for k1, k2 and k3 were set to 3 × 106 g s−2 cm−1, −9 cm−1 and 33.7 × 104 g s−2 cm−1 for systemic arteries and 20 × 106 g s−2 cm−1, −22.5 cm−1 and 86.5 × 104 g s−2 cm−1 for coronary arteries, respectively [9].

Boundary conditions

For the 96-artery model, inflow boundary conditions (prescribed flow rate Q) and outflow boundary conditions (three-element Windkessel models, WK3) and all other parameters were adapted from Mynard & Smolich [9]. For the coronary network, the proximal pressure tracing was prescribed at the aortic root. In contrast to systemic arteries, coronary arteries experience increased impedance during systole due to the contraction and increased pressure in the left ventricle. To account for this effect, a lumped parameter WK model WKcor was used at coronary outlets [23]. A schematic of the model is shown in figure 9 in appendix A.1 and the a priori computed left ventricle pressure waveform is shown in figure 8. The left ventricle pressure waveform was obtained by coupling a varying elastance (VE) heart model to a WK3 model [24], and further by parameter optimization to minimize the discrepancy between Pp and Pao, where Pao is the aortic pressure resulting from the VE-WK3 model. The total arterial resistance, Rtot was estimated from CO, mean arterial pressure, and outflow WK pressure, Pout,WK (5 mmHg) according to Ohm’s Law. Total arterial compliance was estimated from the VE-WK3 model. About 4.5% of CO was assumed to supply coronary arteries and used to estimate total coronary resistance and compliance, and was further distributed among coronary outlets according to Murray’s Law [25]. Simulation of a hyperaemic state is necessary for FFR calculations. Hyperaemia was modelled by reducing the resting resistance of the coronary outlets by a factor α. The value of α was based on the work of Uren et al. [26] who studied myocardial blood flow and resistance in relation to the severity of coronary stenosis, and was set to 3 for ‘healthy’ outlets, and to 1.25 for outlets distal of the coronary stenosis. For details see appendix A.2.1.

Figure 9.

Schematic of the two lumped parameter models used in this work, WK3 model (a) and WKcor (b). R1, R2, Rp, R and Rd are resistance parameters, C, Ca and Cm are compliance parameters and PLV and Pout,WK are left ventricle and outflow windkessel pressures, respectively.

Network reduction

Network reduction involves lumping distributed 1D segments into 0D parameter models, specifically WK models, intended to represent the same physical problem. Each WK model represents all arteries situated distal of the point of interest with resistance elements and capacitors in series and parallel, as visualized in figure 2.

Figure 2.

An example of a reduced network that has been obtained from the baseline model in figure 1 by lumping 1D model segments into WK3 models. Lumped 1D model segments are shown in colour.

Method 1, algebraic estimation of lumped parameters

Here, we present a method for network reduction which was adapted from Epstein et al. [13]. The method was described and applied on a baseline network only including bifurcations. In this work, we have used a different way of estimating the lumped resistance and compliance that can also be applied on networks containing loops and anastomosis. We have also expanded the procedure to account for arterial stenoses.

Estimation of lumped resistance

The linearized version of equations (2.1a) and (2.1b) can be written in terms of the steady-state variables , and : and where l is the length of the segment, and the subscripts ‘in’ and ‘out’ denote variables at the inlet and outlet of the segment, respectively. Equations (2.7a) and (2.7b) may then be combined with equations (2.3a)–(2.3b) and equation (2.5) to form a system of nonlinear algebraic equations. The system was solved iteratively by employing Picard linearization. and is in such an estimate of the time average of P(t) and Q(t), and once solved for, resistance may be estimated anywhere in the network using Ohm’s Law:

Estimation of lumped compliance

We can estimate the compliance (Cv) of a vessel by integrating over the length of the 1D model segment [13]: Furthermore, we estimated the compliance Ct of a terminal vessel (figure 3) coupled with a WK3 with proximal resistance, R1, compliance, C and peripheral resistance, R2 [13]:

Figure 3.

Illustration of an arterial 1D model segment coupled to a WK3 model that may be lumped into an equivalent compliance Ct according to equation (2.10). (Online version in colour.)

Lumped compliance of terminal vessels coupled to WKcor models (see figure 9 in appendix A.2.1.) with compliances Ca and Cm were estimated according to: The total compliance contribution of vessels distal of a point of truncation was then obtained using equation (2.9) for non-terminal vessels and equation (2.10) or equation (2.11) as appropriate for terminal vessels, together with summation rules for compliances/capacitors in series and parallel. See appendix A.2.1. for details.

Lumping vessels distal of a site of truncation

With the lumped resistance, (equation (2.8)) and compliance (equations (2.9)–(2.11)), as defined above we may replace all vessels distal of a point of interest with a WK model. Systemic arteries were replaced by WK3 models in which R1 was set equal to the characteristic impedance, Z: Lumped coronary arteries were replaced by WKcor models and the lumped resistance and compliance were divided among the resistance and compliance parameters of the WKcor model as described in appendix A.2.1.

Method 2, optimization of lumped parameters

Method 1 is based solely on the topology and properties of the baseline model. This means that we can use the method without solving the baseline model. However, the parameters in the WK models that replace the removed vessels are not necessarily the ones that correspond with the least discrepancy between the baseline and reduced networks. This motivates another method which is based on parameter optimization. Since the WK models are lumped models with governing ordinary differential equations (ODEs), we suggest a procedure that treats every truncated site independently. The optimization is thus performed by taking the flow from the 1D solution of the baseline model as given inflow to the WK models, then solving for the unknown pressure. Furthermore, we seek to minimize the error between the pressure obtained by solving the ODE with the corresponding 1D baseline solution. In the following, we explain the procedure for the WK3 model, though it can be easily expanded to other lumped parameter outflow models. Either one, two or all three of R1, C and R2 were allowed to vary to minimize the error. If only one of R1 and R2 was optimized, the total resistance R1 + R2 was found from (Pavg − Pout,WK)/Qavg, where Pavg and Qavg are the time-averaged pressure and flow from the 1D baseline solutions. The method may be summarized in the following steps:

Calculate the flow and pressure waveforms of the 1D baseline model.

Locate the sites where WK3 models will replace distal vessels.

Calculate values of R1 + R2 from Pavg, Qavg, and C using Method 1 (§2.4.1).

Use the flow from the 1D baseline model as given inflow of the WK3 ODE, with parameters R1, C and R2.

Choose parameters to be optimized and use parameters from point 3 otherwise and as initial guess.

Solve the WK3 ODE for the unknown pressure, PWK3.

Find the parameters that minimize the discrepancy between PWK3 and the corresponding pressure waveform from the solution of the 1D baseline model. We used the average relative error, calculated by equation (2.13a) as the measure of discrepancy.

Based on a parameter correlation and identifiability analysis, we chose to optimize on the subset of parameters ([θ1, θ2] = [R1/R2, C]), where R1 + R2 was kept constant. See appendix A.3.2. for details.

Error metrics

The following error metrics were used to compare pressure and flow waveforms obtained from the baseline (B) and reduced (R) models: and where N is the number of time points in a cardiac cycle, i represents a certain time point with corresponding baseline, and reduced, pressure and flow (, ), respectively. ε was normalized by the maximum flow of the baseline model over one cardiac cycle, , to avoid division by numbers close to zero. The maximum (Psys) and minimum pressure (Pdia) was used to calculate the systolic (ε), and diastolic (ε) error, respectively. The pulse pressure, PP is defined as Psys − Pdia. ε is the error in pulse pressure and ε is the error in augmentation pressure, both normalized by the pulse pressure. is the pressure at the inflection point in early systole [27]. εiFR is the difference between predicted iFR from baseline and reduced model.

Application to different physiological and pathological states

The parameters for the baseline 96-artery model were based on data from healthy, young adults [9]. In this part of the study, however, we re-parametrized a series of optimal networks to represent (1) normal ageing, (2) a pathological state of aortic coarctation and (3) states of different heart rate, ejection time and stroke volume. We note that no information from the baseline model was used to re-parametrize the reduced models.

Normal ageing

Normal ageing was simulated by increasing total arterial resistance by a factor of 1.1, and decreasing total arterial compliance by a factor of 2. Arterial stiffening is most marked in the proximal aorta and its major branches—brachiocephalic, carotid, subclavian [28]. The stiffness parameter β for these arterial segments was increased by a factor of 2.5, whereas it was increased by a factor of 1.5 for all other segments. Finally, the compliance of the WK3 models were modified so that the total arterial compliance (sum of WK3 compliance of terminal segments and integrated 1D compliance) was decreased by a factor of 2. The total arterial resistance was modified by increasing the peripheral resistance in all outflow WK3 models. See appendix A.5 for details.

Aortic coarctation

Aortic coarctation was simulated by introducing a 1 cm long, 50% diameter stenosis in the thoracic aorta. This corresponds to segment Id 18 in the electronic supplementary material.

Heart rate, ejection time and stroke volume

Heart rate, ejection time and stroke volume were modified according to the study by Weissler et al. [29]. They studied relationships between left ventricular ejection time, ET, stroke volume, SV and heart rate, HR, in normal individuals. We modified the original aortic inflow curve for the 96-artery model to represent the two extreme cases in terms of HR in their study (HR: 56 bpm, ET: 0.315 s, SV: 106 ml and HR: 120 bpm, ET: 0.2 s, SV: 44 ml). For the latter, total arterial resistance was increased by a factor of 1.67 and compliance halved (effecting the distributed parameters as described for normal ageing), in order to obtain physiological pressure waveforms.

Results

Comparison of Method 1 and Method 2 for network reduction

Figure 4 shows the 96-artery model (black) reduced to a 25-artery model (red). Solution of pressure and flow waveforms at the inlet of the right internal carotid artery, obtained from the baseline model and both methods for network reduction, are also shown. Method 1 overestimated internal carotid pressure in mid systole (ε was 6.2% for Method 1 and 0.2% for Method 2). Furthermore, Method 2 captured the overall shape of pressure and flow waveforms better than Method 1. Average errors, ε between full and reduced models were 1.45% for Method 1 and 0.57% for Method 2. Similarly, ε was 1.47% and 1.16%, respectively. Figure 4 also shows the impedance modulus and angle for the site of interest, calculated in the frequency domain as explained in [30].

Figure 4.

Comparison of Method 1 and Method 2 for network reduction. Baseline model in black and reduced models in red together with pressure and flow waveforms at the inlet of the right internal carotid artery. Impedance modulus and phase angle are also shown.

Framework for optimizing topological complexity

A summary of the quantities of interest, error metrics and values for the network reduction framework applied on the 96-artery model is given in table 1. Here, error metrics are also presented for the cases where parameters were altered to simulate different physiological and pathological states (see §2.6). References to associated figures are also given. In particular, a threshold based on ε + ε at the aorta and brachial artery was used in the top two examples in figure 5. The waveforms for the baseline model and optimal reduced networks are shown in solid lines, and the dashed lines represent the case when the models were altered to represent normal ageing. In the last example, a threshold based on augmentation and pulse pressure was used (). Furthermore, in order to ensure that interaction between different regions in the network and that pressure propagation are correctly captured throughout the larger systemic arteries, a threshold based on pressure waveforms at four locations was used in figure 13 in appendix A.5. Here, the average ε for the aortic root, common carotid, brachial and femoral artery pressure waveforms was required to be less than 0.4%. Additionally, results are shown for ε less than 0.9 and 3.4% for the right common carotid artery and ε less than 0.6 and 1.6% for the middle cerebral artery in figure 6. Method 2 (§2.4.2) was used to reduce the networks in all these cases.

Table 1.

Summary of results from applying the framework outlined in §2.1, on the 96-artery baseline model. For cases where there are more than one quantity of interest, the final error was calculated as the average of the error for the individual quantities. Ref. denotes the reference case, and the threshold used for the optimization is given in brackets. The errors are also shown for states of normal ageing, aortic coarctation (coarc.) and for the two aortic inflow curves as defined in §2.6. All errors are in percentage. The associated figure numbers are referenced below the error, where available.

quantity of interest	no. of arteries	error- metric	ref.	ageing	coarc.	inflow 1	inflow 2
aortic and brachial pressure	29	εP,sys +εP,dia	0.25 (0.3) 5	2.26 5	0.37 —	0.14 —	1.58 —
aortic and brachial pressure	15	εP,sys +εP,dia	0.92 (1.0) 5	1.0 5	0.27 —	1.13 —	1.3 —
aortic pressure	31	εPP +εP,aug	0.68 (0.7) 5	0.97 5	1.91 —	0.41 —	1.27 —
aortic and brachial and carotid and femoral pressure	31	εP,avg	0.33 (0.4) 13	0.47 13	0.24 15	0.41 16	0.85 16
carotid flow	25	εQ,avg	0.87 (0.9) 6	1.53 14	1.09 —	0.81 17	2.97 17
carotid flow	5	εQ,avg	3.35 (3.4) 6	5.2 14	— —	2.93 17	6.3 17
r. middle cerebral flow	38	εQ,avg	0.59 (0.6) 6	1.13 14	0.66 —	0.66 17	2.88 17
r. middle cerebral flow	15	εQ,avg	1.58 (1.6) 6	2.64 14	— —	1.72 17	4.38 17

Figure 5.

Optimal networks and corresponding waveforms (solid lines) obtained from the 96-artery baseline model with pressure at midpoint of ascending aorta (left) and right brachial artery (middle) set as quantities of interest. In the top and middle rows, an averaged (of the two quantities of interest) error threshold of ε + ε less than 0.3% (top row) and 1.0% (middle row) was used. In the last case, an error threshold of ε + ε less than 0.7% was used for the aorta. Dashed waveforms correspond to simulations of normal ageing as described in §2.6.

Figure 13.

Optimal network obtained from the 96-artery baseline model with pressure at midpoint of ascending aorta (top left), right common carotid artery (top middle), right brachial artery (bottom left) and left femoral artery (middle bottom) set as quantities of interest. The average deviation at these four locations was used as threshold, and set to 0.4% average error. Dashed waveforms represent the case when baseline and reduced networks were re-parametrized to simulate normal ageing (see §2.6).

Figure 6.

Optimal networks obtained from the 96-artery baseline model with the flow at the distal end of the right carotid artery set as quantity of interest (left panel) and with 0.9 and 3.4% average error thresholds, respectively. The right panel show results when the flow at the proximal end of the right middle cerebral artery was set as quantity of interest and with 0.6 and 1.6% average error thresholds, respectively.

In the top part of figure 7, εiFR was set to 0.033, which is the standard deviation of repeated iFR measurements, according to the study by Johnson et al. [15]. The results are visualized through the distal pressure waveform, Pd. All side branches except those distal of the measured location can be replaced by lumped WKcor models with no visible effect and with εiFR < 0.000012. If the threshold is increased to 0.04 the network can be reduced to its most simplistic realization, as visualized in the bottom part of the figure. The predicted velocity and the in vivo pressure waveforms are also shown. iFR was measured to 0.40, whereas the predicted value was 0.42 for the baseline network, and 0.42 and 0.38 for the reduced networks, respectively. For FFR, the measured value was 0.52, whereas the predicted value was 0.48 for the baseline network and both of the reduced networks. Method 1 (§2.4.1) was used to reduce the coronary networks.

Figure 7.

Optimal networks obtained from the coronary network with difference in iFR set as quantity of interest, with εiFR < 0.00012 (top) and εiFR < 0.04 (bottom). Results are visualized through the distal pressure waveform. In vivo pressure waveform and predicted velocity waveforms are also shown. The arrow indicates the location of the waveforms, and the red section a significant stenosis.

Discussion

In this study, we have presented a novel approach which optimizes the number of arterial segments for 1D blood flow models. We have illustrated the framework on a 96-artery and a coronary baseline model, and two methods for network reduction have been incorporated: a purely algebraic method (Method 1, §2.4.1) and a novel method based on optimization (Method 2, §2.4.2).

Comparison of methods for network reduction

A major difference in the waveforms obtained from Method 1 and Method 2 may be seen in the systolic part of the cycle, where the pressure obtained using Method 1 was overpredicted. This was observed as a general distinction between the two methods, and is exemplified in figure 4. However, the diastolic phase is very similar, indicating that the discrepancy is not a result of differences in the values of compliance in the WK3 models. The diastolic decay of pressure can be approximated by an exponential function, with an exponent given by the product of the peripheral resistance (R2) and the compliance (C) [31]. Thus changes in the compliance directly effect the diastolic shape, whereas changes in R1 only have secondary effect. On the other hand, R1 has a direct effect on the systolic part of the cycle. Inspection of the values used for the proximal resistance in the WK3 models revealed that Z (Method 1) was in general higher than R1,opt (Method 2) for the larger systemic arteries. The addition of the characteristic impedance to the original two-element WK model was based on frequency analysis of modulus and phase of the input impedance along the aorta. By including the characteristic impedance, the input impedance modulus of the modified WK matched in vivo measurements at high frequencies [2,31]. We also observed (not shown here) better matching of the modulus at the aorta for high frequencies, between baseline and reduced models obtained with Method 1 than with Method 2; however, the same is not true for this more distal location (internal carotid). Impedance phase, on the other hand, was captured better by Method 2 for some frequencies (particularly between 5 and 7 Hz), as can also be seen in the phase of the first minima of the flow waveform (≈ 6 Hz). Minimization of high-frequency oscillations has also been an incentive for using matched (R1 = Zc) WKs as outflow BC’s in 1D blood flow models [12]. However, the price to pay is an overprediction of pressure in systole.

Optimization of topological complexity

Central and larger systemic artery pressure waveforms

Pressure measured with a cuff and sphygmomanometer in the brachial artery is used routinely and accepted as an important predictor of future cardiovascular risk. However, studies indicate that central blood pressure (CBP) relates more strongly to cardiovascular events [32]. Systolic and pulse pressures are amplified as the pulse wave propagates through the larger systemic arteries. This amplification may vary significantly among subjects [32], making it difficult to map measurements of pressure at more peripheral sites directly to CBP.

Although it is still unclear if routine measurement/estimation of CBP will provide significantly improved risk stratification [33], the 1D nonlinear equations for blood flow can be used to investigate pulse wave amplification [8,34,35]. In previous studies, the topology of the 1D model was chosen ad hoc. Our novel framework provides a mathematical approach to determine the optimal topology to study pulse wave amplification from the aortic root to the brachial artery.

The results presented in the first two rows of figure 5 indicate that inclusion of detailed descriptions of upper and lower limbs are not needed in order to study pulse wave amplification from the aortic root to the brachial artery. Moreover, the entire cerebral circulation can be replaced by WK3 models with negligible effects on aortic and brachial pressure waveforms. This is reasonable since these are relatively small and stiff arteries for which the behaviour is well captured by WK3 models [3]; however, it is important to note that the proximal part of the aorta, which accounts for about 50% of total systemic compliance, needs to be kept in the reduced 1D model.

Both pulse pressure and augmentation pressure, and their relation (augmentation index) is associated with cardiovascular risk [36]. Even though the aortic pressure waveforms obtained by the reduced models in the top two examples in figure 5 captured the pulse pressure very well, some subtle deviations are visible in the systolic part of the waveforms. This could have an effect on the calculated augmentation pressure, and thus also on evaluations of cardiovascular risk. In the last example in figure 5, an error threshold of ε + ε of 0.7% at the aorta, was used, and results indicate that this 31-artery model captures the most important features of wave propagation for central aortic pressure. A similar model was found when a combined threshold of average ε of 0.4% was set for four arterial sites; midpoint of ascending aorta, right common carotid artery, right brachial artery and left femoral artery, as illustrated in figure 13 in appendix A.5. This network was also able to capture waveform features with good qualitative and quantitative precision when the model was re-parametrized to model different physiological and pathological states.

Carotid and cerebral circulation

In the study by Reymond et al. they compared carotid flow predictions with and without description of the cerebral circulation and stated that a detailed description was necessary in order to produce physiological correct waveforms. Our results, on the other hand, indicate that the entire cerebral circulation can be appropriately lumped into WK3 models effecting only the diastolic part of the flow waveforms and with , as shown in figure 6. Furthermore, by increasing the threshold to 3.4% the network is reduced to a very simplistic model including only five arterial segments. Though the overall features are represented in this five-artery model, the arterial tree is truncated close to the carotid artery and will thus be more influenced by the WK3 models. High-frequency details are not described well by the three-element WK [31], which in this case is visible through the smoothing of the second and third peaks of the flow waveform. Such errors were magnified when the model was transformed to represent normal ageing, as visualized in figure 14.

Figure 14.

Optimal reduced networks for flow at the distal end of the right carotid artery (left panel) and proximal end of the right middle cerebral artery (right panel). The networks are the same as shown in figure 6; however, parameters were altered in order to represent normal ageing as described in §2.6.

Figure 6 also shows results with flow rate at the inlet of the right middle cerebral artery set as the quantity of interest. This site is located more distal than the other quantities of interest studied in this work, and as can be seen in the case where a threshold of was considered, the circle of Willis can be ‘broken’ and represented by WK3 models without altering the flow waveform significantly. Furthermore, the arterial tree can be truncated in close proximity to the middle cerebral artery without introducing significant constraints on the solution, more so than was the case for the right common carotid artery. This is attributed to the fact that the flow in this region is more dominated by frictional forces resulting in pressure and flow waveforms that are of similar shape and phase and can be more readily described by the WK3 model. Moreover, by increasing the threshold to more of the larger systemic arteries may also be lumped, resulting in very simplistic descriptions of the arterial network that were still able to capture the main features of the flow waveform in the middle cerebral artery. For this model, however, errors were magnified when parameters were altered to represent different physiological states, indicating that having a reasonably complete description of the larger arteries is more important than including the nearby system of 1D model arteries.

Blood flow can be measured non-invasively by ultrasound in both the carotid and middle cerebral arteries; however, there are many sources of uncertainty and standard errors of measurements are normally higher than 10% [37]. In comparison, the modelling errors introduced by applying network reduction to obtain simpler descriptions of the arterial system were smaller.

Coronary pressure waveforms

Figure 7 shows the results from applying our methodology on the patient-specific coronary network. The model can be reduced to its most simplistic realization while still keeping the error for the predicted iFR on a level which is comparable with the standard deviation of repeated iFR measurements. The differences in predictions of FFR between baseline and reduced models were even smaller, and in fact smaller than the significant figures used in clinical decision-making. This is attributed to the fact that, unlike iFR, FFR is a cardiac cycle averaged quantity. Our approach for network reduction maintained the correct resistance throughout the domain, and thus also average flow and pressure distributions. The limited resolution of CCTA imaging contributes a layer of uncertainty since only features larger than approximately 1.0 mm can be resolved [38]. However, our results indicate that one should not necessarily strive to segment arteries down to this limit.

Concluding remarks

Our results have shown that to capture important features of the aortic pressure waveform, such as timing and shape of reflected waves, pressure augmentation and pulse pressure, a model with all aortic segments, but close to minimal description of the head and lower and upper limb arteries is sufficient. Furthermore, a detailed description of the cerebral circulation is not needed in order to capture physiologically correct waveforms in the common carotid and middle cerebral arteries. Even though our framework for network reduction was performed on a single set of parameters representing a normal physiological state, waveform features were also captured with good qualitative and quantitative precision when the models were re-parametrized to simulate different physiological and pathological states.

Our approach is targeted at computational models of the cardiovascular system, however, it should also be useful for the design of in vitro haemodynamic experiments. Such physical models are attractive tools for fundamental research on pulse wave propagation [30,39], and also play a key role in validating computational models [6]. Through further work, one could also imagine the relevance of our approach in the design of multi-scale models of the cardiovascular system, e.g. hybrid 3D–1D–0D models.