Towards Patient-Specific Modeling of Coronary Hemodynamics in Healthy and Diseased State

A model describing the primary relations between the cardiac muscle and coronary circulation might be useful for interpreting coronary hemodynamics in case multiple types of coronary circulatory disease are present. The main contribution of the present study is the coupling of a microstructure-based heart contraction model with a 1D wave propagation model. The 1D representation of the vessels enables patient-specific modeling of the arteries and/or can serve as boundary conditions for detailed 3D models, while the heart model enables the simulation of cardiac disease, with physiology-based parameter changes. Here, the different components of the model are explained and the ability of the model to describe coronary hemodynamics in health and disease is evaluated. Two disease types are modeled: coronary epicardial stenoses and left ventricular hypertrophy with an aortic valve stenosis. In all simulations (healthy and diseased), the dynamics of pressure and flow qualitatively agreed with observations described in literature. We conclude that the model adequately can predict coronary hemodynamics in both normal and diseased state based on patient-specific clinical data.


Introduction
Diagnosis of coronary circulatory disease based on coronary pressure and/or flow has been shown to improve the clinical outcome of treatment [1]. The direct measurement of coronary hemodynamics, however, is still limited to the large epicardial vessels, which means that microvascular disease can only be determined from proximal measurements using an appropriate model of the vessels and their interaction with the cardiac muscle [2].
Due to the location of the coronary arteries, being embedded into the myocardium, the contraction of the heart influences coronary hemodynamics. This results in the unique feature that blood is mainly supplied during diastole, while coronary epicardial pressure is high in systole. Several mathematical models have been proposed to model the effect of cardiac contraction on the coronary vessels. Downey and Kirk [3] proposed the vascular waterfall mechanism, which explained the increased resistance to blood flow by vascular collapse when the intramyocardial pressure, determined by ventricular cavity pressure, exceeds the lumen pressure. Spaan et al. [4] introduced the intramyocardial pump model, which accounted for the role of vascular compliance. The intramyocardial pressure, determined by the ventricular cavity pressure, served as the extravascular pressure. The timevarying elastance concept [5] was first applied to the coronary circulation by Krams et al. [6], in which flow is impeded due to a varying stiffness of the cardiac wall. Other, more elaborate, models also include the effect of the coronary vessels on the cardiac contraction (e.g., [7]).
One-dimensional (1D) wave propagation models are potentially useful in the interpretation of arterial hemodynamics [2,8]. They can provide better insight into the effect of combinations of epicardial and microvascular disease on clinically used indices. Application of 1D wave propagation models to the coronary circulation has been reported in several studies. Smith et al. [9] and Huo and Kassab [10] applied 1D models to vessels representing the coronary arterial anatomy. However, in these studies the interaction with the cardiac muscle was not taken into account. The effect of the cardiac contraction was taken into account in the 1D models of Mynard and Nithiarasu [11], where the coronary vessels were loaded with an approximated left ventricular pressure. Recently, as part of an elaborate arterial tree, Reymond et al. [12] coupled a 1D model to a timevarying elastance model to describe the relation between ventricular pressure, volume, and coronary hemodynamics.
The aim of this study was to construct a model of the human cardiovascular system in which clinical data can be incorporated, to enable patient-specific modeling of coronary hemodynamics. The different submodels were chosen such that the model complexity remains minimal, while still enabling the incorporation of both normal and diseased physiologic and geometric clinical data. To represent the heart, the single-fiber contraction model [13,14] was chosen. This model is based on microstructural material and macrostructural geometrical properties, allowing the simulation of cardiac disease with geometry-based parameter changes. The large vessels are modeled one-dimensionally [15] to enable easy implementation of the geometry of the vessels, whereas the small vessels are represented by lumped elements.
Here, the different components of the model are explained and the ability of the model to describe both normal and pathological coronary pressure and flow dynamics is evaluated via comparison with observations described in literature. Two disease types are modeled: coronary stenoses located in the epicardial vessels and left ventricular hypertrophy with an aortic valve stenosis, affecting the coronary microvasculature.

Materials and Methods
The model consists of three main elements: a heart contraction model, a wave propagation model for the large arteries, and lumped elements to model the coronary and systemic microcirculation. As part of the wave propagation model, Bessems [16] developed an element that can describe the effect of a stenosis on the local hemodynamics. Since the wave propagation model and heart model have already been described in Bessems et al. [15] and Bovendeerd et al. [14], respectively, only a short description of the models is given below.

Heart Contraction
Model. Similar to Bovendeerd et al. [14], the left ventricle is modeled as a thick-walled sphere, consisting of nested spherical shells. When assuming rotational symmetry and homogeneity of mechanical load, the relation between tissue stress and ventricular pressure ( lv ), cavity volume ( lv ), and wall volume ( ) can be described as: Here, is the fiber stress and is the radial wall stress at lv = lv , the shell located at one third of the ventricular wall. This representative shell was chosen because the strain at this location is similar to the fiber strain [14]. At this location, assuming incompressibility of the myocardial tissue, the fiber stretch ( ) and the radial stretch ( ) can be related to the ventricular geometry by [13]: with the instantaneous sarcomere length and 0 the sarcomere length at lv,0 ; the cavity volume at zero transmural pressure.
The myofibers are modeled one-dimensionally, exerting only stress in the fiber direction. The fiber stress consists of an active ( ) and passive ( ) stress component, where only depends on the sarcomere length ( ), while also depends on the sarcomere shortening velocity ( ) and time elapsed since activation ( ): The active stress is modeled according to Kerckhoffs et al. [18], which describes a combination of the contractility ( ) and three functions: 1 relates the active stress to the sarcomere length and is given by: where it is assumed that the fibers can only exert stress in tension. 0 and are a scaling and curvature parameter, respectively. 0 is the sarcomere length at which the stress is zero. The time dependent activation function 2 is defined as: Here, max is the activation duration and and are the activation rise and decay time constant, respectively.
The dependency of the active stress on the sarcomere shortening velocity is modeled hyperbolically: with 0 the unloaded shortening velocity and the curvature of the hyperbolic relation.

Computational and Mathematical Methods in Medicine 3
Womersley velocity profile profile Approximated velocity Figure 1: A schematic representation of the exact velocity profile (left) and the approximation (right). is the approximated core radius and the viscous layer. Adapted from Bessems [16].
The passive stress in the fiber ( ) and radial ( ) direction are modeled in a similar way: The passive stress-length relation is determined by the scaling parameters 0 and 0 and the curvature parameters and , respectively. The intramyocardial pressure ( im ) is used as the extravascular pressure on the coronary circulation and is assumed to be linearly dependent on the radial position in the wall. The shell at lv = lv is also considered representative for im : with ,lv and ,lv the inner and outer radius of the ventricle, respectively.

Valves.
Since the atrial contraction was not taken into account, the behaviour of the mitral valve was simplified and modeled as an ideal diode, where the pressure gradient over the mitral valve (Δ mv ) is determined by Ohm's law: with mv the flow through the mitral valve and mv defined as: mv, and mv, are the resistance of the valve in the open and closed situation, respectively. For the aortic valve the inertia is taken into account and opens due to a positive pressure gradient (Δ av ) and closes when the flow through the valve ( av ) becomes negative. The differential equation relating Δ av to av is defined as: Here, av is the resistance and av the inertia of the valve. av is determined by the cross-sectional area , the effective valvular length av , and blood density ; av = av / av . The value of the resistance parameter av is determined by the state of the valve: Here, av, and av, are the resistance of the aortic valve in the open and closed situation, respectively. The values of the different parameters are listed in Table 1.

Wave Propagation
Model. The governing equations describing the one-dimensional propagation of pressure and flow waves of a Newtonian incompressible fluid are derived from the conservation of mass and momentum and were taken from Bessems et al. [15]. The conservation of mass is derived as: with and representing the axial direction and time, is the local arterial lumen area, is the volumetric flow rate, and the flow per unit length distributed to small side-branches that are not separately modeled by vascular segments. As described by Bessems et al. [15], an appropriate velocity profile function is assumed that describes the frictional forces and non-linear terms in the balance of momentum equation (see Figure 1): Here, is the transmural pressure, is the vessel radius, is the body force, that is, force per unit mass, in the axial direction, and and are the dynamic viscosity and density of the fluid, respectively. The wall shear stress is given by: with = (max[0, (1 − √ 2/ )]) 2 representing the relation between the radius of the inertia dominated core and the thickness of the Stokes layer. is the Womersley parameter according to = √2 / , with the heart rate. The parameter also determines the parameter of the convective term in (15):

Arterial Wall Model.
To solve (14) and (15) a constitutive relation between and is required. In this section we will derive an analytical description of the coronary arterial wall mechanics, with parameters that depend on the geometry of the vessel only and are based on a microstructural constitutive model. The main advantage of this approach is that this   [17]. The fiber orientation is determined by angle . (b) A two-dimensional representation of the stenosis element. is the length of the stenosis, the radius of the vessel at the neck of the stenosis, and 0 the radius of the vessel proximal to the stenosis. Adapted from Bessems [16].
analytical model enables easy implementation into the wave propagation equations, while retaining the microstructural properties of the arterial wall. In van der Horst et al. [20], it was demonstrated that a model based on the two-fiber constitutive model developed by Holzapfel et al. [17] was able to accurately capture the radius-pressure relations of porcine and human coronary arteries. In that model the arterial wall was modeled as a cross-ply of helically wound fibers embedded in a cylinder composed of hyperelastic material. The stress ( )-stretch ( ) behaviour is determined by: with ℎ the hydrostatic pressure, I the unity tensor, the shear modulus, B the Finger tensor, and the fiber stress of fiber . ⃗ is the fiber orientation, which is represented by the angle with the circumference (see Figure 2(a)). It was assumed that in compression no stress can be transmitted by the fibers, with defined as: Here, is the collagen fiber stretch and 1 and 2 are constants determining the stress-stretch relation of the collagen fibers. For coronary arteries the value of these four parameters have been determined [20]. With removal of the outliers, the median of the four parameters are: = 19.3 kPa, comply with this rule. The second rule emerges from the finding that the fiber orientation at physiological loading ( phys ) is almost constant for all arteries. Using this rule, can be directly related to the circumferential stretch at = 13.3 kPa, via two constants: 0 and phys = 36.4 ∘ .
The balance equations resulting from this two-fiber model are solved using numerical integration. Since a numerical integration scheme is also employed to solve (14) and (15), the direct implementation of this two-fiber model will be computationally expensive. Therefore, as an intermediate step, a phenomenological model described by Langewouters et al. [22] is used to analytically relate the compliance ( ) to the pressure ( ): Here, 0 , 1 , , and determine the -relation. To relate this analytical model to the two-fiber model, it is assumed that these four parameters depend on the only clinically measurable quantities: the radius ( ) and wallthickness (ℎ ) at physiological axial stretch and pressure ( = 13.3 kPa). First, the -relation was determined with the two-fiber model (including the optimization rules) for different combinations of and ℎ . ranged from 0.25 to 3 mm and ℎ ranged from 0.025 to 0.3 mm. Then, for each combination of and ℎ within the range 0.05 < < 0.15 ( = ℎ / ), the four parameters of the Langewouters model were fitted with the Gauss-Newton algorithm as implemented in Matlab (R2010a, The Mathworks, Natick, MA). Using the multiple regression function in Statgraphics (Centurion XVI, statpoint technologies, inc. Warrenton, VA) a polynomial was extracted based on the best 2 -adjusted value: The parameters ) are determined by = 2 and and the correlation is good, as indicated by the 2 values. With these relations, the compliance of the coronary arteries as function of the instantaneous diameter could be determined and used in the wave propagation model.
For the systemic arterial walls we use a simple linear elastic model: Here, is the Young's modulus and is Poisson's ratio. For all systemic arteries, incompressibility is assumed ( = 0.5) and = 0.4 MPa [19].

Stenosis Element.
While one-dimensional theory is suitable to model the pressure and flow waves in relatively straight arteries, it may yield unrealistic results in pathological regions like aneurysms and stenoses. In the derivation of the one-dimensional model it is assumed that variations in the cross-sectional area of the vessels is relatively small, so the radial and circumferential blood velocity is negligibly small with respect to the axial component. Considering that epicardial coronary arteries are prone to the development of stenoses, it is necessary to use a specific element that can be incorporated into to the 1D model. Bessems [16] developed a 1D stenosis element, based on the semi-empirical relations obtained by Young and Tsai [23,24] but with an improved contribution of the viscous and unsteady components, to calculate the pressure-drop over an axisymmetric narrowing. The parameters of this model are based on two-dimensional axisymmetrical finite element simulations of stenotic hemodynamics. The axisymmetric stenosis geometry is depicted in Figure 2 From oscillatory flow simulations Bessems [16] derived the following relation for the pressure drop over a stenosis (Δ ): The parameters 0 and are the cross-sectional areas of the vessel proximal to and at the neck of the stenosis, respectively. is the average flow, and , , , and are empirically determined constants obtained by Bessems [16]. They are given by: is the resistance and is inertia across the stenosis: 0 is the radius of the vessel proximal of the stenosis and ( ) the varying radius of the vessel at the site of the stenosis. When assuming that the pressure drop develops linearly over the length of the stenosis, (23) can be written in the following differential form: The conservation of mass in the stenosis is given by (14) and its compliance is assumed to be negligible.

Lumped Elements.
The contribution of the peripheral vasculature at each end-point of the 1D model is lumped with the three-element model depicted in Figure 3 (see [25]). The relation between the pressure and flow for this model can be written as: ] .
The wave impedance is given by: with and the average cross-sectional area and compliance of the connecting vessel. The total resistance ( + ) was determined from the average flow and pressure drop over the element. Finally, is the compliance of the artery defined by a time constant = , with = 2 s.  Table 2. The LMCA has a length of 5 mm and splits into the LAD and the LCx, with a length 7.5 cm and 6 cm, respectively. Side branches are modeled at intervals of 1.5 cm. Each coronary segment is represented by the characters a-f. The radius of segment a is 1 mm and Murray's law is used to determine the radius of segments b-f. All a-segments are connected to the three-element model representing the coronary microvessels. The intramyocardial pressure ( im ) acts on the three capacitors that represent the vessel compliance. When a stenosis is modeled, it is incorporated into the -segment of the LAD. resulting pressure gradient the flow over the aortic valve is calculated using (12) and serves as the input for the aortic wave propagation elements.
To include the effect of the systemic wave reflections, the aorta is modeled with all major side branches, with at each distal end a terminal impedance that is prescribed using the three-element model. The geometrical data are taken from Stergiopulos et al. [19] and are listed in Table 2.
A hypothetical anatomy of the main coronary branches is assumed. The left main coronary artery (LMCA) and right coronary artery (RCA) originate 5 mm distal from the aortic valve. The LMCA splits into the left anterior descending (LAD) and circumflex (LCx) arteries. The LAD has a total length of 7.5 cm with four side branches (representing the diagonal and septal side branches), while the LCx has a length of 6 cm with three side branches (marginal and posteriorlateral branches). The geometry of the RCA is equal to the LAD and it is assumed that it supplies both the left and right ventricle (RV), with a ratio of 0.4. Since the RV Table 2: Geometric and physiological parameters of the arterial vessels. The length ( ), proximal radius ( ), distal radius ( ), and wall thickness (ℎ) of the systemic arteries are based on Stergiopulos et al. [19]. The parameters of the three-element model (Z, , and ; see Figure 3) were determined as described in Section 2.3. The numbers of the vessels correspond to the numbers shown in Figure 4. is not modeled here, it is assumed that the intramyocardial pressure ( im ) of the RV is smaller by a factor proportional to the ratio of maximum pressure in the two ventricles ( im,rv = 0.2 im,lv ). For all terminal coronary 1D vessels (14 in total) a radius of 1 mm was prescribed and for each bifurcation Murray's law [26] relates the diameter of the parent and daughter vessels by a power of 3. Based on van den Broek et al. [27], it is assumed that for all coronary vessels the wall thickness is equal to 10% of the radius ( = 0.1).
The microvasculature is modeled with three serial threeelement models. The total resistance ( ) is determined using Ohm's law by assuming an average pressure of 100 mmHg and prescribing a flow of approximately 20 mL/min through every terminal branch. Following Bovendeerd et al. [14], is distributed over the four resistances according to: art = (7/27) , myo1,2 = (9/27) , and ven = (2/27) . The values of the three capacitors are based on measurements by Spaan [28]: art = 0.2 mm 3 Pa −1 , myo = 0.53 mm 3 Pa −1 , and ven = 0.65 mm 3 Pa −1 . The intramyocardial ( im ) pressure that is generated by the heart contraction model is connected to the three capacitors to model the extravascular pressure on the coronary vessels. Since the circulation is not closed, a constant venous pressure of 700 Pa (±5 mmHg) is prescribed at the output of the model.
To be able to solve the full set of equations, the method of Kroon et al. [29] is used in which the 1D and 0D pressure and flow relations of (14), (15), and (27) are solved fully coupled by writing the differential equations in the same form. The 1D vessels are divided into a number of nonoverlapping elements of 2.5 mm and the temporal discretization is performed using the Euler implicit integration scheme. The final set of equations is solved using a direct solver [30], as implemented in the finite element package Sepran (Ingenieursbureau SEPRA, Leidschendam, The Netherlands).

Simulations and Data Analysis.
To test whether the model is able to describe coronary hemodynamics in both normal and pathological situations, three different simulations are performed. The normal, healthy situation is modeled with the parameters listed in Tables 2 and 1. The pressure, flow, and volumes of the heart and aorta obtained with the model are qualitatively compared to similar signals described by Van De Vosse and Stergiopulos [8]. The modeled coronary pressure and flow in the LMCA and RCA are compared to pressure and velocity measurements performed simultaneously in the LMCA and RCA by Hadjiloizou et al. [21]. Two types of pathological situations are modeled: a stenosis in the LAD (see Figure 4) and left ventricular hypertrophy due to an aortic valve stenosis (LVH-AVS).
Dynamic pressure measurements proximal and distal to a mild stenosis (50% diameter, length 2.65 mm) and severe stenosis (70% diameter, length 7.48 mm), performed at the catheterization laboratory of the Catharina Hospital (Eindhoven, The Netherlands), are used to verify the stenosis element. Since the flow through the stenotic vessels was not measured, quantitative comparison is difficult. Therefore, the dynamics of the measured and modeled pressure signals is only compared qualitatively. Since the pressure measurements are performed during hyperaemia, the flow in the model was increased five-fold, by decreasing the coronary microvascular resistances. Furthermore, the clinically most used index to quantify coronary stenoses, fractional flow reserve (FFR), which is defined as the ratio of the pressure distal and proximal to a stenosis, is determined with both the model and the measurements.
The ability of the model to describe coronary hemodynamics when LVH-AVS is present is verified with clinical measurements performed by Hildick-Smith and Shapiro [31]. With transthoracic Doppler echocardiography, they measured the dynamics of flow in the LAD in LVH-AVS patients before and six months after aortic valve replacement (AVR). While the average left ventricular cavity volume was constant before and six months after AVR, the measured average ventricular mass decreased significantly: from 271 to 226 g. The average aortic valve pressure gradient before AVR was 93 mmHg and systemic pressures were normal with a minimum and maximum pressure of 89 and 134 mmHg. The two situations before and after AVR are modeled with ventricular wall volumes based on the measured ventricular wall mass, assuming a mass density of 1.1 kg/L. The contractility and aortic valve resistance are increased such that the pressure gradient across the aortic valve is approximately 93 mmHg, while the average aortic pressure remains normal. The model parameters of the heart contraction model and valves are listed in Table 1. The main features of the dynamics of the modeled flow in the LAD are qualitatively compared to the measurements by Hildick-Smith and Shapiro [31].
The difference between the arterial wall model derived in Section 2.2.1 and a simple linear elastic model (22), with respect to the hemodynamics, is investigated. As it is expected that the difference between the used coronary arterial wall model and a simple linear elastic model is largest in the low pressure range, the difference between the two models is determined both proximal and distal to the severe stenosis described above. The relative differences between the pressure, flow, cross-sectional area, and wall shear stress calculated with the model are quantified with parameter : Here, is the hemodynamic signal obtained with the microstructural-based coronary arterial model and lin is the signal obtained with the linear elastic model with a Young's modulus of 1.5 MPa and a Poisson ratio of 0.5. Figure 5 shows that the heart and systemic pressures, flows, and volumes obtained with the model, qualitatively agree with values found in literature [8]. With a stroke volume of 70 mL/min, a mean aortic flow of 4.8 L/min, and an aortic mean and pulse pressure of 93 and 47, respectively, the main clinically relevant parameters are within the normal physiological range. At the transition between systole and diastole the effect of the closing of the aortic valve together with a reflection originating from the bifurcation to the iliac arteries is also clearly visible. The dynamics of the different signals are similar, except the timedependent behaviour of the mitral flow, especially at late diastole (Figures 5(f) and 5(c)). This is obviously due to the lack of the atrial contraction in the model. Another clear difference is that modeled aortic pressure in early systole increases faster than found in Van De Vosse and Stergiopulos [8], which is related to the choice of heart activation function.

Normal Hemodynamics.
The modeled pressure and fleow in the left main (LMCA) and proximal right (RCA) coronary artery are depicted in Figure 6. The data are compared to pressure and blood velocity measurements acquired simultaneously in a human LMCA and RCA [21]. It should be noted that the mean and pulse pressure measured by Hadjiloizou et al. [21] were relatively high and there was an average offset of approximately 15 mmHg between the pressure measured in the RCA and LMCA (Figure 6(a)). Although these pressures are not considered to be representative for non-diseased vessels, the data do enable the qualitative comparison between the flow velocities in both the LMCA and RCA and the pressureflow relation. The modeled pressures in the LMCA and RCA were almost identical and are determined by the aortic pressure. The flow in the LMCA was diastolic dominated, with the typical flow impediment during early systole. The ratio of maximum diastolic and systolic flow in the LMCA was 2.1 for both the simulation and measurements. Although it depends to what degree the RCA supplies blood to the left or right ventricle, the flow in the RCA was markedly less dominant in diastole, compared to the LMCA. For the RCA, the ratio of maximum diastolic and systolic flow was 1.2 and 0.9 for the measurements and simulation, respectively. The difference between the LMCA and RCA demonstrates the influence of the intramyocardial pressure on the coronary flow. Besides the pressure in early systole, the main difference between the simulated coronary hemodynamics and the clinical data is that the flow in the early diastole displays a peak in the simulation with subsequently a relatively large decline, whereas this is not the case in the experimental data. Figures  7(a) and 7(b) demonstrate the effect of a mild stenosis (50% diameter) and severe stenosis (70% diameter), respectively. It is obvious that the pressure gradient was much larger for the severe stenosis, especially in diastole when the flow was highest. The pressures determined with the model showed the same behaviour as the measurements, with the largest pressure gradient in diastole (Figures 7(c) and 7(d)). For the modeled mild stenosis this pressure gradient variation between systole and diastole is larger than in the experimental data. A possible explanation for this discrepancy is a less distinct difference between systolic and diastolic flow in that measurement. The FFR values determined with the

Left Ventricular Hypertrophy with an Aortic Valve Stenosis.
In Figure 8 the effect of LVH-AVS on the coronary flow is shown and compared to transthoracic Doppler echocardiography measurements by Hildick-Smith and Shapiro [31], before and six month after AVR. The normal characteristic flow dynamics, with a small positive systolic and large diastole component, was found after AVR, in both the measurements and the simulations. Before the AVR, so when the LVH-AVS is present, the measurements reveal that the positive systolic flow component was replaced by a period of negative flow. The measured maximum diastolic flow decreased slightly after AVR. The ratio of the maximum positive or negative systolic velocity before and after AVR was −1.3, whereas ratio of the maximum diastolic velocity before and after AVR was 1.1. These features were also captured by the model, demonstrating the influence of the increased intramyocardial pressure on the coronary flow dynamics. The ratio of the simulated maximum positive or negative systolic velocity before and after AVR was −0.6, whereas the ratio of the maximum diastolic velocity before and after AVR was 1.3.

Arterial Wall
Model. The effect of using the arterial wall model proposed by Langewouters et al. [22] compared to a linear elastic model will be most apparent in the low pressure range. Therefore, the difference between the two arterial wall models was determined both proximal and distal to the 70% diameter stenosis. In Figure 9 the relative difference between the two models, as defined by (30), on the pressure ( ), flow ( ), cross-sectional ( ), compliance ( ), and wall shear stress ( ) are shown. Proximal to the stenosis, the difference between the two models was small. While the effect on the pressure (max = 5%), flow (max = 2%), and crosssectional area (max = 8%) was rather limited, the wall shear stress changed significantly (max = 17%).

Discussion
In the present study, previously published models of the heart and vessels have been combined to create a model capable of describing coronary hemodynamics in health and disease.
By coupling a heart model to a 1D wave propagation model, the effect of heart disease on both the coronary microvessels and the aortic perfusion pressure could be related to coronary epicardial hemodynamics. With the combination of models, stable solutions were obtained and the waveforms found with the model featured the main characteristics of both systemic and epicardial coronary pressure and flow dynamics. Additionally, by changing a limited amount of parameters, a coronary stenosis and left ventricular hypertrophy with an aortic valve stenosis (LVH-AVS) could be modeled and produced specific hemodynamical features that qualitatively agreed with experimental observations described in literature. The heart mechanics is governed by the single-fiber contraction model developed by Bovendeerd et al. [14]. The main advantage of this particular model over existing models (e.g., the intramyocardial pump model [4] and the timevarying elastance model [5]) is that it is based on geometric data that can be obtained in the clinic in combination with microstructural properties of the myocardium. This enables the simulation of cardiac disease with physiologybased parameter changes, as was shown by simulating LVH-AVS. Being modeled as a sphere with myofibers oriented in the same direction in each shell, the heart model is a simplified representation of the cardiac muscle. Although the validity of this model should be evaluated for each type of cardiac disease, this simplified representation is also the strength of model, since it is able to produce physiological hemodynamics with a limited amount of parameters. Due to the use of a representative intramyocardial pressure, and average values for the coronary compliances and resistances, the coronary flow in the model should be regarded as a mean flow over the myocardium. It therefore cannot describe radial layer-specific differences in coronary perfusion, which can be clinically relevant in relation to ischemia. Although it will increase the number of parameters, these spatial differences can be incorporated by modeling branches at different layers in the myocardium. Even though the main features of coronary hemodynamics were captured, in future studies it would be interesting to incorporate these branches at different transmural positions to get a more physiologic representative flow distribution in the myocardium. The effect of the deformation of the vessels on its compliance and resistance, especially on the venous side, should then be taken into account as well [32,33]. From the comparison between the model and literature it was found that the activation model used results in a rise in pressure that is too fast in early systole. An activation function as proposed in van der Hout-van der Jagt et al. [34], in which the early and late part of the activation function can be tuned separately, might prove to resolve this issue. However, this does increase the number of parameters. The contraction of the left atrium was not modeled, which was clearly reflected in the mitral valve flow. This, however, did not result in an unrealistic pressure-volume relation in the left ventricle. The right ventricle was also not incorporated into the model. Therefore, the right ventricular intramyocardial pressure ( im,rv ) was approximated by a factor (0.2) proportional to the left ventricular intramyocardial pressure ( im,lv ). To get a more realistic measure of im,rv , the right ventricle can also be modeled with a similar heart contraction model as was demonstrated by Cox et al. [35]. The systemic large epicardial coronary arteries are modeled one-dimensionally, which enables the investigation of the propagating pressure and flow waves as was validated by Bessems et al. [15]. This specific model has the advantage that it is time domain-based and has a velocity profile that approximates the actual Womersley profiles. For the relatively small coronary arteries with Womersley numbers of approximately 2 the velocity profiles are almost similar to the Poiseuille profile, whereas in the aorta there is a phase difference between velocities near the wall and in the central core. This similarity to Womersley profiles is also important when the 1D model is used as the boundary condition for a more detailed 3D model. The choice which branches are lumped or modeled individually, depends on the point of interest and the specific disease that being modeled and the available clinical data. In the model presented here the total aorta and the first part of its side-branches were modeled individually to include its main reflection sites. The coronary arteries with a radius smaller than 1 mm are lumped, since intracoronary measurements are mainly limited to the larger coronary vessels. In future studies, the 1D representation of the coronary vessels also enables the analysis of pressure and flow wave patterns, which have been the subject of recent research [21,36].
A three-element model was chosen to represent the coronary microvasculature. While this representation did result in physiological coronary hemodynamics, a four-element model with an inertia term [37] might improve the signal, particulary at the large increase in flow during early diastole where the inertia of the blood will play a role. Furthermore, it was found that the parameters of the first three-element model have a large influence on the dynamics of the coronary flow. A proper sensitivity analysis of the model parameters may be helpful in the correct parameter choice for patientspecific modeling.
The compliance of the arterial wall of the coronary vessels was modeled with the analytical model of Langewouters et al. [22]. The parameters of the model were fitted to the model described in van der Horst et al. [20]. For different radii and wall thicknesses an accurate, polynomial description was found for each parameter of the Langewouters model. The main advantage of this approach is that the analytical description enables easy implementation into the model, at low additional computational cost, while the microstructural

13
properties are taken into account. By comparing the pressure and flow waves obtained with this wall model with a linear elastic model, it was found that the differences where very small, even distal to a severe stenosis where the change in compliance are the largest. The wall shear stress, however, did change significantly distal to this stenosis (17%). This might be of clinical interest, since the wall shear stress has been indicated as a factor involved in the development and destabilization of plaques [38]. The stenosis element has been shown to be compatible with the wave propagation elements and agrees qualitatively with pressure measurements from the clinic. The flow in these stenotic vessels was, however, not measured, which makes a proper quantitative comparison impossible. The stenosis is modeled as being smooth and axisymmetric, whereas in clinical practice stenoses are irregular. This might also be the reason why the measured FFR values were lower than the ones obtained with the model. Although 3D modeling [39] is required to investigate to which extent the stenosis element is capable of describing the hemodynamics of irregular stenoses, it is likely that in most cases the stenosis element cannot adequately describe stenoses found in patients. Besides the shape-dependency of the constants of the stenosis element equation, it is also assumed that these constants depend on the heart frequency only, neglecting the contribution of other frequency components. Although Bessems [16] verified the relation with finite element simulations for physiologic flows and found only small differences, this is obviously a limitation of the stenosis model.
The simulated LVH-AVS also qualitatively agreed with the data described in literature, indicating that model is able to capture the global effect of LVH-AVS on coronary flow. However, for a proper verification of the model, simultaneous measurements of left ventricular pressure and volume and coronary pressure and flow should be performed in both nondiseased and LVH-AVS hearts. Furthermore, a number of case studies should also be performed, in which the effect of for different disease types on coronary hemodynamics is measured under well controlled conditions. An isolated beating heart set-up [40] might be a suitable platform for these studies.
The next step in improving the model would be to include autoregulatory mechanisms. The baroreflex mechanism could be included to regulate the heart rate, as was already incorporated into a similar heart contraction model by Cox et al. [35]. Furthermore, by including the coronary autoregulation, the difference between resting and hyperaemic flow can be simulated [41]. This is valuable since it enables the determination of clinical indices based on the difference between resting and hyperaemic hemodynamics, for example, coronary flow reserve [42] or diastolic coronary vascular reserve [43]. This heart contraction model is suitable to include this mechanism, since the work performed by the heart can be used as a parameter in the autoregulation mechanism.
Besides application of this model to enhance the diagnosis in case of combinations of multiple disease types, the model can also be used to investigate the global effect of an intervention, bypass surgery, or collaterals on coronary epicardial hemodynamics. To be able to use the model for patient-specific modeling of the (diseased) coronary circulation, model parameters need to be fitted to hemodynamical measurements. Besides adequate measuring devices for coronary pressure and flow (and lumen area), this will require a proper parameter sensitivity analysis of the model. Due to the number of parameters involved, a Monte Carlo approach as used by Huberts et al. [44] might be suitable to find the most influential parameters.

Conclusion
We constructed a model of the cardiovascular system, in which physiologic and geometric clinical data can be incorporated for patient-specific modeling of coronary hemodynamics. The modeled pressure and flow dynamics are in qualitative agreement with clinical measurements described in literature, especially with respect to the shape details. Although further research is required to improve and verify the model, we conclude that the model adequately can predict coronary hemodynamics in both normal and diseased state based on patient-specific clinical data.