Analysis of Operating Modes for Left Ventricle Assist Devices via Integrated Models of Blood Circulation

Abstract

Left ventricular assist devices provide circulatory support to patients with end-stage heart failure. The standard operating conditions of the pump imply limitations on the rotation speed of the rotor. In this work we validate a model for three pumps (Sputnik 1, Sputnik 2, Sputnik D) using a mock circulation facility and known data for the pump HeartMate II. We combine this model with a 1D model of haemodynamics in the aorta and a lumped model of the left heart with valves dynamics. The model without pump is validated with known data in normal conditions. Simulations of left ventricular dilated cardiomyopathy show that none of the pumps are capable of reproducing the normal stroke volume in their operating ranges while complying with all criteria of physiologically feasible operation. We also observe that the paediatric pump Sputnik D can operate in the conditions of adult circulation with the same efficiency as the adult LVADs.

1. Introduction

Left ventricular assist devices (LVADs) provide a therapeutic option to treat patients with end-stage heart failure (HF). LVAD connects the left ventricle (LV) and the aortic arch (AA), provides pulsatile or continuous blood flow and maintains circulatory support. The total aortic flow is the sum of the LV and LVAD outflows. Thus, LVAD decreases the work of the LV on ejecting blood to the aorta. The outflow from the LVAD to the aorta depends on the pressure drop over the LVAD. The pressure drop over LVAD is a complex interplay of many factors including LV contraction and ejection, aortic valve function, aorta extensibility and outflow to the distal parts of the systemic arteries. Modern LVADs are the rotary blood pumps (RBPs) which produce continuous flow to maintain temporary and permanent circulatory support [1,2]. The area of the pressure–volume (P-V) loop of the LV represents its stroke work. Dynamic head pressure-bypass flow (H-Q) curves characterise the RBP function during the cardiac cycle. These two curves correlate with each other [3]. They help to analyze dynamic interaction between the LV and RBP [4]. The other useful parameters of such analysis are the stroke work, the hydraulic pump work, and the cardiac mechanical efficiency [5,6].

Clinical efficacy of LVADs has been recently proven [7], although their impact on the cardiovascular system is not always clear. Algorithms of autonomous optimal control for LVADs are still widely discussed. A lot of information on LVADs operation in different regimes is available from mock circulation facilities [8,9,10]. Blood flow near aortic valve after implantation of LVAD was simulated in [11]. However, systematic data of the impact of pump operation in patient’s physiological conditions are limited due to complexity of measurements and a relatively small number of observable cases. In this work we develop an in silico model of the left heart and aorta with LVAD which allows us to simulate the impact of a pump in realistic physiological conditions. Our primary goal is to study behaviour of RBPs Sputnik 1, Sputnik 2 and Sputnik D [8,9,10,12,13] in patient’s physiological conditions. We also compare models of the Sputnik devices with a model of the HeartMate II [14], the well known and widely used pump.

In Section 2 of this work, we present an integrated model of the left heart function with aortic and mitral valves dynamic opening and closing (see Section 2.3). The model includes two segments of the aorta which are simulated by a 1D haemodynamic model (see Section 2.2). We identify the LVAD model for pumps Sputnik 1, Sputnik 2, Sputnik D by fitting parameters with data from mock circulation facilities (see Section 2.1). The parameters of the LVAD HeartMate II model are set according to the literature [14]. The LVAD model is included in the integrated model as a nonlinear lumped compartment which connects the left ventricle and the aorta. In Section 3.1, the heart and aorta model is validated in healthy conditions using data from the literature. Essential conditions of the physiologically feasible pump function include the temporary opening of the aortic valve, the positive direction of the flow through the LVAD as well as absence of ventricular suction and recovery of the typical stroke volume [13]. In Section 3.2, we study the haemodynamic effect of every pump in case of HF accompanied with left ventricular dilated cardiomyopathy in a range of rotational speed which covers ranges defined by manufacturers. We observe that in the considered conditions none of the pumps are capable of restoring the normal stroke volume in the ranges recommended by the manufacturers and at the same time complying with all criteria of physiologically feasible operation. We show that although Sputnik D was initially designed for paediatric patients, it can operate in the conditions of adult circulation at a higher pump speed, with about the same efficiency as adult LVADs. In Section 4 we discuss the results, limitations, conclusions and our future work.

2. Materials and Methods

2.1. Identification of Pump Models

Head pressure–flow rate (H-Q) relationship is a mechanical characteristic of a pump which can be determined from the laboratory tests. It provides a convenient interface for incorporating the pump model to a model of the cardiovascular system as a nonlinear compartment. The two options for deriving the H-Q relationship are the usage of (semi-)empirical formulas and the derivation from the physical principles. In rather general form $H(Q,\omega )$ is a quadratic form which is sometimes extended with the terms of the flow and pump rotation accelerations $H(Q,\omega ,\frac{dQ}{dt},\frac{d\omega }{dt})$. Here H is the head pressure, Q is the flow through the pump, $\omega$ is the rotation speed of the pump rotor. The Euler head equation with added quadratic term $Q^2$ due to experimental evidence and the flow inertia term [15] gives Model 1:

$$H=a{Q}^2+bQ+c{\omega }^2+d\frac{dQ}{dt}.$$

(1)

A similar model [16] with the rotational acceleration of the pump defines Model 2:

$$H=a{Q}^2+bQ+c{\omega }^2+d\frac{dQ}{dt}+e\frac{d\omega }{dt}.$$

(2)

The steady-flow model based on the conservation laws of mass, momentum and energy [17] yields Model 3:

$$H=a{Q}^2+bQ\omega +c{\omega }^2.$$

(3)

An addition to Model 3 of unsteady-flow effects and periphery parts [14] results in Model 4:

$$P_{lv}-P_p=a{Q}^2+bQ\omega +c{\omega }^2+d\frac{dQ}{dt}+P_{rec}-H_{per},$$

(4)

where $P_p$ is the pressure in the junction of the pump outlet and aorta, $P_{lv}$ is the pressure in the LV,

$$P_{rec}=\{\begin{array}{cc}0,\hfill & Q>e\omega \hfill \\ R_{rec}{(Q-e\omega )}^2,\hfill & Q⩽e\omega \hfill \end{array},$$

(5)

$$H_{per}=-L_{per}\frac{dQ}{dt}+R_{per}{Q}_p|Q_p|.$$

(6)

We note that Model 4 is not actually an H-Q relationship. It also includes parameters of the external (periphery) part which connects the pump to the LV and the aorta. The physical background of (4) is discussed in [14]. The theoretical Euler head equation gives the terms proportional to ${\omega }^2$ and $Q\omega$. The fluid friction losses produce quadratic growth $(Q^2)$ with the flow elevation. The flow detachment at the leading and trailing edges of the blade produces eddy and separation losses proportional to ${\omega }^2$, $Q\omega$ and $Q^2$. Part-load recirculation in the blade channels occurs below the design flow rate. It partly blocks the channels, decreases their effective diameter and increases the head pressure introducing ${(Q-e\omega )}^2$ term. The flow inertia term is proportional to $\frac{dQ}{dt}$. Fluid friction and inertia frequency-dependent losses in the peripheral part are included via $H_{per}$ in (6). See [14] and references herein for more details.

All the pump models (1)–(4) have physical interpretation. They were successfully validated with experimental data from different pumps in [14,15,16,17]. We take all these models as possible candidates for the H-Q mathematical relationship of non-pulsatile axial flow LVADs Sputnik D, Sputnik 1 and Sputnik 2. We use data from laboratory experiments with physical models of the paediatric mock circulation with Sputnik D [8,9,10] and the adult mock circulation with Sputnik 1 and Sputnik 2 [8,12] for validation. Sensors are placed as close to the pumps as possible, thus, we exclude peripheral term $H_{per}$ from (4) at the model fitting stage. Experimental setup imitates physiological conditions, including the Frank-Starling autoregulation mechanism of the heart which regulates the cardiac output depending on the ventricle preload. The 32% aqueous glycerol solution was used as the model fluid. Head pressure–flow rate (H-Q) curves for Sputnik D, Sputnik 1 and Sputnik 2 were measured at various constant pump speeds. For Sputnik D, the data from the range $6\times {10}^3$–$12\times {10}^3$ rpm with the step ${10}^3$ rpm were used as the training dataset, and the data from the range $13\times {10}^3$–$15\times {10}^3$ rpm were used as the test dataset. For Sputnik 1 and Sputnik 2 the data from the range $5\times {10}^3$–${10}^4$ rpm with the step 200 rpm and a contractility factor of the artificial LV 0.25 were used as the training dataset and the data from the range $5\times {10}^3$–${10}^4$ rpm and the contractility factor of the artificial LV 0.5 were used as the test dataset. The contractility factor [8] is a coefficient which decreases the end-systolic elasticity.

We set head pressure H as a target variable. The parameters of the models were identified by the damped least-squares method (Levenberg–Marquardt algorithm) [18,19]. We smooth up the raw data by Savitzky–Golay filter [20] for computing time derivatives of the flow and the rotational speed of the pump. The coefficient of determination $R^2$ was used as the best-fit criterion. According to results presented in

Table 1. Coefficients of determination ($R^2$) for models (1)–(4).

	Sputnik D	Sputnik 1	Sputnik 2
Model 1	$-0.2$	0.92	0.93
Model 2	$-0.19$	0.92	0.94
Model 3	0.87	0.37	0.39
Model 4	0.96	0.93	0.97

, Model 4 provides the best fit with experimental data for all Sputnik pumps.

Table 2. Parameters of model (4).

Parameter	Unit	Sputnik D	Sputnik 1	Sputnik 2	HeartMate II [14]
a	mmHg/(L/min)${}^2$	0.48	$3.58\times {10}^{-2}$	$-0.46$	$-0.86$
b	mmHg/(rpm · L/min)	$-1.52\times {10}^{-3}$	$-9.86\times {10}^{-4}$	$-5.64\times {10}^{-4}$	3.21 $\times {10}^{-4}$
c	mmHg/rpm${}^2$	$8.49\times {10}^{-7}$	$1.74\times {10}^{-6}$	$1.73\times {10}^{-6}$	$9.54\times {10}^{-7}$
d	mmHg·s${}^2$/L	$-60.06$	$-83.50$	$-85.91$	$-22.97$
e	(L/min)/rpm	$4.92\times {10}^{-5}$	$-2.18\times {10}^{-4}$	$-3.70\times {10}^{-4}$	$3.59\times {10}^{-4}$
$R_{rec}$	mmHg/(L/min)${}^2$	5.63	4.15	5.59	3.07
$L_{per}$	mmHg·s${}^2$/L	19.33	19.33	19.33	20
$R_{per}$	mmHg/(L/min)${}^2$	0.35	0.35	0.35	0.38

comprises identified parameters of Sputnik pumps for Model 4, as well as Model 4 parameters of the LVAD HeartMate II from [14]. Due to the lack of experimental data for Sputnik pumps periphery, we use mean values of the corresponding parameters for HeartMate II from [14].

In the following sections, we incorporate Model 4 into a lumped model of the heart coupled with a 1D model of the aorta.

2.2. 1D Mathematical Model of the Blood Flow in Aorta Segments

The blood flow in the aorta is simulated by a 1D reduced-order model of unsteady flow of viscous incompressible fluid in elastic tubes. The aorta is divided into two segments. The 1D model of the aorta is connected to the LV at the inlet, to the Windkessel compartment at the outlet and to the pump compartment between its segments I and II (see Figure 1).

Reviews and details of 1D haemodynamic models can be found in [21,22,23,24,25]. Algorithms of patient-specific parameter identification of such models were suggested in [26,27]. In this section, we briefly present this approach. We consider two 1D segments of the aorta which correspond to two parts of the ascending aorta (I and II in Figure 1). We assume that the pump is connected to the aorta at the aortic arch before the carotid arteries. The flow in every vessel is described by mass and momentum conservation equations

$$\frac{\partial \mathbf{V}}{\partial t}+\frac{\partial \mathbf{F}(\mathbf{V})}{\partial x}=\mathbf{G}(\mathbf{V}),$$

(7)

$$\mathbf{V}=\left(\begin{array}{c}A\\ u\end{array}\right),\mathbf{F}(\mathbf{V})=\left(\begin{array}{c}Au\\ u^2/2+p(A)/\rho \end{array}\right),\mathbf{G}(\mathbf{V})=\left(\begin{array}{c}0\\ \psi \end{array}\right),$$

where t is the time, x is the distance along the vessel counted from the vessel junction point, $\rho$ is the blood density (constant), $A(t,x)$ is the vessel cross-section area, p is the blood pressure, $u(t,x)$ is the linear velocity averaged over the cross-section, $\psi$ is the friction force

$$\psi =-8\pi \mu \frac{u}{\rho A},$$

(8)

$\mu$ is the dynamic viscosity of the blood. The elasticity of the vessel wall material is characterised by the $p(A)$ relationship

$$p(A)=\rho c_0^{2}f(A),$$

(9)

where $c_0$ is the velocity of small disturbances propagation in the vessel wall, $f(A)$ is the monotone S-like function (see [28] for the review of other options)

$$f(A)=\{\begin{array}{cc}\exp (\eta -1)-1,\hfill & \eta >1\hfill \\ \ln \eta ,\hfill & \eta ⩽1\hfill \end{array},\eta =\frac{A}{A_0},$$

(10)

$A_0$ is the cross-sectional area of the unstressed vessel.

The mass conservation condition at the aortic root includes the blood flow through the aortic root $Q_{ao}$ which is also a variable of the heart model from Section 2.3:

$$u_I(t,0)A_I(t,0)=Q_{ao}(t).$$

(11)

Boundary conditions at the connection of aorta and the pump include mass conservation condition

$$u_I(t,L_I)A_I(t,L_I)+Q_{pump}=u_{II}(t,0)A_{II}(t,0)$$

(12)

and the continuity of the total pressure

$$p_I(A_I(t,L_I))+\frac{\rho u_I^2(t,L_I)}{2}=p_{II}(A_{II}(t,0))+\frac{\rho u_{II}^2(t,0)}{2}=p_p+\frac{\rho }{2}{\left(\frac{Q_p}{S_p}\right)}^2,$$

(13)

where $p_p$ is the static pressure at the output of the pump, $Q_p$ is the flow through the pump contributing to (4), $S_p$ is the cross-section area of the tube which connects the output of the pump and the aorta.

The outflow boundary conditions assume that the terminal part of the aorta is connected to the Windkessel compartment which describes the rest of the systemic circulation

$$\frac{dQ}{dt}=\frac{1}{R_1}(\frac{d{p}_{II}(A_{II}(t,L_{II}))}{dt}-\frac{d{p}_{WK}}{dt}),$$

(14)

$$\frac{d{p}_{WK}}{dt}=\frac{Q}{C}(1+\frac{R_1}{R_2})-\frac{d{p}_{II}(A_{II}(t,L_{II}))-p_{\infty }}{R_{2}C},$$

(15)

$$Q=u_{II}(t,L_{II})A_{II}(t,L_{II}),$$

(16)

where $R_1,R_2,C,p_{\infty }$ are parameters presented in

Table 3. Parameters of the 1D model of haemodynamics in the segments of the aorta.

Parameter	Unit	Value	Parameter	Unit	Value
$A_{0,I}$	cm${}^2$	7.1	$A_{0,II}$	cm${}^2$	5.7
$c_{0,I}$	cm/s	700	$c_{0,I}$	cm/s	700
$L_I$	cm	4.4	$L_{II}$	cm	3
$R_1$	Ba·s/mL	60	$R_2$	Ba·s/mL	1500
C	mL/Ba	${10}^{-3}$	$p_{\infty }$	Ba	7000
$\rho$	g/cm3	1.04	$\mu$	cP	4

, $p_{WK}$ is pressure in the Windkessel compartment.

The formulations of boundary conditions at the aortic root (11), at the connection of the aorta and the pump (12), (13) and at the terminal part of the aorta (14)–(16) include a numerical discretisation of compatibility condition along the characteristic curve of the hyperbolic system (7) which leaves the integration domain for every incoming and/or outgoing segment of the aorta (see [21,22,23] for details). The systems of nonlinear algebraic equations, which represent boundary conditions with the time discretisation of the differential part are numerically solved by the Newton’s method.

The hyperbolic system (7) inside every segment is numerically solved by the second order grid-characteristic method (see [29] for the details of the method and [21,22] for the features of its implementation to the 1D model of the blood flow). The analysis of the characteristic curves of (7) and similar formulations of 1D blood flow model allows implementing discontinuous Galerkin method [30,31]. The deep analysis of the quasilinear effects in a hyperbolic model blood flow through compliant axi-symmetric vessels can be found in [32]. The generalised approach to the numerical implementation of the models describing various nonlinear wave process on graphs is described in [33].

The parameters of the 1D model are given in Table 3. The cross-sectional area and the length of the aortic segments I and II are set according to [34]. The parameters of the Windkessel compartment are set manually. These values allow us to achieve the well known systolic and diastolic aortal pressures in the normal conditions (rf. Section 3.1). For $\rho$ and $\mu$ we use the well known values [35].

2.3. Integrated Mathematical Model of the Heart Function, Pump and Aortic Flow

The two chamber model of the heart comprises the LV and the left atrium (LA), the mitral and aortic valves. It connects the pulmonary veins (PV) with the aorta. The nonlinear LVAD compartment connects the LV with the aorta (see Figure 1). The variable elasticity concept of the heart contractions [36,37] allows describing the heart chambers dynamics by the following lumped compartment model

$$I_k\frac{d^2{V}_k}{d{t}^2}+R_k{P}_k\frac{d{V}_k}{dt}+E_k(t)(V_k-V_k^0)+P_k^0=P_k,$$

(17)

where $k\in \{lv,la\}$, indices $lv$ and $la$ refer to the LV and the LA, respectively, $V_k(t)$ is the volume of the chamber, $V_k^0$ is the reference volume of the chamber, $P_k(t)$ is the pressure in the chamber, $P_k^0$ is the reference pressure in the chamber, $I_k$ is the inertia coefficient of the chamber, $R_k$ is the hydraulic resistance coefficient of the chamber, $E_k(t)$ is variable elasticity which is approximated by

$$E_k(t)=E_{k,d}+\frac{E_{k,s}-E_{k,d}}{2}{e}_k(t),$$

(18)

$E_{k,d}$ and $E_{k,s}$ are elasticity constants related to the end diastolic and end systolic states of chamber k (rf.

Table 4. Parameters of the lumped model of the left heart. ${}^{\star }$ Parameter is set manually.

Parameter	Unit	Value	Reference	Parameter	Unit	Value	Reference
$E_{lv,s}$	$\frac{\mathrm{mm}\mathrm{Hg}}{\mathrm{mL}}$	$4.0$	${}^{\star }$	${\theta }_{ao}^{min}$		$0^{°}$	[38]
$E_{lv,d}$	$\frac{\mathrm{mm}\mathrm{Hg}}{\mathrm{mL}}$	$0.09$	${}^{\star }$	${\theta }_{ao}^{max}$		${75}^{°}$	[43]
$I_{lv}$	$\frac{\mathrm{mm}\mathrm{Hg}·{\mathrm{s}}^2}{\mathrm{mL}}$	${10}^{-7}$	${}^{\star }$	${\theta }_{mi}^{min}$		$0^{°}$	[38]
$R_{lv}$	$\frac{\mathrm{s}}{\mathrm{mL}}$	$1.5\times {10}^{-3}$	[41]	${\theta }_{mi}^{max}$		${75}^{°}$	[43]
$E_{la,s}$	$\frac{\mathrm{mm}\mathrm{Hg}}{\mathrm{mL}}$	$1.2$	${}^{\star }$	$V_{lv}^0$	mL	5	${}^{\star }$
$E_{la,d}$	$\frac{\mathrm{mm}\mathrm{Hg}}{\mathrm{mL}}$	$0.3$	${}^{\star }$	$V_{la}^0$	mL	4	${}^{\star }$
$I_{la}$	$\frac{\mathrm{mm}\mathrm{Hg}·{\mathrm{s}}^2}{\mathrm{mL}}$	${10}^{-7}$	${}^{\star }$	$T_{s1}$	s	$0.3$	[43]
$R_{la}$	$\frac{\mathrm{s}}{\mathrm{mL}}$	$1.5\times {10}^{-3}$	[41]	$T_{s2}$	s	$0.35$	[43]
$T_{pw}$	s	0.1	[43]	$T_{pb}$	s	$0.9$	[43]
$K^p$	$\frac{\mathrm{rad}}{{\mathrm{s}}^2·\mathrm{mm}\mathrm{Hg}}$	${10}^4$	${}^{\star }$	$K^f$	$\frac{\mathrm{rad}}{\mathrm{s}}$	50	[43]
$P_{pv}$	mm Hg	13	${}^{\star }$	$S_p$	cm${}^2$	1.1	${}^{\star }$
$L_{pv}$	$\frac{\mathrm{mm}\mathrm{Hg}·{\mathrm{s}}^2}{\mathrm{mL}}$	${10}^{-2}$	${}^{\star }$	${\beta }_{pv}$	$\frac{\mathrm{mm}\mathrm{Hg}·{\mathrm{s}}^2}{{\mathrm{mL}}^2}$	$4\times {10}^{-4}$	${}^{\star }$
$L_{mi}$	$\frac{\mathrm{mm}\mathrm{Hg}·{\mathrm{s}}^2}{\mathrm{mL}}$	$5\times {10}^{-10}$	${}^{\star }$	$B_{mi}$		300	${}^{\star }$
$L_{ao}$	$\frac{\mathrm{mm}\mathrm{Hg}·{\mathrm{s}}^2}{\mathrm{mL}}$	$5\times {10}^{-5}$	${}^{\star }$	$B_{ao}$		500	${}^{\star }$
${\tilde{A}}_{mi}$	cm${}^2$	5	${}^{\star }$	$A_{ao}^{max}$	cm${}^2$	4	${}^{\star }$
T	s	1	${}^{\star }$	$A_{mi}^{max}$	cm${}^2$	4	${}^{\star }$

). For the LV we set

$$e_{lv}(t)=\{\begin{array}{cc}1-\cos \left(\frac{t}{T_{s1}}\pi \right),\hfill & 0⩽t⩽T_{s1},\hfill \\ 1+\cos \left(\frac{t-T_{s1}}{T_{s2}-T_{s1}}\pi \right),\hfill & T_{s1}<t<T_{s2},\hfill \\ 0,\hfill & T_{s2}⩽t⩽T,\hfill \end{array}$$

(19)

whereas for the LA

$$e_{la}(t)=\{\begin{array}{cc}0,\hfill & 0⩽t⩽T_{pb},\hfill \\ 1-\cos \left(\frac{t-T_{pb}}{T_{pw}}2\pi \right),\hfill & T_{pb}<t<T.\hfill \end{array}$$

(20)

Here we modify the model [38] by adding to (17) the term proportional to $P_k\frac{d{V}_k}{dt}$, which accounts for viscoelasticity of the myocardium [39,40,41]. The values of constants $T_{s1}$, $T_{s2}$, $T_{pb}$, $T_{pw}$ are presented in Table 4.

The mass conservation law for the LV and LA states

$$\begin{array}{c}\frac{d{V}_{lv}}{dt}=Q_{mi}-Q_{ao}-Q_p,\hfill \\ \frac{d{V}_{la}}{dt}=Q_{pv}-Q_{mi},\hfill \end{array}$$

(21)

where $Q_{mi}$ is the flow through the mitral valve, $Q_{ao}$ is the flow through the aortic valve, $Q_p$ is the flow through the pump, $Q_{pv}$ is the flow from the PV.

We set the pressure drop $\Delta P=P_{pv}-P_{la}$ for PV – LA connection, $\Delta P=P_{la}-P_{lv}$ for LA – LV connection, and $\Delta P=P_{lv}-p(A_I(t,0))$ for LV – AA connection. For unsteady flow in a channel with a variable cross-section, the pressure drop satisfies the relation [39,42]

$$\Delta P=L(g)\frac{dQ}{dt}+\alpha (g)Q+\beta (g)Q|Q|,$$

(22)

where $g(\theta )=\{{\theta }^{min}⩽\theta ⩽{\theta }^{max},0⩽g(\theta )⩽1\}$ is a smooth monotone function of the angle of a valve opening $\theta$ [43]:

$$g(\theta )=\{\begin{array}{cc}\frac{{(1-\cos {\theta }^{min})}^2}{{(1-\cos {\theta }^{max})}^2},\hfill & \theta <{\theta }^{min},\hfill \\ \frac{{(1-\cos \theta )}^2}{{(1-\cos {\theta }^{max})}^2},\hfill & {\theta }^{min}⩽\theta ⩽{\theta }^{max},\hfill \\ 1,\hfill & \theta >{\theta }^{max}.\hfill \end{array}$$

(23)

The value $g({\theta }^{min})$ corresponds to the closed valve, while the value $g({\theta }^{max})=1$ corresponds to the opened valve. For $L=0$, $\beta =0$ we have linear Poiseuille pressure drop condition which also accounts for the viscous friction losses. By analogy with [44,45] we neglect this term and set $\alpha =0$ for all cases. For $L=0$, $\alpha =0$ we have the orifice pressure drop condition. The first term in (22) accounts for the inertia of non-stationary flow. The coefficient $\beta$ is defined as [42,46,47]

$$\beta (A_{*})=\frac{\rho }{2B_{*}}{(\frac{1}{{\tilde{A}}_{*}}-\frac{1}{A_{*}})}^2,$$

(24)

where parameters ${\tilde{A}}_{mi}$, $B_{ao}$ and $B_{mi}$ are defined in Table 4 whereas ${\tilde{A}}_{ao}=A_I(t,0)$. For the PV – LA connection $\beta =const$. For both mitral and aortic valves, their cross-section $A_{*}$ depends on the angle of the valve opening, $A_{*}(\theta )=A_{*}^{max}g(\theta )$.

The dynamics of the aortic and mitral valves is governed by the second Newton law. Pressure gradient across the valve, vorticity generation and shear forces acting on the valve leaflets [48] have to be accounted in the model, rf. [43,49]. In this work we set valve dynamics equations as

$$\begin{array}{c}\frac{d^2{\theta }_{ao}}{d{t}^2}=-K_{ao}^f\frac{d{\theta }_{ao}}{dt}+(P_{lv}-P_{ao})K_{ao}^p\cos {\theta }_{ao}-F_{ao}^r({\theta }_{ao}),\hfill \\ \hfill \frac{d^2{\theta }_{mi}}{d{t}^2}=-K_{mi}^f\frac{d{\theta }_{mi}}{dt}+(P_{la}-P_{lv})K_{mi}^p\cos {\theta }_{mi}-F_{mi}^r({\theta }_{mi}).\end{array}$$

(25)

where ${\theta }_{ao}(t)$ is the angle of the aortic valve opening, ${\theta }_{mi}(t)$ is the angle of the mitral valve opening, $K_{ao}^f$, $K_{ao}^p$, $K_{mi}^f$, $K_{mi}^p$ are parameters presented in Table 4, the first term at the right-hand side corresponds to the friction force, the second term corresponds to the pressure force driving the valve motion, $F^r$ is the force which helps to avoid physiologically abnormal valve positions ($\theta <{\theta }^{min}$ and $\theta >{\theta }^{max}$)

$$F^r(\theta )=\{\begin{array}{cc}0,\hfill & {\theta }^{min}⩽\theta ⩽{\theta }^{max},\hfill \\ e^{{10}^3(\theta -{\theta }^{max})}-1,\hfill & \theta >{\theta }^{max},\hfill \\ 1-e^{{10}^3({\theta }^{min}-\theta )},\hfill & \theta <{\theta }^{min}.\hfill \end{array}$$

(26)

The other forces are neglected.

Parameters of the lumped model of the left heart are summarized in Table 4. We take some values from [38,41,43] and set other values manually basing on the values from [38,43] and keeping them in the physiological range.

3. Results

Validation of the pump model is addressed in Section 2.1. In Section 3.1, we validate the integrated model by comparing simulations of the normal heart function with known physiological data from the literature. In Section 3.2, we analyze the effect of four LVADs in patients with the end-stage HF associated with the LV dilated cardiomyopathy (DCM).

3.1. Validation of the Model

The parameters of the integrated model for healthy conditions are shown in Table 2, Table 3 and Table 4. Simulations with these parameters without LVAD produce values of stroke volume of LV, systolic and diastolic pressures in the aortic root which are in a good agreement with the well-known physiological data [35,50] (rf.

Table 5. Validation of the integrated model for the heart and the aorta in healthy conditions without LVAD.

Parameter	Unit	[35,50]	Model
End systolic LV volume	mL	50–70	42
End diastolic LV volume	mL	130	121
Stroke volume of LV	mL	60–80	79
Systolic pressure in the aortic root	mm Hg	130	124
Diastolic pressure in the aortic root	mm Hg	78	76

). We observe remarkable difference in the end systolic and end diastolic volumes of the LV. We note, that these parameters are highly individual. Even in healthy cases they depend on many factors including age, sex, sports lifestyle etc. [50,51]. Our values fall in the physiological range: they are typical for men of the age 70–80, women of the age about 40–70 [51] and other individual cases.

The PV diagram of the LV and time curves of the LV volume, the aortic flow and the aortic pressure in the terminal point of segment II (rf. Figure 1) in healthy conditions without a pump are shown in Figure 2. The loop in the lower right part of the PV diagram (rf. Figure 2a) accounts for the backflow from LV to LA through the mitral valve in early systole. This backflow is a result of non-instant closing of the mitral valve. Such loop is observed in critically ill patients. It is typical both for the right ventricle and the atria, but also it may be monitored in the left ventricle [52]. We also observe aortic regurgitation at the end of systole which accounts for the non-instant closing of the aortic valve [35,50] (rf. Figure 2c for $t\approx 0.3$ s).

3.2. Haemodynamic Simulations in the Aorta for LV DCM with LVAD

We simulate the haemodynamic characteristics in the left heart and the aorta in the presence of one of four LVADs (Sputnik 1, Sputnik 2, Sputnik D and HeartMate II) operating at various rotation speeds under HF conditions. RBPs are applied as long term circulatory support systems in patients with end-stage HF both as a bridge to heart transplantation and as an alternative to heart transplantation [53,54]. Some types of HF are accompanied with LV DCM which is the common indication for the long term LVAD installation. LV DCM is characterized by decreased LV contractility, thinning of the LV wall and increased cavity volume of the LV. These changes produce substantial decrease in the cardiac output and related cardiovascular dysfunction. LVAD unloads the LV and decreases its volume by pumping a portion of blood to the aorta. Thus, it supports the heart and sometimes may produce conditions for LV wall recovery. We update some parameters of the heart with LV DCM as shown in

Table 6. Parameters of the heart model with LV DCM.

Parameter	Unit	Value
$P_{pv}$	mm Hg	10
$E_{lv,d}$	mm Hg/mL	0.04
$E_{lv,s}$	mm Hg/mL	0.44
$V_{0,lv}$	mL	20
$E_{la,s}$	mm Hg/mL	1.1
$R_{lv}$	s/mL	$5\times {10}^{-4}$
$R_{la}$	s/mL	$5\times {10}^{-4}$

. The other parameters from Table 4 remain intact. The simulations with these parameters produce values which correlate with physiological data from the literature [53,54] (see

Table 7. Comparison of LV DCM simulations without LVAD with the data from the literature [53,54].

Parameter	Unit	[13,53,54]	Model
End systolic LV volume	mL	215	227
End diastolic LV volume	mL	259	275
Stroke volume of LV	mL	44	48
Systolic pressure in the aortic root	mm Hg	83	81
Diastolic pressure in the aortic root	mm Hg	55	47

for comparison).

We compare the function of the four pumps by setting the same parameters of the heart function and the aorta for all cases. The normal operating conditions of the pump are defined according to [13]. In the context of our model we formulate them as follows:

• The aortic valve should be opened within a part of cardiac cycle, i.e., the LV should eject some portion of blood to the aorta.
• The flow through the pump should be positive, i.e., it always should be directed from the LV to the aorta.
• The ventricular suction is not admitted.
• The total ejected volume per cardiac cycle in the aorta should be possibly close to the physiological value in normal conditions (see Section 3.1).

Figure 3 shows the duration of the opening of the aortic and mitral valves. Duration of the opening of the aortic and mitral valves at zero speed of the pumps are in agreement with the data on valve function in the normal conditions [43,49]. Decrease of time of the aortic valve opening starts at $5\times {10}^3$ rpm for all four pumps. The permanent closure of the aortic valve is observed at $7\times {10}^3$ rpm for Sputnik 1 and Sputnik 2, at $12\times {10}^3$ rpm for Sputnik D and at $8\times {10}^3$ rpm for HeartMate II which produce the upper bounds for the settings of these devices in the considered conditions. The noticeable increase of the mitral valve opening time starts at $9.5\times {10}^3$ rpm for all four pumps. The permanent opening of the mitral valve is observed beyond the values $12\times {10}^3$ rpm for Sputnik 1, $11.5\times {10}^3$ rpm for Sputnik 2, $13.5\times {10}^3$ rpm for HeartMate II and is not observed for Sputnik D in the whole range.

Figure 4 shows the minimum flow through the pump over the cardiac cycle. We observe that a negative flow through the pump disappears at $5.75\times {10}^3$ rpm for Sputnik 1 and Sputnik 2, at $8.5\times {10}^3$ rpm for Sputnik D and at $7\times {10}^3$ rpm for HeartMate II which sets the lower bound for the settings of these devices in the considered conditions.

Figure 5 shows the volume of the blood which is ejected through the aortic valve, through the pump and the total volume ejected to the ascending aorta per cardiac cycle. These volumes are calculated as the time integral of the corresponding flow. Thus, the negative value of the volume ejected through the pump means that the pump takes the blood from aorta back to the LV (see Figure 4). The total volume ejected to the ascending aorta per cardiac cycle is an analog of the stroke volume (SV). In the rest of the paper we refer to it as SV. We observe zero value of the volume ejected through the aortic valve for the values of the LVAD rotation speed which correspond to the permanent closure of the aortic valve (see Figure 3). The normal physiological value of the SV is observed at $8.5\times {10}^3$ rpm for Sputnik 1 and Sputnik 2, at $15\times {10}^3$ rpm for Sputnik D and at ${10}^4$ rpm for HeartMate II. This means that none of the four pumps can produce the normal SV in the considered conditions jointly with the opening of the aortic valve (see Figure 3).

It is interesting to observe (see Figure 5) that the volume of blood passing through the aortic valve is close to the normal value at zero rotation speed of the pump rotor. This effect is due to the negative flow through the pump from the aorta back to the LV (rf. Figure 4). From Figure 5 we notice that the pumps with speed below $5\times {10}^3$ rpm produce the same SV as LV DCM model without a pump. The effect of the pumps on the SV becomes significant only after decrease of the aortic valve opening time (rf. Figure 3).

Figure 6 shows the systolic and diastolic pressures in the aortic root. We observe that in a particular range (up to $7\times {10}^3$ rpm for Sputnik 1 and Sputnik 2, up to $12\times {10}^3$ rpm for Sputnik D and up to $8\times {10}^3$ rpm for HeartMate II ) the pump increases the diastolic pressure. The systolic pressure remains the same in this range. At the upper bound of this range, the systolic and diastolic pressures are almost equal, and the pulse pressure tends to zero. The elevation of the systolic pressure starts at the rotation speed corresponding to the permanent closure of the aortic valve (rf. Figure 3). The pulse pressure remains almost zero. The flow becomes non-pulsatile.

Figure 7 shows the work of the LV. It decreases with the increase of the pump speed and, thus, with the rise of the pump work. We observe three specific values of the pump speed. For values below $5\times {10}^3$ rpm, the LV work is almost constant due to the decreasing time of the aortic valve opening (rf. Figure 3). For values above $5\times {10}^3$ rpm, the LV work decreases. The kink of the work curve at $7\times {10}^3$ rpm for Sputnik 1 and Sputnik 2, $12\times {10}^3$ rpm for Sputnik D and $8\times {10}^3$ rpm for HeartMate II occurs due to the permanent closure of the aortic valve (rf. Figure 3). The LV performs almost zero work at $12\times {10}^3$ rpm for Sputnik 1 and Sputnik 2, $14\times {10}^3$ rpm for HeartMate II. Zero LV work is not observed for Sputnik D due to the permanent opening of the mitral valve (rf. Figure 3).

4. Discussion

In this work we study the impact of LVADs Sputnik 1, Sputnik 2, Sputnik D and HeartMate II on haemodynamics in the left heart and the aorta at various pump speeds. For every pump we validate the model of pressure – flow relationship using data from physical experiments. The pump model is combined with the 1D model of the aorta and the lumped model of the left heart with valves dynamics. The model without the pump reproduces successfully the known physiological characteristics of the heart in the healthy conditions: PV diagram of the LV, the LV volume, the aortic flow and the aortic pressure. We fit the model parameters to the LV DCM conditions and perform haemodynamic simulations with the pumps.

We observe different regimes of the heart function depending on the value of the pump speed. At low rotation speeds, the work of the pump is insufficient, and we observe reverse blood flow from the aorta to the LV through the pump. Formally, the volume of the blood ejected through the aortic valve is close to the standard value, but the pump backflow decreases the total amount of blood ejected to segment II of the aorta and makes it close to the value of LV DCM case without the pump. Increase of the rotation speed reduces the time of the aortic valve opening until valve’s permanent closure. Three of the four operating conditions of the normal pump functioning (see Section 3.2) hold within this range. Unfortunately, none of the pumps are capable of recovering the standard SV within this range. Further increase of the rotation speed produces non-pulsatile flow and permanent opening of the mitral valve. The technical characteristics and clinical restrictions limit the in vivo change of the rotation speed. For Sputnik 1 and Sputnik 2 the range is $5\times {10}^3$–${10}^4$ rpm, for Sputnik D the range is $6\times {10}^3$–$2\times {10}^4$ rpm, for HeartMate II the range is $6\times {10}^3$–$15\times {10}^3$ rpm.

In general, we observe the slight difference between Sputnik 1 and Sputnik 2 in all parameters (see Figure 3, Figure 4, Figure 5, Figure 6 and Figure 7). These pumps produce similar impacts on the haemodynamics, although they have different technical characteristics. Therefore, Sputnik 2 should be preferred as it provides such benefits as lesser weight, size, etc. [8,9,10]. The HeartMate II has a wider range of the rotation speed with acceptable operating conditions and provides flexibility in tuning the device settings. However, this LVAD recovers the normal SV at a higher rotation speed (see Figure 5a,b,d, which may require more energy. Sputnik D was initially designed for paediatric patients. It has lesser weight and size, but it also has less power. In our simulations, we test this LVAD in the adult LV DCM conditions. We show that it is possible to achieve heart operating conditions which are similar to the conditions with the adult LVADs. Sputnik D produces these conditions at a higher pump speed.

In this work we validate the model using the general characteristics of the heart function in the healthy and LV DCM conditions and data from the test facility. Personalized parameters may produce a more accurate assessment of LVAD impact.

The following parameters cause a substantial impact on the haemodynamics: elasticity of the aorta and the heart chambers, the central venous pressure, insufficiency of the heart valves, parameters of the peripheral circulation (Windkessel model). The dependency of LVAD function on these parameters should be studied in a future work.

In Section 3.2 we observe that the normal SV is achieved at zero pulse pressure and the closed aortic valve. Healthy blood flow in systemic arteries is a nonlinear wave phenomenon. The absence of pulsations may decrease the blood velocity and promote conditions for blood coagulation and thrombus formation in distal arteries. Thus, the loss of pressure and flow pulsatility in the systemic circulation in the presence of LVAD should be analyzed.

The variation of the pump speed changes the haemodynamic parameters in the heart and the aorta. These changes may activate regulatory mechanisms of the heart function. For instance, the change of the heart rate is associated with a modified duration of the systole and the value of the LV output through the aortic valve. These effects are beyond the scope of this work.

The changes in the SV and the aortic pressure modify the central venous pressure which, in turn, changes the LA filling conditions and adjusts the other haemodynamic parameters of the heart function. A closed model of the cardiovascular system is needed to study such feedback system.