Lumped parameter model for hemodynamic simulation of congenital heart diseases

Abstract

The recent development of computer technology has made it possible to simulate the hemodynamics of congenital heart diseases on a desktop computer. However, multi-scale modeling of the cardiovascular system based on computed tomographic and magnetic resonance images still requires long simulation times. The lumped parameter model is potentially beneficial for real-time bedside simulation of congenital heart diseases. In this review, we introduce the basics of the lumped parameter model (time-varying elastance chamber model combined with modified Windkessel vasculature model) and illustrate its usage in hemodynamic simulation of congenital heart diseases using examples such as hypoplastic left heart syndrome and Fontan circulation. We also discuss the advantages of the lumped parameter model and the problems for clinical use.

Introduction

Hemodynamic management of patients with congenital heart diseases remains a challenge for pediatric cardiologists and pediatric cardiac surgeons, and the number of specialists for pediatric hemodynamic management is limited. Recent development of computer technology has allowed simulations of hemodynamics of congenital heart diseases on desktop computers. The performance of a recent desktop computer has already reached a sufficient stage to perform some kind of real-time simulation.

There are two different approaches to simulating cardiovascular systems of congenital heart diseases. One approach is multi-scale computational fluid dynamics (CFD) modeling, which is often combined with lumped parameter models [1]. Modern imaging technologies such as three-dimensional computed tomography (CT) and magnetic resonance imaging (MRI) provide specific anatomical information of complex heart anomalies [2]. Hemodynamic simulation based on these anatomical data will be helpful for surgical decision making [3]. However, because large-scale computations are required to use these data, execution of the simulation process is too long to meet the requirement of real-time clinical applications [4]. Another approach is the simulation with a lumped parameter model, which can be performed even on a small bedside computer [5]. Although this type of hemodynamic simulation does not include specific anatomical information, the greatest advantage of this approach is that real-time simulation can be performed while acquiring clinical data from patients. This approach is likely to be more helpful for perioperative hemodynamic management.

In this review, we introduce computational modeling of congenital heart diseases using the lumped parameter model and discuss its usefulness and issues in clinical use.

Lumped parameter model

The lumped parameter model of the cardiovascular system consists of two components; a time-varying elastance chamber component and a (modified) Windkessel vascular component. The time-varying chamber component describes instantaneous ventricular and atrial pressure–volume relations. The Windkessel vascular component describes hemodynamics in arterial and venous trees. The combination of these two components enables us to simulate the cardiovascular system; the left and right heart, and systemic and pulmonary circulations.

Time-varying elastance chamber component

In the lumped parameter model of cardiovascular system, each chamber (ventricle or atrium) is represented by a time-varying elastance model [6]. The basis of time-varying elastance chamber model has been reported by Suga et al. [7].

The chamber pressure at time (t) after onset of contraction is approximated as:

$$P(t) \cong E(t)[V(t) - V_{0} ],$$

(1)

where t is the time after onset of contraction, P(t) and V(t) are instantaneous chamber pressure and volume, respectively, and V ₀ is the volume axis intercept of the end-systolic pressure–volume relationship (ESPVR). E(t) is an instantaneous elastance of the chamber, which is independent of the loading conditions. At end-systole (t _es), E(t) reaches a maximal value (E _es), which is the slope of the ESPVR. Hence, ESPVR is described as:

$$P_{\text{es}} (t) = E_{\text{es}} [V_{\text{es}} (t) - V_{0} ],$$

(2)

where P _es and V _es are end-systolic pressure and volume, respectively. In contrast with ESPVR, end-diastolic pressure–volume relationship (EDPVR) is often represented by an exponential curve:

$$P_{\text{ed}} (t) = A\{ \exp [B(V_{\text{ed}} (t) - V_{0} )] - 1\} ,$$

(3)

where P _ed and V _ed are end-diastolic pressure and volume, respectively. A and B are constants. The parameter B is called as a stiffness constant.

Instantaneous chamber pressure is described by the sum of P _ed and the developed pressure (difference between P _es and P _ed) scaled by normalized elastance curve [e(t)]:

$$P(t) = [P_{\text{es}} (V(t)) - P_{\text{ed}} (V(t))]e(t) + P_{\text{ed}} (V(t)).$$

(4)

Parameters E _es, V ₀, A, and B are different among the four chambers (left and right ventricles and atria).

Elastance curve

In the real world, elastance curve (E(t)) looks like a skewed sine curve [8]. However, in computational simulation, the skewed sine curve is often replaced by a simple sine (or cosine) curve (Fig. 1a) [9, 10]:

$$\begin{array}{l} \,e(t) = 0.5[1 - \cos (\pi t/T_{\text{es}} )]\,(0 \le t < 2T_{\text{es}} ) \\ e(t) = 0\,(2T_{es} \le t < T), \\ \end{array}$$

(5)

where t is the time from the start of systole, T _es is the duration of systole, and T is the duration of cardiac cycle.

Fig. 1

Simulated normalized elastance curves (T _es = 175 ms). a A sine curve. b A sine curve combined with exponential curve. Time constant of relaxation (τ) is 25 ms for the solid line and 100 ms for the dotted line. An increase in τ prolongs the tail

A sine curve combined with an exponential curve is also used as a normalized elastance curve (Fig. 1b) [11]:

$$\begin{aligned} & e(t) = 0.5[1 + \sin (\pi t/T_{\text{es}} - \pi /2)]\,(0 < t < 3T_{\text{es}} /2) \\ & e(t) = 0.5\exp [(3T_{\text{es}} /2 - t)/\tau ]\,(3T_{\text{es}} /2 \le t), \\ \end{aligned}$$

(6)

where τ is the time constant of relaxation. An increase in τ prolongs the tail of the normalized elastance curve (Fig. 1b, dotted line). Because e(t) is discontinuous from the end of one cardiac cycle to the beginning of the next, this may cause an error in the simulation result.

Windkessel vascular components

In the lumped parameter model of cardiovascular system, vasculatures are simulated as Windkessel vasculature models. There are three major vasculature models: two-, three-, and four-element Windkessel models [12]. The two-element Windkessel model consists of a resistor and a capacitor (Fig. 2a). The resistor represents vascular resistance and the capacitor represents vascular volume capacitance. In the three-element Windkessel model, another resistor corresponding to characteristic impedance (R _c) is serially added to the two-element Windkessel resistor and capacitor (Fig. 2b). In the four-element Windkessel model, an inductor is added to the three-element Windkessel model (Fig. 2c). There are several variations in the position of the inductor in the four-element Windkessel model [13].

Fig. 2

Windkessel vasculature models. a Two-element model; b Three-element model; c Four-element model. R resistance, C capacitance, R _c characteristic impedance, L inductance

When we consider the input impedance, the two-element Windkessel model is less accurate than the three- and four-element Windkessel models especially in the high-frequency range corresponding to the characteristic impedance. Since aortic and/or pulmonary characteristic impedances in patients with congenital anomaly dramatically change between before and after the surgical treatment, three- or four-element Windkessel model may be suitable for hemodynamic simulation of congenital heart diseases.

At each resistance (R), a linear relation between pressure drop (ΔP) and flow (Q) applies, according to the Ohm’s law:

$$\Delta P = Q \cdot R.$$

(7)

At each capacitance (C), the relation between pressure (P) and volume (V), and the change in volume [dV(t)/dt] calculated by the difference between total inflow (∑Q _inflow) and total outflow (∑Q _outflow) are as follows:

$$P = \frac{V}{C},$$

(8)

$$\frac{{{\text{d}}V(t)}}{{{\text{d}}t}} = \sum Q_{\text{inflow}} (t) - \sum Q_{\text{outflow}} (t).$$

(9)

At each inductance (L), the relation between pressure drop ($\Delta P$) and change in flow [dQ(t)/dt] is as shown below:

$$\Delta P = L\frac{{{\text{d}}Q(t)}}{{{\text{d}}t}}.$$

(10)

Valves

Each valve (aortic, pulmonary, mitral, or tricuspid valve) is often represented as an ideal diode connected serially to a small resistor [14]. Several papers also used a diode connected with a non-linear resistance as a valve model [15], as described below:

$$\Delta P = k \cdot Q(t)^{2} .$$

(11)

Valve regurgitation is often simulated by introducing a second diode in the opposite direction with a small resistor [14].

Electrical analogs of cardiovascular systems

Using the framework of the time-varying elastance chamber model combined with the modified three-element Windkessel vasculature model, cardiovascular systems are described in terms of electrical analogs. Figure 3 shows the electrical analogs of the cardiovascular systems of a normal heart (Fig. 3a), Fontan circulation (total cavopulmonary connection, Fig. 3b), atrial septal defect (ASD, Fig. 3c), and ventricular septal defect (VSD, Fig. 3d).

Fig. 3

Electrical analogs of cardiovascular systems of a normal heart, b Fontan circulation, c patient with atrial septal defect (ASD), and d patient with ventricular septal defect (VSD). LA left atrium, LV left ventricle, RA right atrium, RV right ventricle, SA single atrium, SV single ventricle, MV mitral valve, TV tricuspid valve, AV aortic valve, PV pulmonary valve, AVV atrioventricular valve. R _CS, R _AS, and R _VS denote systemic characteristic impedance, arterial resistance, and venous resistance, respectively. R _CP, R _AP, and R _VP denote pulmonary characteristic impedance, arterial resistance, and venous resistance, respectively. C _AS and C _VS denote pulmonary arterial and venous capacitances, respectively. C _AP and C _VP denote pulmonary arterial and venous capacitances, respectively. R _ASD and R _VSD denote the resistances across ASD and VSD, respectively

In these electrical analogs, no non-linear components are used. However, non-linear components may sometimes be useful for hemodynamic simulation of congenital heart diseases as discussed below.

Non-linear components

In cardiovascular simulations of normal subjects and simple congenital heart diseases such as ASD and VSD, use of non-linear components may be unnecessary. However, hemodynamic simulation of some congenital heart diseases may require the use of non-linear components for modeling because shunt and stenosis behave non-linearly in clinical settings. To simulate these non-linear components, the Bernoulli’s principle for the total hydraulic energy of an ideal fluid is important [16]. The general form of the Bernoulli’s equation is;

$$P + \frac{1}{2}\rho v^{2} + \rho gh = {\text{constant,}}$$

(12)

where P is pressure, v is blood velocity, ρ is blood density, g is acceleration of gravity, and h is height. The term ρgh can be neglected in many cases of hemodynamic simulation because the height difference across a shunt or stenosis is small. Addition of this equation in modeling may improve the accuracy of hemodynamic simulation of congenital heart diseases.

However, to simulate the fluid dynamics in non-linear components more accurately, multi-scale CFD modeling and solution of the Navier–Stokes equations are needed [17]. Because the numerical method (for example, the finite element method) is necessary to solve Navier–Stokes equations, the simulation may take time and is therefore unfavorable for bedside simulation. On the other hand, the Navier–Stokes equation may be required when we want to couple a structural model with a fluid dynamics model. As an example, in the case of the aortic arch reconstruction of the Norwood procedure, which accompanies structural changes in the aortic arch, we may not be able to predict postoperative hemodynamics from a preoperative lumped model. In such a case, a multi-scale CFD model with the Navier–Stokes equation would provide more accurate prediction of postoperative hemodynamics.

Solution of simultaneous equations

To simulate hemodynamics, solution of simultaneous equations (Eqs. 2, 3, 4, one of Eqs. 5 and 6, Eqs. 7, 8, 9, 10, Eqs. 11 and/or 12 where necessary) with given initial parameters is necessary. The required initial parameters of chamber components are E _es, V ₀, A and B for each chamber, T _es, and a delay between atrial and ventricular systole, which can be calculated from heart rate [18]. Initial parameters of vascular components are R _c, R, L, and C for each of the four vasculatures: systemic arterial, systemic venous, pulmonary arterial, and pulmonary venous vasculatures. In addition, the initial stressed blood volume, which is the sum of the initial volumes of all chambers and capacitors, is necessary. If non-linear components are used, these initial parameters are also needed to start the solution.

The simplest method to solve simultaneous equations is the Euler method, which is the first-order numerical procedure for solving differential equations with given initial values. When the first derivative of y(t) (dy/dt) is described as a function of time (t) and y(t),

$$\frac{{{\text{d}}y}}{{{\text{d}}t}} = f(t,y(t)),$$

$$y(t_{0} ) = y_{0} ,$$

where t ₀ is the time of start and y ₀ is the value of y at t ₀. The time after n steps (t _n ) is described using the step size (h):

$$t_{n} = t_{0} + nh.$$

For the Euler method, one step from t _n to t _n+1 is defined as:

$$\begin{aligned} & y_{n + 1} = y_{n} + hf(t_{n} ,y_{n} ) \\ & t_{n + 1} = t_{n} + h. \\ \end{aligned}$$

(13)

Although this method is easy for programming, the accuracy is relatively low.

Local truncation error for the Euler method is proportional to h ², while that for the fourth-order Runge–Kutta method is h ⁵. Therefore, fourth or higher order Runge–Kutta method may be favorable for the solution of simultaneous equations including non-linear components. Commercial software such as MATLAB/Simulink (The MathWorks, Inc. MA, USA) may be useful for solution using the high-order methods.

Hemodynamic simulation of congenital heart diseases

In this section, we introduce several examples of lumped parameter models used to simulate hemodynamics of congenital heart diseases such as hypoplastic left heart syndrome and Fontan circulation.

Modeling of hypoplastic left heart syndrome

Since hypoplastic left heart syndrome is one of the most serious congenital heart anomalies, many researchers are interested in hemodynamic simulation of hypoplastic left heart syndrome. Hemodynamic simulations of the Norwood procedure, a stage I palliation for hypoplastic left heart syndrome, are widely performed because the mortality of the procedure remains high (approximately 16%) even in 2017 [19] and hemodynamic assessments using computational models may improve its clinical outcome.

Migliavacca et al. [18] modeled the cardiovascular system after the original Norwood procedure with the modified Blalock–Taussig (BT) shunt, using a lumped parameter model combined with a non-linear BT shunt model. They examined the effects of shunt size, vascular resistances, and heart rate on the hemodynamics and oxygenation after the original Norwood procedure. Their simulation demonstrates that (1) larger shunts divert an increased proportion of cardiac output to the lungs, away from systemic perfusion, resulting in poorer O₂ delivery, (2) systemic vascular resistance exerts greater effect on hemodynamics than pulmonary vascular resistance, and (3) systemic arterial oxygenation is minimally influenced by heart rate changes.

In the Sano modification of the Norwood procedure, the BT shunt is replaced by a right ventricle-to-pulmonary artery shunt [20]. The Sano modification was also simulated using a lumped parameter model [21]. The electrical analog and simulated pressures of the Sano modification are shown in Fig. 4. In this simulation, the Sano modification, even using a non-valved conduit, is preferable to the original Norwood procedure from the viewpoint of ventricular energetics.

Fig. 4

a The electrical analog of the modified Norwood procedure with right ventricle (RV) to pulmonary artery (PA) shunt (Sano modification). The pressure gradient across the RV–PA shunt is described as a non-linear function. b Simulated aortic (black solid line), right ventricular (gray solid line) and pulmonary artery pressures (gray dashed line) of the Sano modification. AP aortic pressure, RVP right ventricular pressure, PAP pulmonary artery pressure

Researchers also have an interest in another type of stage I palliation; bilateral pulmonary artery (PA) bandings combined with ductal stenting (the so-called hybrid procedure). Young et al. [22] developed a lumped parameter model of the hybrid procedure and examined the effects of diameters of PA bandings and ductal stent on ventricular workload. In their model, larger PA banding diameter or a small stent diameter below 7 mm substantially increases ventricle workload and reduces systemic perfusion. Simulations were also used to compare the hybrid procedure and the Norwood procedures from the viewpoint of ventricular energetics [23]. Mechanical efficiency of the hybrid procedure is equivalent to that of the original Norwood procedure, but inferior to that of the Sano modification.

Modeling of Fontan circulation

Fontan operation is a surgical goal for patients with single ventricle physiology [24]. Although heart failure after the Fontan operation requires heart transplantation or ventricular assist device (VAD) implantation, VAD implantation in Fontan patients remains a challenging issue in clinical settings [25]. Therefore, several researchers have investigated the effects of VAD on failing Fontan circulation [26,27,28]. Furthermore, reduced exercise capacity in Fontan patients is also a clinical interest, although the mechanism is not fully understood. Therefore, researchers also focus on exercise physiology in Fontan patients [29, 30].

Di Molfetta et al. [26] investigated the effects of cavopulmonary assistance (right-sided VAD) on the Fontan circulation using a lumped parameter model. They also examined the effects of combined use of continuous flow and pulsatile VAD on the Fontan circulation, and reported that left-sided continuous flow VAD concurrent with right-sided pulsatile VAD maximizes the hemodynamic benefits [27].

The effects of partial cavopulmonary assistance have also been examined using a lumped parameter model combined with a non-linear VAD model [28]. Partial cavopulmonary assistance from the inferior vena cava (IVC) to the pulmonary artery maintains cardiac index at lower IVC pressure but increases the superior vena cava pressure substantially, when compared with total cavopulmonary assistance from inferior and superior vena cavae to pulmonary artery.

Kung et al. [29] developed a closed loop lumped parameter model for simulation of exercise physiology in Fontan patients. Their model successfully reproduced the average exercise response of a cohort of typical Fontan patients. Based on the simulation results using a lumped parameter computational model, Koeken et al. [30] reported that exercise capacity in Fontan patients is limited due to an increase of central venous pressure and the incapability to further reduce systemic resistance, consequently restricting systemic flow.

Modeling of other congenital heart diseases

Effectiveness of the one and a half ventricle repair for hypoplastic right ventricle remains controversial in clinical settings. Several criteria of patient selection are based on anatomical factors such as Z score of the tricuspid valve and the right ventricular volume [31], but there are no accurate criteria based on physiological factors. A simulation study using a lumped parameter model reveals that the right ventricular stiffness constant (B) may predict the effectiveness of one and a half ventricular repair in improving postoperative hemodynamics [32].

The tetralogy of Fallot is one of the most common congenital heart anomalies, and is characterized by a ventricular septal defect, pulmonary stenosis, an overriding aorta and right ventricular hypertrophy. Kilner et al. [33] used a lumped parameter model to examine the factors influencing pulmonary regurgitation after repair of tetralogy of Fallot. They reported that pulmonary regurgitation was exacerbated by increased pulmonary artery compliance and by elevated pulmonary arteriolar resistance.

Advantages of hemodynamic simulation using the lumped parameter model

The major advantage of the lumped parameter model is that it requires shorter execution time even using a desktop computer. Although recent development of computer technology have enabled us to simulate hemodynamics of congenital heart diseases on desktop computers, simulation using a multi-scale CFD model is still time-consuming, and could take several days to complete [4]. In contrast to the multi-scale CFD model, the lumped parameter model may allow real-time simulation. This means that the lumped parameter model can be used to simulate patients’ specific hemodynamics at the bedside under real-time data acquisition from the patients. Broomé et al. [34] reported the possibility of building a clinically relevant real-time computer simulation model of the normal adult cardiovascular system. Although several problems of applying a real-time computer simulation model to complex congenital anomalies remain to be solved, the usefulness of lumped parameter model in supporting clinical decision-making has been advocated [35].

Problems of hemodynamic simulation using the lumped parameter model

The lumped parameter model has the potential of simulating patient-specific hemodynamics at the bedside. However, there are several problems in clinical application. One of the most critical problems for bedside simulation is how to identify patient-specific (initial) parameters such as end-systolic elastance, resistances, inductances and compliances in clinical settings.

To identify systolic and diastolic chamber properties such as end-systolic elastance (E _es) and the stiffness constant (B), simultaneous monitoring of chamber pressure and volume is necessary. However, simultaneous measurement of chamber pressure and volume in clinical settings may be difficult even in adult patients without congenital heart anomaly. On the other hand, the recent development of imaging technologies such as knowledge-based reconstruction of three-dimensional echocardiographic images with MRI images may allow monitoring of chamber volume [36]. Combined with pressure data in the catheter laboratory, these volume data may be used in the estimation of chamber properties.

The three-element Windkessel parameters can be estimated in a time-domain and/or a frequency-domain manner from vascular pressure and flow waveforms [37, 38]. Furthermore, the method of beat-to-beat estimation of peripheral resistance (R) and arterial compliance (C) during transient conditions has been reported [39]. A technique to estimate the four-element Windkessel parameters for transient and steady beats has also been developed [40]. However, in contrast to the three-element Windkessel model, it may be difficult to identify the four-element Windkessel parameters, especially the inductances (L) [41]. Segers et al. reported that the four-element Windkessel model yielded the best quality of fit, but model parameters reached physically impossible values for L in about 12% of all cases [42]. Although an iterated unscented Kalman filter algorithm may improve estimation accuracy of four-element Windkessel parameters [43], further investigations are necessary to identify L with sufficient accuracy in every patients.

The direct estimation of venous capacitances and resistances is difficult in clinical settings as well as animal experiments. Therefore, several assumption of these values may be required at the start of the simulation and an algorithm to tune these values to physiologically proper values may be necessary.

In addition, estimation of all the parameters may take much time and may be invasive for patients with congenital heart diseases. Therefore, an algorithm that automatically estimates lumped parameters would be beneficial for clinical application of lumped parameter models. Schiavazzi et al. [44] attempted to develop a framework for automated tuning of lumped model parameters to match clinical data of patients after the Norwood procedure.

Instantaneous autonomic outflow at the estimation may affect lumped model parameters. Therefore, the results of hemodynamic simulation may depend on baseline autonomic nervous activities. In addition, global impairment of cardiac autonomic nervous activity is reported in patients after the operation of congenital heart diseases [45, 46]. How to incorporate the influence of autonomic nervous activities into the model remains an important matter. Furthermore, in children, body growth may also affect lumped model parameters. Therefore, re-identification of the parameters is necessary along with body growth. However, the model of cardiovascular system itself may be reusable.

Conclusions

The lumped parameter model is potentially beneficial for hemodynamic simulation of congenital heart diseases at the bedside because it allows real-time simulation. However, how to obtain patients’ specific parameters remains a large obstacle for bedside simulation. Bedside simulation of congenital heart diseases can be realized if the problem of parameter identification can be solved.