General overview of the model
The model (Figure 1) is based on the following compartments: the four cardiac chambers and corresponding valves, the pericardium, the systemic circulation consisting of the aortic root, ascending aorta, proximal aortic arch, right carotid/subclavian artery, distal aortic arch, left carotid/subclavian artery, descending aorta and the representation of a peripheral arterial system (See also electrical analogue in Additional file 1: Figure S1). The peripheral arterial system contains resistance arteries (arteriolae), capillaries and small veins emptying into the inferior caval vein, while the two carotid/subclavian arteries in a similar way end up in the superior caval vein. The pulmonary circulation is represented by the main pulmonary trunk, pulmonary resistance arteries (arteriolae), capillaries, small veins and pulmonary veins emptying into the left atrium. Two coronary vessels are also included. Input to the vascular model are physical dimensions [12, 13], tissue properties and number of parallel vessels in each compartment modified to fit the anatomical structure of the present vascular model; while pressure, flow, volume and oxygen saturation are simulated outputs and can be displayed. The cardiac, valvular, pericardial and vascular properties can be modified in detail by changing the parameters presented in the supplement (Additional file 2: Table S1-5). Model complexity is chosen to enable future simulations of central hemodynamics in major congenital and acquired circulatory pathologies and extra-corporeal circulatory support.
Cardiac model
The heart chambers
The four cardiac chambers passive and active properties are modeled with separate time-varying elastance functions based on the “Double-Hill equation” as previously described [11, 14] including a slightly modified Starling mechanism (Figure 2) [8]. The model includes internal chamber flow resistance and viscous wall properties of the heart, in accordance with previous models [9, 11]. Details are shown in the cardiac model section of the supplement.
By using a continuous function through the cardiac cycle manual setting of valvular timing [9] is avoided. Each cardiac chamber has its own elastance function of realistic shape (Additional file 3: Figure S3). Timing and duration of contractions in the cardiac chambers can be set independently, in order to create variable atrioventricular (AV) and interventricular (VV) intervals as well as realistic timing of valve opening and closure. Thereby the model offers additional opportunities to simulate arrhythmias as well as interventricular dyssynchronicity. As an example atrial fibrillation can be simulated by turning off atrial contractions and randomizing the time between ventricular contractions.
The electrocardiogram (ECG) (Additional file 4: Figure S4; Additional file 5: Figure S5) is simulated by summing the four cardiac elastance functions with a linear weight factor, based on the approximate distance from each cardiac chamber to the chosen surface ECG lead assuming an electromechanical delay of 60 ms. The ECG does not by itself influence the simulated physiology and is included as an output of the model only to add clinical realism during educational sessions, therefore the detailed algorithms are not described.
The heart valves
Valve pressure gradients are composed of a Bernoulli resistance and an inertial term (Additional file 2: Equation S4 and Figure S6 of Additional file 6) as earlier described [8, 9]. Minimum and maximum valvular areas can be set to simulate valvular stenosis and regurgitation. The valve areas change gradually (Additional file 4: Figure S4) as a function of pressure gradients and valve inertial constants [11] as is shown in the supplement (Additional file 2: Table S4). Blood flow inertance L is determined by density ρ, the variable inflow length l and the valve area A (Equation 1). Inflow length l is set to a value identical to the instantaneous diameter of the valve in contrast to the constant value in Mynard et al.[11].
(1)
The atrial and ventricular septum
The interventricular septal displacement is simulated using pressure transmission through the muscular septum as previously described [9, 15] (Equation 2). As a further development the ventricular septal stiffness E
sv
is modeled in order to increase in proportion to the left ventricular systolic elastance e
lv
(t) (Equation 3), thereby stiffening during contraction, in accordance with known physiological characteristics. Additional file 5: Figure S5 shows the ventricular septal shift at different stiffness.
(2)
(3)
Atrial septal interaction is modeled in a similar way. The septal interaction does not give any significant contribution during normal physiology in apnea, but is important when simulating interventricular dyssynchrony, variations in intrathoracic pressures including artificial ventilation and changes in ventricular loading conditions.
The pericardium
In order to resemble pericardial function in health and disease a user-defined exponential function relating intra-pericardial pressure to total heart volume is adopted from the literature [9] (Additional file 7: Figure S7). As a further development of the previous model the minimal pericardial pressure is allowed to be a negative value in agreement with what is found experimentally when measured during intrathoracic pressure changes or hypovolemia [16]. Further details can be found in the supplement.
Intracardiac and arteriovenous shunts
Atrial septal defects, ventricular septal defects and persistent ductus arteriosus flows are calculated according to Equation 4, with instantaneous flow Q proportional to the area A of the defect, the square-root of the pressure difference ∆P and the Gorlin empirical constant C assumed to be 1 [17]. Acceleration g due to gravity is set to 980 cm/s2.
(4)
Vascular model
The systemic circulation
Transmural pressure p in each vascular compartment is related to the variable segmental vascular volume v according to an exponential relation (Additional file 2: Equation S7, Additional file 8: Figure S12) [9] creating an increase in stiffness with progressive distension of the vascular wall.
Resistance R
0
, inertia I
0
, and elastance E
0
at normal mean pressure P
0
are calculated from published and estimated vascular properties [12, 13], and relations [18] described in Equations 5–7. These vascular properties are all approximately inversely proportional to n, which is the number of parallel vessels, lumped into each compartment [18]. The parameter η is blood viscosity, l, r
0
and h is the length, radius and thickness of a single vessel respectively and ρ is blood density.
(5)
(6)
(7)
The same value of Young’s modulus Y
inc
is used in every vascular segment, implying that vascular stiffness properties are solely explained by this constant and vascular dimensions. The actual volume-dependent segmental elastances, resistances and inertias are updated in each calculation step based on the volume and radius assuming constant vessel length and thickness. Details can be found in the blood vessels section of the supplement (Additional file 2: Table S5-7).
The resistance Ω, in series with the capacitor situated in parallel to the blood flow in every vascular segment, corresponds to the viscous property of the vascular wall, dampening the flow pulse. The parameter value assigned to this component is calculated as in Equation 8 with addition of a scaling factor λ as compared to the characteristic impedance in a classical 3-component Windkessel model [19, 20].
The same formula was used to calculate Ω in every vascular segment as previously suggested [19], although the name characteristic impedance is less well suited in a distributed model. A scaling factor λ between 0.0 and 1.0 is used to tune the vascular damping properties. The effects of changing the scaling factor λ can be seen in Figure 3.
The systemic capillaries and veins are described by the same principal vascular model as the arteries although inertial effects are of less importance due to lower pulsatility [21] and a larger cross-sectional area. Every single vessel has high resistance and stiffness; however the large amounts of vessels in parallel lumped into each compartment result in both a low resistance (Equation 5) and a low elastance (Equation 7), corresponding to high compliance and being related to their physiological volume reservoir function. The distribution of the blood volume between arteries, capillaries and veins in the systemic and pulmonary circulation in the model is shown in Additional file 2: Table S7.
The pulmonary circulation
The pulmonary artery, resistance vessels, capillaries and veins are modeled as in the systemic circulation. An adjustable pulmonary shunt mimics the circulation without gas exchange due to anatomical anomalies or ventilation-perfusion disturbances as in atelectatic lung tissue.
The coronary circulation
Coronary circulation is simulated with a left and right coronary artery emptying in the right atrium as described in the supplement (See also Additional file 9: Figure S8). Myocardial oxygen consumption is calculated according to the pressure-volume area concept of Suga et al.[22] on a beat-by-beat basis and estimations of myocardial oxygen balance can therefore be calculated based on both supply and demand.
Oxygen transport
Blood oxygen content in every compartment is calculated as a constant multiplied by the hemoglobin (Hb ), the oxygen saturation (Sat ) and the volume of the compartment (Volume ) as in clinical routine and most catheterization laboratories [23] (Equation 9).
(9)
Oxygen saturation within each compartment is considered homogenous and oxygen transport between compartments is proportional to blood flow, hemoglobin level and saturation. Calculations concerning blood mixing are simplified as effects of blood jets, incomplete mixing and physically dissolved oxygen not are taken into account.
Baroreceptor reflex
A baroreceptor reflex can be activated in the model that affects heart rate, cardiac contractility (maximum elastance) and arterial vascular resistance as in Sun et al.[9] and described in the baroreceptor section of the supplement.
Simulation
The simulator was compiled as a stand-alone software developed in Visual Basic.NET 2010 (Microsoft Corporation, Redmond, WA, USA) and .NET Framework 4.0 (Microsoft Corporation, Redmond, WA, USA). A time step of 0.25 ms was used to assure stability. The implemented 62 differential equations were solved with an implicit Euler numerical method [24]. Mean values were calculated as running arithmetical means to simplify calculations and conserve memory resources. The program version used was Aplysia CorVascSim 4.4.0.62 (Aplysia Medical AB, Stockholm, Sweden). The model can be run in real-time on a standard personal computer (example: Intel processor Core i3 3.06 GHz, RAM 3 Gb, Graphics ATI Radeon HD4550, Windows 7 32/64-bit).
Parameter settings in simulations of normal physiology, systolic heart failure, diastolic heart failure, aortic stenosis, aortic regurgitation, a Valsalva maneuver, exercise and progressive arteriosclerosis
All hemodynamic data were presented with zero intrathoracic pressure in order to resemble end-expiratory values. With the basic set of parameters presented in the supplement an adult man weighing 70 kg is described with a blood volume of 5600 ml (80 ml/kg), blood viscosity of 0.00024 mmHg · s and blood density of 1.060 g/cm3. Systolic left ventricular heart failure was simulated by decreasing the maximum elastance from 2.8 to 1.0 mmHg/ml. Diastolic left ventricular heart failure was simulated by increasing the basal passive diastolic elastance E
min
from 0.05 to 0.12 mmHg/ml. Aortic stenosis was simulated by decreasing the open aortic valve area from 5.0 to 0.7 cm2. Aortic regurgitation was simulated by increasing the closed aortic valve area from 0.0 to 0.2 cm2. Simulation of a Valsalva maneuver was accomplished by increasing intrathoracic pressure from 0 to 10 mmHg for approximately 20 seconds with and without baroreceptors activated. All other parameters and settings were identical.
The stepwise simulation of exercise was accomplished by increasing heart rate from 72 to 144 min-1, increased left and right ventricular contractility by 50% and increased systemic and pulmonary resistance arterial diameters by 50%. The cardiac elastance parameters α
1
and α
2
were scaled with the square-root of the factor changing the heart rate to preserve a realistic relation between systolic and diastolic times. The vascular damping scaling factor λ was set to 0.5 in all vascular segments and presented simulations. The oxygen saturation in the pulmonary capillaries was set to 99.4%, corresponding to a normal alveolar pO2 of 13.3 kPa and the pulmonary shunt was set to 10% of cardiac output in all presented simulations. Hemoglobin level was set to 140 g/l and systemic oxygen consumption to 250 ml/min in addition to myocardial oxygen extraction. A Young’s modulus Y
inc
of 3000 mmHg (= 0.40 MPa) similar to the value in Wang et al.[12] was used to calculate the vascular elastance E
0
at the vessel-specific normal mean pressure P
0
(see Additional file 2: Table S5) in all simulations above. A progressive arteriosclerotic process was simulated by changing vascular stiffness in all vascular segments through variations in Young’s modulus Y
inc
between 2000 and 6000 mmHg. Note that comparisons with real clinical cases will not be performed due to unknown patient-specific secondary adaptations.
Sensitivity analysis
The relation between changes in input parameters x and model output variables y was studied with a sensitivity analysis during simulation of normal physiology (Additional file 10). Input parameters were increased 10% one by one. Seven model output variables [left ventricular end-systolic pressure (LVESP), right ventricular end-systolic pressure (RVESP), left atrial pressure (LAP), right atrial pressure (RAP), left ventricular stroke work (LVSW), right ventricular stroke work (RVSW) and cardiac output (CO)] focusing on cardiac function were evaluated after 60 seconds assuring a hemodynamic steady-state. Sensitivity S was calculated according to Equation 10.