# Stochastic Modeling and Dynamic Analysis of the Cardiovascular System with Rotary Left Ventricular Assist Devices

Authors: Jeongeun Son; Dongping Du; Yuncheng Du
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

Left ventricular assist devices (LVADs) have been used for end-stage heart failure patients as a therapeutic option. The aortic valve plays a critical role in heart failure and its treatment with a LVAD. The cardiovascular-LVAD model is often used to investigate the physiological demands required by patients and predict the hemodynamic of the native heart supported with a LVAD. As it is a “bridge-to-recovery” treatment, it is important to maintain appropriate and active dynamics of the aortic valve and the cardiac output of the native heart, which requires that the LVAD pump be adjusted so that a proper balance between the blood contributed through the aortic valve and the pump is maintained. In this paper, we investigate how the pump power of the LVAD pump can affect the dynamic behaviors of the aortic valve for different levels of activity and different severities of heart failure. Our objective is to identify a critical value of the pump power (i.e.,breakpoint) to ensure that the LVAD pump does not take over the pumping function in the cardiovascular-pump system and share the ejected blood with the left ventricle to help the heart to recover. In addition, the hemodynamic often involves variability due to patients’ heterogeneity and the stochastic nature of the cardiovascular system. The variability poses significant challenges to understanding dynamic behaviors of the aortic valve and cardiac output. A generalized polynomial chaos (gPC) expansion is used in this work to develop a stochastic cardiovascular-pump model for efficient uncertainty propagation, from which it is possible to rapidly calculate the variance in the aortic valve opening duration and the cardiac output in the presence of variability. The simulation results show that the gPC-based cardiovascular-pump model is a reliable platform that can provide useful information to understand the effect of the LVAD pump on the hemodynamic of the heart.


---

## Body

## 1. Introduction

Cardiovascular disease is one of the major causes of death in the United States. Approximately 5.7 million adults in the USA suffer from heart failure (HF). HF occurs when the heart fails to maintain appropriate circulation to support the physiological demands of the patient’s body [1]. Heart transplantation is the well-recognized treatment for end-stage HF. However, only a few patients are eligible for transplantation, due to the limited organ donors as well as the physical limitations such as age, health condition, or other health issues (i.e., impaired renal function, other comorbidities, or a high pulmonary vascular resistance) [2]. To overcome this limitation, an alternative treatment is to implant a ventricular assist device (VAD) to help unload the ventricles. VADs are mechanical pumps, which are designed to assist either the right ventricle or the left ventricle, or both ventricles in some cases, to eject the blood into the arterial system and further into the peripheral and end-organ [2].The left ventricular assist device (LVAD) is the most commonly used device for HF patients, since the right side of the heart can often make use of the heavily increased blood flow from the LVAD. The LVAD can partially replace the mechanical work of the failing left ventricle to maintain a desired blood flow between the left ventricle and the aorta. For example, it has been used to support an ailing heart as a “bridge to transplant” until a suitable donor heart is available. In addition, it is considered as a “destination therapy” for HF patients who are not eligible for heart transplantation. Recently, LVADs have been proposed as a “bridge-to-recovery” therapeutic option to help patients recover normal heart function [3, 4]. It was previously reported that the native heart function of patients can be improved with the support of LVADs. The reverse of HF can possibly allow patients to return to their normal life without the LVADs and potentially improve the quality of life for HF patients [5–7].LVADs can be generally divided into two generations by the pump types, i.e., pulsatile pumps and rotary pumps. Pulsatile pumps generate pulsatile blood flow close to the native heart in abeat-like fashion, whereas rotary pumps generate a continuous blood flow [8]. A rotary pump-based LVAD has more advantages over pulsatile pumps in terms of size, efficiency, durability, noise, and weight [9]. However, an important issue with a rotary LVAD is the optimal control of pump speed, which limits its extensive use. When the pump is operated at a lower rotational speed, it can induce regurgitation (i.e., backflow) from the aorta to the left ventricle [10]. In contrast, ventricle suction can happen when the pump speed is too high, which can lead to ventricular collapse, because the pump draws more blood from the left ventricle than available. In addition, the tuning of the pump speed can inevitably affect the function of the aortic valve. For example, inappropriate selection of pump speed may lead to permanent closure of the aortic valve, which is detrimental to cardiac recovery.The LVAD pump is normally set at a constant speed by physicians during the implantation surgery and cannot be adjusted freely. However, the heart function changes over time and patients may engage in a time-varying activity such as sleep, rest, and wild exercise. For both cases, the LVADs should be able to adjust the pump speed to meet different physiologic demands without inducing ventricular suction or regurgitation. The control of a LVAD is a difficult problem to formulate, since physiological variables of patients with the implanted LVADs have not been well studied and the effect of changes in control variables on the cardiovascular-pump system is not well understood either. Another challenge associated with LVAD control is to assess the aortic valve dynamics, while adjusting a control variable such as pump speed or pump power. To maintain normal operation of the aortic valve for the “bridge to recovery,” it is important to balance the amount of blood ejected through the pump and the aortic valve in order to avoid the situation in which the LVAD dominates blood circulation and takes over heart function.The aortic valve opens or closes periodically to allow the blood flows from the left ventricle to the aorta in each cardiac cycle. The aortic valve opens when the left ventricular pressure (LVP) is larger than the aortic pressure (AoP). As the blood flows out of the left ventricle, theLVP decreases and the aortic valve closes. However, when the pump speed is too high, the LVAD will provide the majority of left ventricle unloading; thus the left ventricle cannot generate a sufficient pressure to open the aortic valve [11]. Consequently, the LVAD will bypass the left ventricle and the aortic valve will be closed permanently. This can significantly change circulation physiology and introduce complications such as thrombosis and commissural fusion [12, 13]. The complications are fatal to HF patients, especially when the LVAD is used as a “bridge to recovery.” It is important to ensure that the aortic valve can remain active, when the pump speed of the LVAD is adjusted.Mathematical models of the human cardiovascular circulatory system have been developed to understand the cardiac hemodynamic. Most of these models integrate the left and right ventricles and atria with the systemic and pulmonary arterial and venous system to predict heart functions [14–17]. Although these models have the potential to predict dynamic changes of the physiological states, such as the left ventricular pressure, their clinical applications are still restricted, since these models fail to consider the interacting subsystems and networks in the cardiovascular system. In addition, cardiac function varies between different individuals and is different for the same patient over time. These variations, i.e., theinter- and/orintrapatient variability, defined as uncertainty, pose a major challenge to the development of an accurate model for the cardiovascular-pump system [18].To improve the reliability and credibility of the models, it is important to consider the uncertainty in the cardiovascular circulatory system. It should be noted that sampling-based techniques such as Monte Carlo (MC) simulations are the most popular method for uncertainty analysis [19]. However, MC-based methods can be computationally demanding for modeling and control of the cardiovascular circulatory system with an implanted LVAD, since a larger number of simulations are often required to obtain accurate results [20]. Recently, generalized polynomial chaos (gPC) expansion-based uncertainty quantification and propagation has been studied in different modeling, optimization, and control problems [21–25]. The gPC-based method can propagate a probabilistic uncertainty onto model predictions in a real-time manner, from which the uncertainty in model predictions can be easily estimated from gPC coefficients [22]. Due to the computational efficiency of the gPC, it is chosen for the uncertainty analysis in the cardiovascular-pump system in this work.Following the discussion above, a stochastic model of the human cardiovascular-pump system is developed in this work, using the gPC theory. In the presence of uncertainty, the dynamic behaviors of the aortic valve will be investigated for different electric powers. Note that the electric powers can be adjusted to vary the pump speed to meet various physiologic demands. The main contribution of this work is to efficiently quantify the uncertainty in cardiac outputs and dynamic behaviors of the aortic valve in each cardiac cycle of the cardiovascular-pump system. The uncertainty represents a time-varying physiologic change of patients in this work. Specifically, the aortic valve opening duration will be studied for different levels of physical activity and for different severities of HF. In addition, a probability description of the cardiac output such as the mean and the variance can be rapidly calculated using the gPC model, while taking uncertainty into account. It is important to note that the cardiac output can be determined by heart rate (HR) and stroke volume [10]. The stroke volume depends on the preload, contractility, and afterload of the heart. It is recognized that the rotary pump has poor sensitivity to the preload of the ventricle that is related to ventricular filling with venous blood and high sensitivity to the afterload of the ventricle that is the resistance to systolic ejection of blood [26, 27]. Therefore, we mainly focus on the analysis of the aortic valve dynamics and cardiac output in the presence of uncertainty in the afterload, i.e., systemic vascular resistance (SVR). The stochastic model can take into account the uncertainty among patients and individual patients’ physiological activities, which lays a firm foundation for control design of LVADs.This paper is organized as follows. Section2 presents a deterministic cardiovascular-LVAD model and the theoretical background of the gPC theory. The results of the computer simulations for the stochastic cardiovascular-pump model are presented in Section 3, which is followed by conclusions in Section 4.


## 2. Mathematical Background and Model

### 2.1. Cardiovascular-LVAD Model

The deterministic model of the cardiovascular-LVAD system in this work was experimentally validated by comparing the hemodynamic waveforms with data of patients [16, 28, 29]. It is assumed that the patient has a healthy and normal right ventricle and pulmonary system for simplicity. Thus, its effect on the LVAD is negligible [13, 16]. Figure 1 shows a schematic of the circuit model of the cardiovascular-LVAD system. The model parameters and their corresponding values are listed in Table 1.Table 1
Model parameters for the cardiovascular-LVAD system.
Circuit parameters
Value
Physiological Meaning
Resistance (mmHg∙s/ml)
R
s
1.0000
Systemic Vascular Resistance
R
M
0.0050
Mitral Valve Resistance
R
A
0.0010
Aortic Valve Resistance
R
c
0.0398
Characteristic Resistance
R
i
0.0677
Inlet Pump Resistance
R
o
0.0677
Outlet Pump Resistance
R
p
0.17070
Pump Resistance
R
k
See (1)
Suction Resistance
α = -3.5 s/ml and x
-
1 = 1 mmHg
Compliance (ml/mmHg)
C(t)
Time-varying
Left Ventricular Compliance
C
R
4.4000
Left Atrial Compliance
C
s
1.3300
Systemic Compliance
C
A
0.0800
Aortic Compliance
Inertances (mmHg∙s2/ml)
L
s
0.0005
Inertance of Blood in Aorta
L
i
0.0127
Inlet Inertance
L
o
0.0127
Outlet Inertance
L
p
0.02177
Pump Inertance
Valves (no units)
D
M
/
Mitral Valve
D
A
/
Aortic ValveFigure 1
Cardiovascular-LVAD circuit model.As seen in Figure1, the complianceC
R represents the preload and pulmonary circulation and the resistorsR
M andR
A are used to define the resistance related to the mitral valve and aortic valve, respectively.D
M andD
A are two ideal diodes to describe different phases in a cardiac cycle, i.e., (i) isovolumic relaxation, (ii) filling, (iii) isovolumic contraction, and (iv) ejection. For example, the mitral and aortic valves of a healthy heart are closed during isovolumetric relaxation or isovolumetric contraction.C
A is the aortic compliance, and parametersR
c,L
s,C
s, andR
s are used to describe the afterload in a four-element Windkessel model. It is important to note thatR
s, i.e., systemic vascular resistance (SVR), varies with respect to the level of physical activity of the patient. For example, when a person starts to do mild exercise from resting,R
s will decrease and the cardiac output will increase. Similarly, the value ofR
M, i.e., mitral valve resistance in Figure 1, will change with respect to the preload condition of the heart. The inlet and outlet resistances and inertances of the pump cannulae are described byR
i,R
o,L
i, andL
o, respectively. In addition, the pump resistance is defined asR
p, and the suction resistanceR
k, illustrating the phenomenon of suction, can be defined with two parameters as below:(1)
R
k
=
0
,
if
x
1
t
>
x
-
1
α
x
1
t
-
x
-
1
,
if
x
1
t
≤
x
-
1
,where α indicates a cannula dependent weighting parameter and x
-
1 is predetermined threshold pressure of the left ventricle.For the cardiovascular-pump system, the left ventricle complianceC(t) is a time-varying parameter, which is used to describe the contractibility of the left ventricle. In addition, the inverse ofC(t) is the elastic function of the left ventricle, i.e.,E(t), which can be described by the pressure and volume of the left ventricle as(2)
E
t
=
1
C
t
=
L
V
P
t
L
V
V
t
-
V
0
,whereLVP(t) andLVV(t) represent the pressure and volume of the left ventricle, respectively.V
0 is the theoretical volume at zero pressure defined as a reference volume of the ventricle. The mathematical expression of the elastance function E
t used in this work is described by a “double Hill” function E
n
t
n as below [16, 30]:(3)
E
t
=
E
m
a
x
-
E
m
i
n
E
n
t
n
+
E
m
i
n
(4)
E
n
t
n
=
1.55
t
n
/
0.7
1.9
1
+
t
n
/
0.7
1.9
1
1
+
t
n
/
1.17
21.9
.It is important to note that t
n is defined as follows: t
n =t/ (0.2+0.15t
c) in (4), and t
c represents the cardiac cycle, which is a sequence of events that occurs in each heart beat and is related to the heart rate (HR). Additionally, the elastance functionE(t) has different values corresponding to different heart conditions or severities of heart failure [13]. For example, the maximum value of the elastance function,E
max, can be set to 2 mmHg/ml, which means that the elastance and the compliance function can describe the dynamic behavior of a healthy heart. In contrast, if the value ofE
max is less than 2, it represents an unhealthy heart. Different values ofE
max can be used, depending on the severity of heart failure (HF). For clarity, Figure 2 shows the simulation results of the elastance function E
t for three values ofE
max, i.e., a healthy heart with a value of 2 mmHg/ml, mild heart failure (HF) with a value of 1 mmHg/ml, and a severely failing heart with a value of 0.5 mmHg/ml.Figure 2
Elastance function of the left ventricle with respect to different severities of heart condition (cardiac cycle = 60/HR, where HR is the heart rate).As mentioned above, the right ventricle and pulmonary system of the patient with the LVAD are assumed to be healthy in this work. Thus, the blood returning from the pulmonary system enters the left atrium and subsequently flows back into the left ventricle when the mitral valve opens. In addition, it is assumed that some portion of the blood in the left ventricle is ejected by the LVAD and the rest of the blood is pumped by the native heart through the aortic valve to the aorta and arterial system. In order to describe the cardiovascular-pump system, six state space variables, as shown in the circuit model in Figure1, are defined in Table 2.Table 2
State variables used in the cardiovascular-LVAD model.
Circuit Variables
Physiological variables
Physiological Meaning
Units
x
1(t)
LVP(t)
Left Ventricular Pressure
mmHg
x
2(t)
LAP(t)
Left Atrial Pressure
mmHg
x
3(t)
AP(t)
Arterial Pressure
mmHg
x
4(t)
AoP(t)
Aortic Pressure
mmHg
x
5(t)
Q
T(t)
Total Flow Rate
ml/s
x
6(t)
Qp(t)
Pump Flow Rate
ml/sBased on the basic circuit analysis, a6
t
h-order state space cardiovascular-pump model can be formulated to describe the hemodynamic of a HF patient implanted with a LVAD as(5)
x
˙
=
A
t
x
+
P
t
p
x
+
b
u
t
,wherex is a vector of the 6 state variables given in Table 2 and A
t and P
t are 6 × 6 and 6 × 2 time-varying matrices, respectively, which can be defined as(6)
A
t
=
-
C
˙
t
C
t
0
0
0
0
-
1
C
t
0
-
1
R
S
C
R
1
R
S
C
R
0
0
0
0
1
R
S
C
S
-
1
R
S
C
S
0
1
C
S
S
0
0
0
0
0
-
1
C
A
A
1
C
A
0
0
-
1
L
S
1
L
S
-
R
C
L
S
0
1
L
∗
0
0
-
1
L
∗
0
-
R
∗
L
∗
(7)
P
t
=
1
C
t
-
1
C
R
0
0
0
0
-
1
C
t
0
0
1
C
A
0
0
T
.It is important to note that a 2 × 1 vector p
x in (5) is used to mimic the nonlinear behaviors of the mitral and aortic valves, represented by two diodes in the circuit model in Figure 1. The details about this vector are given in (8) and (9). In addition, b in the last term in (5) is expressed by a 6 × 1 vector as given in (10). The control variable in (5) isu(t) = P
E(t), where P
E(t) is the electric pump motor power, which can be regulated to adjust the pump speed of the LVAD to meet the physiologic demands (e.g., blood) with respect to different levels of activity. The autonomous tuning of the pump speed is not discussed in this work, since our objective is to evaluate the distribution of the blood ejected through the LVAD and the aortic valve as a function of the electric power in the presence of uncertainty, e.g., variations in systemic vascular resistance (SVR).(8)
p
x
=
1
R
M
M
r
x
2
-
x
1
1
R
A
r
x
1
-
x
4
(9)
r
ζ
=
ξ
,
if
ζ
≥
0
0
,
if
ζ
<
0
(10)
b
x
=
0
0
0
0
0
δ
L
∗
x
6
T
.In addition, the pressure gain (or pump head) H
P (see Figure 1) can be defined as a function of the pressure difference across the pump and pump flow x
6 as shown in(11)
x
1
-
x
4
=
R
∗
x
6
+
L
∗
d
x
6
d
t
-
H
P
,where the pump head H
P can be approximated using the pump speed ω as follows [13]:(12)
H
P
=
β
ω
2
.Note that the coefficient β in (12) is set to 9.9025
×
10
-
7 mmHg/(rpm)2 [13]. In (11), the total resistance R
∗ and inductance L
∗ are two parameters related to the LVAD, which can be given by the following equations as(13)
R
∗
=
R
i
+
R
o
+
R
p
+
R
k
(14)
L
∗
=
L
i
+
L
o
+
L
p
.The pressure gain across the pump H
P has the direct relation to the electric power delivered to the pump motor P
E. Defining the pump efficiency as η, the electric power P
E can be related to the hydrodynamic power P
P as shown in (15). Further, P
P can be approximated with the density of the reference fluid ρ, gravitational acceleration g, the pump flow rate x
6, and the pump head H
P, as shown in (16).(15)
P
P
=
η
P
E
(16)
P
P
=
ρ
g
H
P
x
6
.Rearranging both equations above, the relation between the pressure gain H
P and the electric power P
E can be defined as follows [31]:(17)
H
p
=
δ
P
E
x
6
,where the constant δ is set to 7495 mmHg∙ml/s∙W, which is a function of the density of the reference fluid ρ, the gravity acceleration g, and the pump efficiency η [26]. For simplicity, η is set to 100% in this work.Based on the definition of pump headH
p in (12) and (17), the corresponding rotational pump speed ω can be defined in terms of the pump electric power P
E as(18)
ω
=
δ
P
E
β
x
6
.To study the effect of uncertainty in aortic valve dynamics and cardiac outputs, the operating power range of the LVAD pump is set to 0.12 to 1.56W in this work [32]. The corresponding rotational pump speeds for different electric powers are shown in Appendix A.


### 2.2. Generalized Polynomial Chaos (gPC) Expansion

The generalized polynomial chaos (gPC) expansion generally approximates a random variable using an arbitrary probability density function (PDF) defined by another random variable (e.g.,ξ) witha prior distribution (PDF) in the Wiener-Askey framework [22]. In this work, the gPC expansion will be used to approximate the uncertainty such as systemic vascular resistance (SVR) (i.e.,R
s) to study and quantify how uncertainty can affect the model predictions, e.g., the aortic valve dynamics and cardiac outputs. Note that SVR in the model is used to describe the different levels of activity, which can vary over time within the same HF patient. Thus, it is assumed that the exact value of SVR is unknown at each time; however, the PDF of SVR over a period of time is available, which can be determined by physicians or estimated through offline estimation techniques. The rationale of choosing SVR as the uncertain source in this work will be explained in the results session through the sensitivity analysis.Suppose that the cardiovascular-LVAD system in (5), described by a set of nonlinear ordinary differential equations (ODEs), can be simplified as follows:(19)
x
˙
=
f
t
,
x
,
v
,
g
;
u
,where the vector x consists of the 6 state variables defined in Table 2 with initial values x
0 att= 0. v is a vector of deterministic model parameters in the cardiovascular-pump system, which are fixed constants. In contrast, a vector of parametric uncertainties is defined by g, which will be approximated with their PDFs instead of using fixed constants, such as the SVR explained above, representing the level of activity of HF patients. In addition, u is the control variable, e.g., electric power P
E
(
t
), which can be adjusted to meet the physiologic demands.To evaluate the effect of uncertainty on the model predictionsx, each parameter g
i (i=1,2, …,n
g) in g will be approximated with a gPC model as a function of a set of independent random variables ξ
=
ξ
i as(20)
g
i
=
g
i
ξ
i
,where ξ
i is the i
t
h random variable used to approximate g
i that followsa prior PDF defined in the Wiener-Askey framework. Since each parametric uncertainty is approximated with a gPC model, the model predictions x can be also defined with random variables ξ, which can assess the effect of uncertainty on the model predictions of the cardiovascular-pump system. Using the orthogonal polynomial basis functions defined in the Wiener-Askey framework, the gPC approximations for both uncertainty and model predictions can be defined as follows:(21)
g
i
ξ
=
∑
k
=
0
∞
g
^
i
,
k
φ
k
ξ
i
≈
∑
k
=
0
q
g
^
i
,
k
φ
k
ξ
i
(22)
x
j
t
,
ξ
=
∑
k
=
0
∞
x
^
j
,
k
t
ψ
k
ξ
≈
∑
k
=
0
Q
x
^
j
,
k
t
ψ
k
ξ
,where g
^
i
,
k and x
^
j
,
k
(
t
) are the gPC coefficients of thei
th parametric uncertainty and thej
th model prediction (state variable) and φ
k
(
ξ
i
) and ψ
k
(
ξ
) are multidimensional polynomial basis functions. When the PDF of g
i is given or can be estimated so that g
i
ξ followsa prior distribution, the gPC coefficients g
^
i
,
k in (21) can be subsequently determined. As compared to g
^
i
,
k, the gPC coefficients of x
^
j
,
k
(
t
) will be calculated by substituting (21) and (22) into (5), which is followed by using a Galerkin projection. The Galerkin projection between a state variable x
j and a polynomial basis function can be defined as(23)
x
j
˙
t
,
ξ
,
ψ
k
ξ
=
f
t
,
x
j
t
,
ξ
,
v
,
g
ξ
;
u
,
ψ
k
ξ
.For each state variable x
j, the Galerkin projection will produce a system with a set of coupled deterministic equations, in which the first equation provides the mean value of x
j at each time instant, while the rest of the equations can be used to estimate the variation resulting from uncertainty. In addition, as seen in (21) and (22), the infinite number of terms is often truncated into a finite number of terms for practical application, i.e.,q andQ, respectively. The number of terms in (21), i.e.,q, can be optimally selected, so that the gPC model of uncertainty can approximatea prior known PDF of uncertainty. The total number of termsQ in (22) can be computed with a heuristic formula, which can be defined by the polynomial orderq and the total number of parametric uncertainties n
g as(24)
Q
=
n
g
+
q
!
n
g
!
q
!
–
1
.In addition, the inner product between any two vectors in (23) can be determined as(25)
ϕ
ξ
,
ϕ
′
ξ
=
∫
ϕ
ξ
ϕ
′
ξ
W
ξ
d
ξ
,where the integration on the right-hand side of (25) is performed over the entire domain defined by random variables ξ and W
(
ξ
) is the weighting function, i.e., the probability density function of ξ, which is selected according to the polynomial basis function in the Wiener-Askey framework. For example, a normally distributed random variable ξ should be used, when the uncertainty follows a normal distribution [22]. Thus, Hermite polynomial basis functions are the best choice of the weighting function.Using the gPC coefficients of state variablesx in (22), the statistical moments such as mean and variance of x at a given time intervalt can be quickly computed as follows:(26)
E
x
j
t
=
E
∑
i
=
0
Q
x
^
j
,
i
t
ψ
i
=
x
^
j
,
i
t
E
ψ
i
+
∑
i
=
1
Q
E
ψ
i
=
x
^
j
,
0
t
(27)
var
⁡
x
j
t
=
E
x
t
-
E
x
j
t
2
=
E
∑
i
=
0
Q
x
^
j
,
i
t
ψ
i
-
x
^
j
,
i
=
0
t
2
=
E
∑
i
=
1
Q
x
^
j
,
i
t
ψ
i
2
=
∑
i
=
1
Q
x
^
j
,
i
t
2
E
ψ
i
2
.As seen, the first statistical moment of x, i.e., mean value, can be approximated with the first gPC coefficient x
^
j
,
i
=
0, while the second statistical moment, i.e., variance of x, can be calculated using the higher-order gPC coefficients and the mean of squared orthogonal polynomial basis function ψ
i
(
ξ
). Once again, the variance in state variables x originates from uncertainty g, as defined in (20).The gPC provides analytical formulas to calculate the statistical moments of model predictions, from which the PDF profiles can be rapidly estimated. This is the main rationale to use the gPC approximation in this work. The fast calculation of uncertainty in hemodynamic waveforms, resulting from the time-varying physiologic change such as SVR (R
s) in this work, can provide useful information to better understand the physiologic demands. Specifically, the PDF profiles of x, approximated from a gPC model, will be used to account for the effect of uncertainty inR
s on the aortic valve opening duration and the cardiac outputs, which will be discussed in detail in Section 3.


## 2.1. Cardiovascular-LVAD Model

The deterministic model of the cardiovascular-LVAD system in this work was experimentally validated by comparing the hemodynamic waveforms with data of patients [16, 28, 29]. It is assumed that the patient has a healthy and normal right ventricle and pulmonary system for simplicity. Thus, its effect on the LVAD is negligible [13, 16]. Figure 1 shows a schematic of the circuit model of the cardiovascular-LVAD system. The model parameters and their corresponding values are listed in Table 1.Table 1
Model parameters for the cardiovascular-LVAD system.
Circuit parameters
Value
Physiological Meaning
Resistance (mmHg∙s/ml)
R
s
1.0000
Systemic Vascular Resistance
R
M
0.0050
Mitral Valve Resistance
R
A
0.0010
Aortic Valve Resistance
R
c
0.0398
Characteristic Resistance
R
i
0.0677
Inlet Pump Resistance
R
o
0.0677
Outlet Pump Resistance
R
p
0.17070
Pump Resistance
R
k
See (1)
Suction Resistance
α = -3.5 s/ml and x
-
1 = 1 mmHg
Compliance (ml/mmHg)
C(t)
Time-varying
Left Ventricular Compliance
C
R
4.4000
Left Atrial Compliance
C
s
1.3300
Systemic Compliance
C
A
0.0800
Aortic Compliance
Inertances (mmHg∙s2/ml)
L
s
0.0005
Inertance of Blood in Aorta
L
i
0.0127
Inlet Inertance
L
o
0.0127
Outlet Inertance
L
p
0.02177
Pump Inertance
Valves (no units)
D
M
/
Mitral Valve
D
A
/
Aortic ValveFigure 1
Cardiovascular-LVAD circuit model.As seen in Figure1, the complianceC
R represents the preload and pulmonary circulation and the resistorsR
M andR
A are used to define the resistance related to the mitral valve and aortic valve, respectively.D
M andD
A are two ideal diodes to describe different phases in a cardiac cycle, i.e., (i) isovolumic relaxation, (ii) filling, (iii) isovolumic contraction, and (iv) ejection. For example, the mitral and aortic valves of a healthy heart are closed during isovolumetric relaxation or isovolumetric contraction.C
A is the aortic compliance, and parametersR
c,L
s,C
s, andR
s are used to describe the afterload in a four-element Windkessel model. It is important to note thatR
s, i.e., systemic vascular resistance (SVR), varies with respect to the level of physical activity of the patient. For example, when a person starts to do mild exercise from resting,R
s will decrease and the cardiac output will increase. Similarly, the value ofR
M, i.e., mitral valve resistance in Figure 1, will change with respect to the preload condition of the heart. The inlet and outlet resistances and inertances of the pump cannulae are described byR
i,R
o,L
i, andL
o, respectively. In addition, the pump resistance is defined asR
p, and the suction resistanceR
k, illustrating the phenomenon of suction, can be defined with two parameters as below:(1)
R
k
=
0
,
if
x
1
t
>
x
-
1
α
x
1
t
-
x
-
1
,
if
x
1
t
≤
x
-
1
,where α indicates a cannula dependent weighting parameter and x
-
1 is predetermined threshold pressure of the left ventricle.For the cardiovascular-pump system, the left ventricle complianceC(t) is a time-varying parameter, which is used to describe the contractibility of the left ventricle. In addition, the inverse ofC(t) is the elastic function of the left ventricle, i.e.,E(t), which can be described by the pressure and volume of the left ventricle as(2)
E
t
=
1
C
t
=
L
V
P
t
L
V
V
t
-
V
0
,whereLVP(t) andLVV(t) represent the pressure and volume of the left ventricle, respectively.V
0 is the theoretical volume at zero pressure defined as a reference volume of the ventricle. The mathematical expression of the elastance function E
t used in this work is described by a “double Hill” function E
n
t
n as below [16, 30]:(3)
E
t
=
E
m
a
x
-
E
m
i
n
E
n
t
n
+
E
m
i
n
(4)
E
n
t
n
=
1.55
t
n
/
0.7
1.9
1
+
t
n
/
0.7
1.9
1
1
+
t
n
/
1.17
21.9
.It is important to note that t
n is defined as follows: t
n =t/ (0.2+0.15t
c) in (4), and t
c represents the cardiac cycle, which is a sequence of events that occurs in each heart beat and is related to the heart rate (HR). Additionally, the elastance functionE(t) has different values corresponding to different heart conditions or severities of heart failure [13]. For example, the maximum value of the elastance function,E
max, can be set to 2 mmHg/ml, which means that the elastance and the compliance function can describe the dynamic behavior of a healthy heart. In contrast, if the value ofE
max is less than 2, it represents an unhealthy heart. Different values ofE
max can be used, depending on the severity of heart failure (HF). For clarity, Figure 2 shows the simulation results of the elastance function E
t for three values ofE
max, i.e., a healthy heart with a value of 2 mmHg/ml, mild heart failure (HF) with a value of 1 mmHg/ml, and a severely failing heart with a value of 0.5 mmHg/ml.Figure 2
Elastance function of the left ventricle with respect to different severities of heart condition (cardiac cycle = 60/HR, where HR is the heart rate).As mentioned above, the right ventricle and pulmonary system of the patient with the LVAD are assumed to be healthy in this work. Thus, the blood returning from the pulmonary system enters the left atrium and subsequently flows back into the left ventricle when the mitral valve opens. In addition, it is assumed that some portion of the blood in the left ventricle is ejected by the LVAD and the rest of the blood is pumped by the native heart through the aortic valve to the aorta and arterial system. In order to describe the cardiovascular-pump system, six state space variables, as shown in the circuit model in Figure1, are defined in Table 2.Table 2
State variables used in the cardiovascular-LVAD model.
Circuit Variables
Physiological variables
Physiological Meaning
Units
x
1(t)
LVP(t)
Left Ventricular Pressure
mmHg
x
2(t)
LAP(t)
Left Atrial Pressure
mmHg
x
3(t)
AP(t)
Arterial Pressure
mmHg
x
4(t)
AoP(t)
Aortic Pressure
mmHg
x
5(t)
Q
T(t)
Total Flow Rate
ml/s
x
6(t)
Qp(t)
Pump Flow Rate
ml/sBased on the basic circuit analysis, a6
t
h-order state space cardiovascular-pump model can be formulated to describe the hemodynamic of a HF patient implanted with a LVAD as(5)
x
˙
=
A
t
x
+
P
t
p
x
+
b
u
t
,wherex is a vector of the 6 state variables given in Table 2 and A
t and P
t are 6 × 6 and 6 × 2 time-varying matrices, respectively, which can be defined as(6)
A
t
=
-
C
˙
t
C
t
0
0
0
0
-
1
C
t
0
-
1
R
S
C
R
1
R
S
C
R
0
0
0
0
1
R
S
C
S
-
1
R
S
C
S
0
1
C
S
S
0
0
0
0
0
-
1
C
A
A
1
C
A
0
0
-
1
L
S
1
L
S
-
R
C
L
S
0
1
L
∗
0
0
-
1
L
∗
0
-
R
∗
L
∗
(7)
P
t
=
1
C
t
-
1
C
R
0
0
0
0
-
1
C
t
0
0
1
C
A
0
0
T
.It is important to note that a 2 × 1 vector p
x in (5) is used to mimic the nonlinear behaviors of the mitral and aortic valves, represented by two diodes in the circuit model in Figure 1. The details about this vector are given in (8) and (9). In addition, b in the last term in (5) is expressed by a 6 × 1 vector as given in (10). The control variable in (5) isu(t) = P
E(t), where P
E(t) is the electric pump motor power, which can be regulated to adjust the pump speed of the LVAD to meet the physiologic demands (e.g., blood) with respect to different levels of activity. The autonomous tuning of the pump speed is not discussed in this work, since our objective is to evaluate the distribution of the blood ejected through the LVAD and the aortic valve as a function of the electric power in the presence of uncertainty, e.g., variations in systemic vascular resistance (SVR).(8)
p
x
=
1
R
M
M
r
x
2
-
x
1
1
R
A
r
x
1
-
x
4
(9)
r
ζ
=
ξ
,
if
ζ
≥
0
0
,
if
ζ
<
0
(10)
b
x
=
0
0
0
0
0
δ
L
∗
x
6
T
.In addition, the pressure gain (or pump head) H
P (see Figure 1) can be defined as a function of the pressure difference across the pump and pump flow x
6 as shown in(11)
x
1
-
x
4
=
R
∗
x
6
+
L
∗
d
x
6
d
t
-
H
P
,where the pump head H
P can be approximated using the pump speed ω as follows [13]:(12)
H
P
=
β
ω
2
.Note that the coefficient β in (12) is set to 9.9025
×
10
-
7 mmHg/(rpm)2 [13]. In (11), the total resistance R
∗ and inductance L
∗ are two parameters related to the LVAD, which can be given by the following equations as(13)
R
∗
=
R
i
+
R
o
+
R
p
+
R
k
(14)
L
∗
=
L
i
+
L
o
+
L
p
.The pressure gain across the pump H
P has the direct relation to the electric power delivered to the pump motor P
E. Defining the pump efficiency as η, the electric power P
E can be related to the hydrodynamic power P
P as shown in (15). Further, P
P can be approximated with the density of the reference fluid ρ, gravitational acceleration g, the pump flow rate x
6, and the pump head H
P, as shown in (16).(15)
P
P
=
η
P
E
(16)
P
P
=
ρ
g
H
P
x
6
.Rearranging both equations above, the relation between the pressure gain H
P and the electric power P
E can be defined as follows [31]:(17)
H
p
=
δ
P
E
x
6
,where the constant δ is set to 7495 mmHg∙ml/s∙W, which is a function of the density of the reference fluid ρ, the gravity acceleration g, and the pump efficiency η [26]. For simplicity, η is set to 100% in this work.Based on the definition of pump headH
p in (12) and (17), the corresponding rotational pump speed ω can be defined in terms of the pump electric power P
E as(18)
ω
=
δ
P
E
β
x
6
.To study the effect of uncertainty in aortic valve dynamics and cardiac outputs, the operating power range of the LVAD pump is set to 0.12 to 1.56W in this work [32]. The corresponding rotational pump speeds for different electric powers are shown in Appendix A.


## 2.2. Generalized Polynomial Chaos (gPC) Expansion

The generalized polynomial chaos (gPC) expansion generally approximates a random variable using an arbitrary probability density function (PDF) defined by another random variable (e.g.,ξ) witha prior distribution (PDF) in the Wiener-Askey framework [22]. In this work, the gPC expansion will be used to approximate the uncertainty such as systemic vascular resistance (SVR) (i.e.,R
s) to study and quantify how uncertainty can affect the model predictions, e.g., the aortic valve dynamics and cardiac outputs. Note that SVR in the model is used to describe the different levels of activity, which can vary over time within the same HF patient. Thus, it is assumed that the exact value of SVR is unknown at each time; however, the PDF of SVR over a period of time is available, which can be determined by physicians or estimated through offline estimation techniques. The rationale of choosing SVR as the uncertain source in this work will be explained in the results session through the sensitivity analysis.Suppose that the cardiovascular-LVAD system in (5), described by a set of nonlinear ordinary differential equations (ODEs), can be simplified as follows:(19)
x
˙
=
f
t
,
x
,
v
,
g
;
u
,where the vector x consists of the 6 state variables defined in Table 2 with initial values x
0 att= 0. v is a vector of deterministic model parameters in the cardiovascular-pump system, which are fixed constants. In contrast, a vector of parametric uncertainties is defined by g, which will be approximated with their PDFs instead of using fixed constants, such as the SVR explained above, representing the level of activity of HF patients. In addition, u is the control variable, e.g., electric power P
E
(
t
), which can be adjusted to meet the physiologic demands.To evaluate the effect of uncertainty on the model predictionsx, each parameter g
i (i=1,2, …,n
g) in g will be approximated with a gPC model as a function of a set of independent random variables ξ
=
ξ
i as(20)
g
i
=
g
i
ξ
i
,where ξ
i is the i
t
h random variable used to approximate g
i that followsa prior PDF defined in the Wiener-Askey framework. Since each parametric uncertainty is approximated with a gPC model, the model predictions x can be also defined with random variables ξ, which can assess the effect of uncertainty on the model predictions of the cardiovascular-pump system. Using the orthogonal polynomial basis functions defined in the Wiener-Askey framework, the gPC approximations for both uncertainty and model predictions can be defined as follows:(21)
g
i
ξ
=
∑
k
=
0
∞
g
^
i
,
k
φ
k
ξ
i
≈
∑
k
=
0
q
g
^
i
,
k
φ
k
ξ
i
(22)
x
j
t
,
ξ
=
∑
k
=
0
∞
x
^
j
,
k
t
ψ
k
ξ
≈
∑
k
=
0
Q
x
^
j
,
k
t
ψ
k
ξ
,where g
^
i
,
k and x
^
j
,
k
(
t
) are the gPC coefficients of thei
th parametric uncertainty and thej
th model prediction (state variable) and φ
k
(
ξ
i
) and ψ
k
(
ξ
) are multidimensional polynomial basis functions. When the PDF of g
i is given or can be estimated so that g
i
ξ followsa prior distribution, the gPC coefficients g
^
i
,
k in (21) can be subsequently determined. As compared to g
^
i
,
k, the gPC coefficients of x
^
j
,
k
(
t
) will be calculated by substituting (21) and (22) into (5), which is followed by using a Galerkin projection. The Galerkin projection between a state variable x
j and a polynomial basis function can be defined as(23)
x
j
˙
t
,
ξ
,
ψ
k
ξ
=
f
t
,
x
j
t
,
ξ
,
v
,
g
ξ
;
u
,
ψ
k
ξ
.For each state variable x
j, the Galerkin projection will produce a system with a set of coupled deterministic equations, in which the first equation provides the mean value of x
j at each time instant, while the rest of the equations can be used to estimate the variation resulting from uncertainty. In addition, as seen in (21) and (22), the infinite number of terms is often truncated into a finite number of terms for practical application, i.e.,q andQ, respectively. The number of terms in (21), i.e.,q, can be optimally selected, so that the gPC model of uncertainty can approximatea prior known PDF of uncertainty. The total number of termsQ in (22) can be computed with a heuristic formula, which can be defined by the polynomial orderq and the total number of parametric uncertainties n
g as(24)
Q
=
n
g
+
q
!
n
g
!
q
!
–
1
.In addition, the inner product between any two vectors in (23) can be determined as(25)
ϕ
ξ
,
ϕ
′
ξ
=
∫
ϕ
ξ
ϕ
′
ξ
W
ξ
d
ξ
,where the integration on the right-hand side of (25) is performed over the entire domain defined by random variables ξ and W
(
ξ
) is the weighting function, i.e., the probability density function of ξ, which is selected according to the polynomial basis function in the Wiener-Askey framework. For example, a normally distributed random variable ξ should be used, when the uncertainty follows a normal distribution [22]. Thus, Hermite polynomial basis functions are the best choice of the weighting function.Using the gPC coefficients of state variablesx in (22), the statistical moments such as mean and variance of x at a given time intervalt can be quickly computed as follows:(26)
E
x
j
t
=
E
∑
i
=
0
Q
x
^
j
,
i
t
ψ
i
=
x
^
j
,
i
t
E
ψ
i
+
∑
i
=
1
Q
E
ψ
i
=
x
^
j
,
0
t
(27)
var
⁡
x
j
t
=
E
x
t
-
E
x
j
t
2
=
E
∑
i
=
0
Q
x
^
j
,
i
t
ψ
i
-
x
^
j
,
i
=
0
t
2
=
E
∑
i
=
1
Q
x
^
j
,
i
t
ψ
i
2
=
∑
i
=
1
Q
x
^
j
,
i
t
2
E
ψ
i
2
.As seen, the first statistical moment of x, i.e., mean value, can be approximated with the first gPC coefficient x
^
j
,
i
=
0, while the second statistical moment, i.e., variance of x, can be calculated using the higher-order gPC coefficients and the mean of squared orthogonal polynomial basis function ψ
i
(
ξ
). Once again, the variance in state variables x originates from uncertainty g, as defined in (20).The gPC provides analytical formulas to calculate the statistical moments of model predictions, from which the PDF profiles can be rapidly estimated. This is the main rationale to use the gPC approximation in this work. The fast calculation of uncertainty in hemodynamic waveforms, resulting from the time-varying physiologic change such as SVR (R
s) in this work, can provide useful information to better understand the physiologic demands. Specifically, the PDF profiles of x, approximated from a gPC model, will be used to account for the effect of uncertainty inR
s on the aortic valve opening duration and the cardiac outputs, which will be discussed in detail in Section 3.


## 3. Results and Discussion

### 3.1. Cardiac Hemodynamic

The cardiovascular-pump model as explained in Section2.1 involves a few model parameters, which can be used to represent the physiological dynamics of the heart supported with a LVAD. As previously reported, parametersR
s,R
M, and heart rate (HR) can govern the behavior of the cardiovascular-pump system [13]. Thus, we propose to investigate their effect on the dynamics of the aortic valve and the cardiac output.Since the gPC-based stochastic model in this work is developed based on a deterministic model in (5), the first step is to validate the ability of the deterministic model in order to mimic the hemodynamic of the cardiovascular-pump system. The linear relationship between the end-systolic pressure and the left ventricle volume in the presence of perturbations inR
s andR
M is used. Note that the perturbations inR
s andR
M represent the stochastic changes during the preload and afterload of the heart. A total of 4 preload and 4 afterload changes were simulated as shown in Figure 3, for which the left ventricle parameters such as E
m
a
x, E
m
i
n, andV
0 are set to constant and the electric pump power is set to P
E = 0.12W.Figure 3
PV-loops for different values of (a) the resistanceR
s (afterload) and (b)R
M (preload).
(a)
(b)The left ventricle pressure and the left ventricle volume graph (PV-loop) are shown in Figure3. The PV-loops in Figure 3(a) show the hemodynamic changes resulting from variations in the afterload, i.e., different values ofR
s (SVR) that represent the changes in the level of activity. In contrast, the PV-loop in Figure 3(b) shows the result by altering the preload, i.e.,R
M representing mitral valve resistance. For both graphs, a linear relationship between the left ventricle pressure and the left ventricle volume is observed, thus confirming the ability of the model to describe the hemodynamic of the cardiovascular-pump system. In addition, it was found that any variations in the cardiovascular-pump system can significantly affect the PV-loop such as the area that represents the stroke work.Note that the dynamic behaviors of the aortic valve will be used in this work as a criterion to evaluate the effect of uncertainty on the cardiovascular-pump system. The rationale is that any changes in physiological parameters such asR
s can affect left ventricular pressure as seen in Figure 3(a). Constant low ventricular pressure can result in the permanent closure of the aortic valve, which is detrimental to cardiac recovery. Specifically, the aortic valve opening duration is used to quantitatively assess the effect of uncertainty on the dynamic behaviors of the aortic valve.The calculation of the aortic valve opening duration proceeds as follows. In a cardiac cycle, the period of time,t
open, during which the left ventricular pressure (LVP) is larger than the aortic pressure (AoP) will be first determined. The aortic valve opening duration can then be obtained by calculating the ratio betweent
open and the time of a cardiac cyclet
c. A schematic of the calculation of the aortic valve opening duration is given in Figure 4. Let us suppose there are 9 data points of bothLVP andAoP in a cardiac cycle, and the time interval between any two data points ist
int. As seen, two data points ofLVP are larger than theAoP, indicating thatt
open is 2×t
int. Similarly, the cardiac cycle can be calculated ast
c = 9×t
int. Thus, the aortic valve opening duration can be determined with 2/9 ≈ 0.222, representing that the aortic valve remains open during approximately 22.2% of a cardiac cycle.Figure 4
Schematic for the calculation of the aortic valve opening duration.The cardiac output (CO) is another key physical property, which can be calculated as a product of stroke volume and HR. The cardiac output represents the total amount of blood pumped from the ventricle in a minute. The stroke volume is the blood volume pumped from the left ventricle in a beat and is determined by the contractility, preload, and afterload. The area of the PV-loop is often calculated to represent the stroke work. As shown in Figure3, any changes in R
s andR
M can affect the total amount of blood entering the ventricle, thus affecting the stroke volume and the cardiac output. Therefore, the cardiac output is used as a second criterion to evaluate the effect of uncertainty on the cardiovascular-pump system in this work.


### 3.2. Sensitivity Analysis

As mentioned above,R
s,R
M, and HR can affect the physiological dynamics of the heart such as the aortic valve opening and the cardiac output. However, each factor may have different effect on the hemodynamic of the failing heart. Thus, a sensitivity analysis is first performed to identify the most significant physiological factor. The effect of variations inR
s,R
M, and HR on the aortic valve and the cardiac output is first investigated, using the cardiovascular-LVAD model with a mild HF patient (i.e.,E
max = 1.0) and electric pump power P
E = 0.12 ~ 0.6W. Similar results were found for other heart conditions and pump powers, but they are not shown for brevity.It is assumed that each parameter can vary randomly between two levels, i.e., +1 and -1, which correspond to a +10% change and a -10% change with respect to its nominal values. Note that the nominal values of each parameter areR
s = 1.0 mmHg∙s/ml,R
M = 0.005 mmHg∙s/ml, and HR = 75 bpm. For each parameter, let us define the aortic valve opening duration (or the cardiac output) as w
g
i
+ and w
g
i
- for two different levels and the corresponding result with the nominal values of model parameters as w
g
i
0. Then, the effect of uncertainty in the i
t
h parameter (R
s,R
M, and HR) on either the aortic valve opening duration or the cardiac output can be defined with a sensitivity index as follows:(28)
δ
w
p
i
=
w
g
i
+
-
w
g
i
0
w
g
i
0
+
w
g
i
-
-
w
g
i
0
w
g
i
0
,where δ
w
p
i is a sensitivity index that can be used to decide the significance of the i
t
h parameter.Based on the sensitivity index, a half-normal probability diagram is used [33] to identify the significant parameters; i.e., the factors induce significant changes in the aortic valve opening duration and the cardiac output. The half-normal probability is defined as(29)
Φ
-
1
0.5
+
0.5
i
-
0.5
k
,
δ
w
p
i
,wherei = 1, …,k is the i
t
h parameter and Φ
-
1 represents the cumulative distribution function (CDF) with respect to a standard normal distribution.For clarity, Figure5 shows the sensitivity analysis results for different electric pump powersP
E. Note that Figures 5(a)–5(c) are the results for the aortic valve opening duration, while Figures 5(d)–5(f) are the results for the cardiac output. The electric pump power in Figures 5(a) and 5(d) was set to 0.12W, and 0.36W was used in Figures 5(b) and 5(e). For Figures 5(c) and 5(f), the pump power was set to a relatively larger value of 0.6W. It is worth mentioning that the pump motor power was chosen to avoid the permanent closure of the aortic valve.Figure 5
Half-normal probability plots for sensitivity analysis: (a) results of aortic valve opening duration withP
E = 0.12W; (b) results of aortic valve opening duration with P
E = 0.36W; (c) results of aortic valve opening duration with P
E = 0.6W; (d) results of cardiac output with P
E = 0.12W; (e) results of cardiac output with P
E = 0.36W; (f) results of cardiac output with P
E = 0.6W.
(a)
(b)
(c)
(d)
(e)
(f)As seen in Figure5, it was found that systemic vascular resistance (SVR), i.e.,R
s, is the most significant uncertainty for all the case studies except in Figure 5(a). To validate the results, additional case studies were performed to evaluate the effect of variations inR
s and HR on the aortic valve opening duration with respect to a relatively larger range of pump powers, and the results were summarized in Figure 6.Figure 6
Sensitivity index ofR
s and HR corresponding to pump motor power.As seen in Figure6, it was found that the sensitivity index ofR
s increases as the pump power P
E increases, whereas the sensitivity index of HR decreases slowly when P
E is increased. This indicates that the hemodynamic such as the aortic valve opening is more sensitive to variations in the level of activity of a HF patient, i.e.,R
s, especially when the electric pump power is higher. Since the pump is often operated above the minimum value of the pump power, we will focus on the level of activity R
s in this work for the rest of the study. It is important to note that similar results were previously reported by clinical studies; i.e., a rotary LVAD pump is more sensitive to the afterload of the left ventricle [26, 27].Based on discussion,R
s is identified as the most sensitive factor in the cardiovascular-LVAD system that can affect the cardiac hemodynamic. Following the procedures as discussed in Section 2.2, a gPC model is developed to study the effect of the variations in R
s on the aortic valve opening duration and the cardiac output.


### 3.3. Formulation of the Stochastic gPC Model

It is assumed that systemic vascular resistance (SVR), i.e.,R
s in this work, follows a normal distribution. Three mean values ofR
s are used for algorithm verification, i.e., 0.5 mmHg/ml/s, 1 mmHg/ml/s, and 2 mmHg/ml/s, which represent different levels of physical activity, i.e., very active (0.5 mmHg/ml/s), moderately active (1 mmHg/ml/s), and inactive (2 mmHg/ml/s). For example, the very active state represents that the patient is climbing stairs, while the inactive state means the patient is resting or sleeping. To build a gPC model of SVR, the order of the polynomial chaos expansion is set to 1, i.e.,q= 1 in (21), since the uncertainty is normally distributed. Using (24), the expansion of each state variable in Table 2 would involve 2 terms, i.e.,Q = 1, since there is one uncertainty n
g
=
1 and the highest order of polynomial expansion of uncertainty is 1 (q=1).To introduce perturbations in SVR, a 10% variation around each mean value is used, which can be used to determine the gPC coefficients ofR
s in (21). The gPC coefficients of the hemodynamic variables (x
1-x
6) can be calculated by substituting the gPC models ofR
s and each state variable into (5) and by using a Galerkin projection, from which a stochastic model can be formulated. This model can describe dynamic behaviors of the cardiovascular-pump system in the presence of uncertainty in SVR (R
s). For brevity, the stochastic model is presented in Appendix B. The simulation results of the stochastic cardiovascular-pump model are shown in Figure 7.Figure 7
Hemodynamic waveforms in the presence of uncertainty in SVR. For the first and second column, the solid blue lines represent the mean value of each state variable in (5), while the red dotted line represents the 2
n
d-order gPC coefficients that can be used to calculate the standard deviation (σ in the last column) at each time interval. The last column shows the mean, maximum, and minimum value of each state variable, i.e., left ventricular pressure, aortic pressure, aortic flow, and pump flow.The first and second columns in Figure7 show the gPC coefficients of the state variables in (5), i.e., the left ventricle pressure (x
1), the aortic pressure (x
4), and the flows of the pump and aorta (x
5 andx
6), respectively. For the simulations of 20 cardiac cycles, the first 3 cardiac cycles are given in the first column, whereas the second column shows the last two cardiac cycles. The third column shows the variations in these states resulting from perturbations in SVR (R
s). For the results shown in Figure 7, the HR was set to 75 bpm, and the mean value ofR
s was 1 mmHg/ml/s. The maximum elastane (E
max) was set to 1 mmHg/ml, which represents a native heart with mild heart failure, and the pump motor power P
E used in this case study was 0.12W to avoid permanent closure of the aortic valve. It is important to note that the difference between the waveforms of the pump flow (x
6) and the aortic flow (x
5) in each cardiac cycle can be used to estimate the blood flow ejected by a native heart.As seen in the third column of Figures7(c) and 7(f), the left ventricular pressure (x
1) and the aortic pressure (x
4) can be affected by perturbations in SVR (R
s). In addition, it was found that the resulting variation in the aortic pressure is more significant than the left ventricular pressure. Note that the variations inx
1 increase as the electric power increases; however, this is not shown for brevity. Using the gPC coefficients, it is possible to estimate the upper and lower bounds of all state variables at each time interval. For example, Figure 7(f) shows the results of aortic pressure, where σ is the standard deviation calculated with the higher-order gPC coefficients in (27), i.e., 0
<
k
≤
Q. The range defined by the upper and lower bounds can quantify the dynamic values of aortic pressure within three standard deviations of the mean value of aortic pressure, which shows the 99% confidence interval of the aortic pressure at a particular time interval.As explained in Section3.1, the aortic valve opens when the left ventricular pressure is higher than the aortic pressure. Using the upper and lower bounds of the left ventricular pressure and aortic pressure, the effect of perturbations in R
s on the aortic valve opening duration can be quickly estimated. In addition, the mean and the variance in the pump flow can be rapidly calculated with (26) and (27) as shown in Figures 7(i) and 7(l), from which the variation in the cardiac output can be estimated. The calculation of the variance in the aortic valve opening duration and the cardiac output will be further discussed in the following sections.


### 3.4. Dynamics of Aortic Valve Flow and Pump Flow

The effect of perturbations in the SVR (R
s) on the state variables has been discussed in the previous section, where the electric pump power P
E was fixed. In this case study, the main objective is to investigate the effect of uncertainty on the aortic valve dynamics of the heart supported by a LVAD with a time-varying pump power. For this purpose, the electric pump power P
E was changed from 0.12 to 1.56W in order to better study the effect of the pump speed on the aortic valve over a wide range of operating conditions. The aortic valve opening duration can be determined by examining the difference between the left ventricular pressure and aortic pressure, as explained above [13]. When the aortic valve opens, a certain portion of the blood from the left ventricle flows through the aortic valve. As the pump motor power increases, the portion of blood flowing through the aortic valve will be decreased. It was found that when the pump power reaches a certain level, the LVAD takes over heart function and the aortic valve can be fully bypassed, i.e., permanently closed. Figures 8 and 9 show the changes in aortic valve flow and pump flow for different pump power values (P
E) in the presence of perturbations in SVR (R
s).Figure 8
Hemodynamic waveforms of aortic valve flow in the presence of uncertainty in SVR (R
s).
(a)
(b)Figure 9
Hemodynamic waveforms of pump flow in the presence of uncertainty in SVR (R
s).
(a)
(b)Due to the pulsatility in hemodynamic waveforms, Figures8(a) and 9(a) show the maximum, minimum, and mean values of aortic valve flow and the pump flow corresponding to different pump powers, respectively. It is worth mentioning that the maximum, minimum, and mean values of the state variables, i.e., aortic valve flow and pump flow in Figures 8(a) and 9(a), are calculated from the first gPC coefficient in (22) of each state variable, i.e., x
^
j
,
k
=
0. In contrast, the bar plot, representing the confidence interval, is calculated with other higher-order gPC coefficients, i.e., x
^
j
,
k (0
<
k
≤
Q), using (27).Figures8(b) and 9(b) show the simulation results of the hemodynamic waveforms of the aortic valve flow and the pump flow for a specific cardiac cycle, when different pump powers P
E were used. The variances around the maximum and minimum values in both flows through the aortic valve and pump are calculated by using the high order of gPC coefficients corresponding to the maximum values, while the variances of mean values for both flows are obtained based on the upper and lower bounds quantified with standard deviation calculated by gPC coefficients. It is important to note that the pump power P
E can be automatically selected according to the level of activity in order to meet different physiological demands, but this is not discussed in this work, as our objective is to evaluate the dynamic behaviors of the aortic valve over a wide range of pump power.As seen in Figure8(a), it was found that the maximum and mean values of the aortic valve flow decrease with the increase of the pump motor power P
E. Note that the aortic valve was completely closed when the pump power P
E was increased to 0.84W, which can be defined as abreakpoint. In addition, as seen in Figure 8(a), when the pump power P
E is below thebreakpoint, negative values of the minimum aortic valve flow were observed. This can possibly be attributed to regurgitation of flow through the aortic valve, because of aortic compliance. Physiologically, this phenomenon can be caused by the adverse pressure gradient, which can be developed as the aortic flow starts to decelerate quickly after reaching its maximum. This can eventually affect the low momentum fluid near the wall of the aorta, thus inducing reverse flow in the sinus region, which was previously reported [34]. It is important to note that the maximum and mean values of aortic valve flows remain positive. The reverse flow only exists in the minimum values of aortic valve flow, which is negligible and can be decreased as the pump power is increased. After thebreakpoint, i.e., P
E = 0.84W, the aortic flow is reduced to approximately 0 L/min, which means that the LVAD takes over the left ventricle function for the entire cardiac cycle and there is negligible blood that can flow through the aortic valve. Further, as can be seen in Figure 8(b), the pulsatility of the aortic flow decreases as the pump power is increased. When the pump power P
E reaches thebreakpoint, the aortic valve flow can completely lose the pulsatility. Such information can be useful for the controller design to adjust the pump speed to meet the different physiological demands.As shown in Figure9(a), the maximum, minimum, and mean values of pump flow rate increase when the pump power is increased. The maximum values of the pump flow increase until the pump power reaches thebreakpoint and then converges to a constant after thebreakpoint
P
E = 0.84W, while the minimum and mean values of the pump flow keep increasing as the pump power increases. In addition, it is important to note that the dashed lines (purple and red) in Figure 9(b) represent the variations in the pump flow, resulting from the perturbations in SVR (R
s). As seen, when the pump takes over the native heart, the variation in pump flows becomes larger, as compared to the cases when the aortic valve operates normally (i.e., before reaching thebreakpoint). Further, it is worth mentioning that the variations of the aortic valve flow and pump flow shown in Figures 8 and 9 are correlated to the pump motor power. It was found that the aortic valve flow can be significantly affected by perturbations inR
s before thebreakpoint. In contrast, the variation in pump flow is larger when the LVAD pump begins to take over the left ventricular function after thebreakpoint. This observation can be used as a tuning constraint of pump power in the controller design, since larger variations in the pump flow and the aortic flow may weaken the myocardium, which is detrimental to cardiac recovery and can be fatal to HF patients with LVADs. For example, constraints can be used to confine the allowable tuning range of the pump power to avoid inducing larger variations in pump flow, while taking perturbations in SVR (R
s) into account.


### 3.5. Aortic Valve Opening Duration

In Section3.4, it was found that the aortic valve flow rate varies with respect to the electric pump power. The aortic valve can be closed when the pump power reaches a certain level (e.g.,breakpoint), at which the LVAD pump takes over heart function. Based on this observation, the objective in this case study is to examine the aortic valve opening duration in a cardiac cycle, while considering the perturbations in SVR (R
s).As explained in Section3.1, the aortic valve opening duration is measured by calculating the ratio between the time that the aortic valve opens during a cardiac cycle and the duration of the cardiac cycle. For clarity, the mean values of the aortic valve opening duration for a mild HF patient (E
max = 1) and a severe HF patient (E
max = 0.5) with different electric pump powers are shown in Figures 10 and 11, respectively. In addition, the variations around the mean values of aortic valve opening duration are summarized in Tables 3 and 4 for brevity. In this study, HR is set to 75 bpm for both cases, and R
s follows the same probabilistic description as done in Section 3.3, i.e., three different levels of physical activity of a patient.Table 3
Aortic valve opening duration calculated from gPC coefficients for a mild heart failure patient.
Pump power,P
E (watts)
0.12
0.36
0.6
0.84
1.08
1.32
1.56
R
s = 0.5
Opening duration (%)
29.70±0.80
26.47±0.94
23.51±1.08
20.40±1.23
16.88±1.41
12.55±1.70
0±9.20
R
s = 1.0
Opening duration (%)
23.39±0.74
18.23±0.91
13.15±1.24
-
-
-
-
R
s = 2.0
Opening duration (%)
16.04±0.66
0±4.66
-
-
-
-
-Table 4
Aortic valve opening duration calculated from gPC coefficients for a severe heart failure patient.
Pump power,P
E (watts)
0.12
0.36
0.6
0.84
1.08
1.32
1.56
R
s = 0.5
Opening duration (%)
26.27±0.85
19.89±0.74
12.06±1.03
-
-
-
-
R
s = 1.0
Opening duration (%)
18.54±0.57
0±1.62
-
-
-
-
-
R
s = 2.0
Opening duration (%)
8.65±0.57
-
-
-
-
-Figure 10
Aortic valve opening duration with respect to the level of physical activity and pump power for a mild heart failure patient (E
max = 1).Figure 11
Aortic valve opening duration with respect to the level of physical activity and pump power for a severe heart failure patient (E
max = 0.5).As seen in Figures10 and 11, aortic valve opening duration highly depends on the level of activity for a patient and the severity level of heart failure. For example, as seen in Figure 10, when the patient is very active, the aortic valve can remain open in a cardiac cycle over a wide range of pump motor power. In contrast, the aortic valve can be easily taken over by the pump when a patient is in an inactive state and when the pump power is higher than 0.36W. In addition, it was found that the operating range of the pump is much smaller, when the level of heart failure is severe. For example, as seen in Figure 11, for a very active patient with severe HF, the aortic valve will be closed when the pump power is larger than about 0.7W. However, the aortic valve can remain open for the whole range of the pump power for a very active mild HF patient (see Figure 10). Thus, thebreakpoint, where the aortic valve closes permanently, can be determined with respect to the level of physical activity as well as the severity of HF. These observations can provide useful information when the pump motor power is used as a control variable. For example, it is difficult to control the pump power to maintain the aortic valve open for a severe heart failure patient, since the allowable operating range of the pump power to keep the aortic valve open is small. This may reduce the cardiac perfusion required to meet the physiological demands. For severe HF patients, an appropriate controller can be developed by considering the severity of HF and by taking into account possible levels of activity, which finds a trade-off between the desired cardiac output and the aortic valve opening duration.To validate the results, the aortic valve opening duration was also compared with previously reported work. As mentioned in [35], the average opening duration of the aortic valve is 30.5% for a native healthy heart without a LVAD, 27% for a mild HF patient, and 25% for a severe HF patient. This clearly shows that our results are in good agreement with these clinical observations. In addition, it was found that in this work the aortic valve opening duration is highly related to the level of activity and the severity of the native heart.In addition, Tables3 and 4 show the confidence interval of the aortic valve opening duration. It was found that the variance in the aortic valve opening duration increases, as the mean value ofR
s decreases and as the pump power increases. For example, as seen in Table 3, for an active and mild heart failure patient (i.e., R
s = 0.5 mmHg/ml/s), the variations in the aortic valve opening duration increase as the pump power increases. It was found that the variation in aortic valve opening duration is about 2 percent point of the mean values on average for these possible pump powers shown in Table 3. However, for an active and severe heart failure patient as seen Table 4 (R
s = 0.5 mmHg/ml/s), the average of the variation in the aortic valve opening duration is approximately 0.8%.


### 3.6. Cardiac Output

Cardiac output (CO) is conventionally calculated as the product of the stroke volume (V) and the HR asC
O
=
V
×
H
R. However, for HF patients with implanted LVADs, the calculation of cardiac output needs to consider the blood ejected by both the native heart and the LVAD pump. To accomplish this, the cardiac output can be calculated by integrating the aortic flow (x
5) [8]. For a cardiac cycle, (30)~(32) can be used to estimate the cardiac output as (30)
V
T
=
V
P
+
V
h
(31)
∫
t
t
+
t
c
x
5
ξ
=
∫
t
t
+
t
c
x
6
ξ
+
V
h
(32)
C
O
T
=
C
O
P
+
C
O
h
,where V
T in (30) is the total blood volume pumped into aorta and V
P and V
h are the blood volume ejected by the LVAD pump and the native heart, respectively. The blood volume can be calculated with (31) for a cardiac cyclet
c, from which the total cardiac output of the aorta can be approximated with the summation of two cardiac outputs generated by the pump C
O
P and the native heart C
O
h, by multiplying both sides of (31) with the HR. Note that since the perturbations in SVR (R
s) are considered in this work, both the pump flow and the aortic flow are functions of a random event ξ that can be used to approximate the variation in SVR. Thus, the integration in (31) is calculated over the domain defined by ξ. To solve (31), the trapezoidal rule is used in this work. Using the gPC model of state variable such asx
5 andx
6, it is possible to provide a measure of confidence interval in the cardiac output prediction.Note thatV
h in (31) is the blood volume of the native heart, which can be estimated from a PV-loop using the gPC models. Figure 12 shows the PV-loops of a mild HF patient with different levels of activity. In this case, the pump power is set to 0.12Wand the HR is 75 bpm. As seen in Figure 12, any perturbations inR
s can affect the stroke work, thus affecting the cardiac output.Figure 12
PV-loops generated by the developed stochastic model using gPC theory (E
max = 1.0, P
E = 0.12W, and HR = 75 bpm).In addition, the results of the total cardiac output with respect to different mean values ofR
s and pump motor powers are shown in Figures 13 and 14 for a mild HF patient and a severe HF patient, respectively. For simplicity, only three pump powers were investigated, i.e., 0.12, 0.84, and 1.56W, respectively. These values were chosen according to the dynamic behavior of the aortic valve. For example, the smallest value of the pump power, i.e., 0.12W, is used to ensure the aortic valve can operate normally. The medium value of the pump power, i.e., 0.84W, is used to study the cardiac output when the pump is operated at thebreakpoint as discussed before. In contrast, the largest value of pump power, i.e., 1.56W, is used to study the effect of perturbations on cardiac output, when the pump takes over heart function.Figure 13
Simulated cardiac output with respect to mean values ofR
s and pump motor power for mild heart failure (E
max = 1.0).Figure 14
Simulated cardiac output with respect to mean values ofR
s and pump motor power for severe heart failure (E
max = 0.5).In Figures13 and 14, the verticalbar represents the mean value of the cardiac output, while theerror-plot in each vertical bar represents the variation around the mean value. As seen, the cardiac output and the variation in the cardiac output decrease as the mean value ofR
s increases, since the HF patient is less active, and less blood is required. In addition, it was found that the cardiac output increases as the pump power increases. Further, it was found that the variation in the cardiac output decreases as heart function gets weaker. For example, as seen in Figure 13, the variation in the cardiac output for an active and mild heart failure patient is about 6.1 L/min, when the pump power is set to 0.12W. In contrast, the variation in cardiac output for an active and severe heart failure patient is approximately 4.2 L/min with the same pump power as seen in Figure 14.


## 3.1. Cardiac Hemodynamic

The cardiovascular-pump model as explained in Section2.1 involves a few model parameters, which can be used to represent the physiological dynamics of the heart supported with a LVAD. As previously reported, parametersR
s,R
M, and heart rate (HR) can govern the behavior of the cardiovascular-pump system [13]. Thus, we propose to investigate their effect on the dynamics of the aortic valve and the cardiac output.Since the gPC-based stochastic model in this work is developed based on a deterministic model in (5), the first step is to validate the ability of the deterministic model in order to mimic the hemodynamic of the cardiovascular-pump system. The linear relationship between the end-systolic pressure and the left ventricle volume in the presence of perturbations inR
s andR
M is used. Note that the perturbations inR
s andR
M represent the stochastic changes during the preload and afterload of the heart. A total of 4 preload and 4 afterload changes were simulated as shown in Figure 3, for which the left ventricle parameters such as E
m
a
x, E
m
i
n, andV
0 are set to constant and the electric pump power is set to P
E = 0.12W.Figure 3
PV-loops for different values of (a) the resistanceR
s (afterload) and (b)R
M (preload).
(a)
(b)The left ventricle pressure and the left ventricle volume graph (PV-loop) are shown in Figure3. The PV-loops in Figure 3(a) show the hemodynamic changes resulting from variations in the afterload, i.e., different values ofR
s (SVR) that represent the changes in the level of activity. In contrast, the PV-loop in Figure 3(b) shows the result by altering the preload, i.e.,R
M representing mitral valve resistance. For both graphs, a linear relationship between the left ventricle pressure and the left ventricle volume is observed, thus confirming the ability of the model to describe the hemodynamic of the cardiovascular-pump system. In addition, it was found that any variations in the cardiovascular-pump system can significantly affect the PV-loop such as the area that represents the stroke work.Note that the dynamic behaviors of the aortic valve will be used in this work as a criterion to evaluate the effect of uncertainty on the cardiovascular-pump system. The rationale is that any changes in physiological parameters such asR
s can affect left ventricular pressure as seen in Figure 3(a). Constant low ventricular pressure can result in the permanent closure of the aortic valve, which is detrimental to cardiac recovery. Specifically, the aortic valve opening duration is used to quantitatively assess the effect of uncertainty on the dynamic behaviors of the aortic valve.The calculation of the aortic valve opening duration proceeds as follows. In a cardiac cycle, the period of time,t
open, during which the left ventricular pressure (LVP) is larger than the aortic pressure (AoP) will be first determined. The aortic valve opening duration can then be obtained by calculating the ratio betweent
open and the time of a cardiac cyclet
c. A schematic of the calculation of the aortic valve opening duration is given in Figure 4. Let us suppose there are 9 data points of bothLVP andAoP in a cardiac cycle, and the time interval between any two data points ist
int. As seen, two data points ofLVP are larger than theAoP, indicating thatt
open is 2×t
int. Similarly, the cardiac cycle can be calculated ast
c = 9×t
int. Thus, the aortic valve opening duration can be determined with 2/9 ≈ 0.222, representing that the aortic valve remains open during approximately 22.2% of a cardiac cycle.Figure 4
Schematic for the calculation of the aortic valve opening duration.The cardiac output (CO) is another key physical property, which can be calculated as a product of stroke volume and HR. The cardiac output represents the total amount of blood pumped from the ventricle in a minute. The stroke volume is the blood volume pumped from the left ventricle in a beat and is determined by the contractility, preload, and afterload. The area of the PV-loop is often calculated to represent the stroke work. As shown in Figure3, any changes in R
s andR
M can affect the total amount of blood entering the ventricle, thus affecting the stroke volume and the cardiac output. Therefore, the cardiac output is used as a second criterion to evaluate the effect of uncertainty on the cardiovascular-pump system in this work.


## 3.2. Sensitivity Analysis

As mentioned above,R
s,R
M, and HR can affect the physiological dynamics of the heart such as the aortic valve opening and the cardiac output. However, each factor may have different effect on the hemodynamic of the failing heart. Thus, a sensitivity analysis is first performed to identify the most significant physiological factor. The effect of variations inR
s,R
M, and HR on the aortic valve and the cardiac output is first investigated, using the cardiovascular-LVAD model with a mild HF patient (i.e.,E
max = 1.0) and electric pump power P
E = 0.12 ~ 0.6W. Similar results were found for other heart conditions and pump powers, but they are not shown for brevity.It is assumed that each parameter can vary randomly between two levels, i.e., +1 and -1, which correspond to a +10% change and a -10% change with respect to its nominal values. Note that the nominal values of each parameter areR
s = 1.0 mmHg∙s/ml,R
M = 0.005 mmHg∙s/ml, and HR = 75 bpm. For each parameter, let us define the aortic valve opening duration (or the cardiac output) as w
g
i
+ and w
g
i
- for two different levels and the corresponding result with the nominal values of model parameters as w
g
i
0. Then, the effect of uncertainty in the i
t
h parameter (R
s,R
M, and HR) on either the aortic valve opening duration or the cardiac output can be defined with a sensitivity index as follows:(28)
δ
w
p
i
=
w
g
i
+
-
w
g
i
0
w
g
i
0
+
w
g
i
-
-
w
g
i
0
w
g
i
0
,where δ
w
p
i is a sensitivity index that can be used to decide the significance of the i
t
h parameter.Based on the sensitivity index, a half-normal probability diagram is used [33] to identify the significant parameters; i.e., the factors induce significant changes in the aortic valve opening duration and the cardiac output. The half-normal probability is defined as(29)
Φ
-
1
0.5
+
0.5
i
-
0.5
k
,
δ
w
p
i
,wherei = 1, …,k is the i
t
h parameter and Φ
-
1 represents the cumulative distribution function (CDF) with respect to a standard normal distribution.For clarity, Figure5 shows the sensitivity analysis results for different electric pump powersP
E. Note that Figures 5(a)–5(c) are the results for the aortic valve opening duration, while Figures 5(d)–5(f) are the results for the cardiac output. The electric pump power in Figures 5(a) and 5(d) was set to 0.12W, and 0.36W was used in Figures 5(b) and 5(e). For Figures 5(c) and 5(f), the pump power was set to a relatively larger value of 0.6W. It is worth mentioning that the pump motor power was chosen to avoid the permanent closure of the aortic valve.Figure 5
Half-normal probability plots for sensitivity analysis: (a) results of aortic valve opening duration withP
E = 0.12W; (b) results of aortic valve opening duration with P
E = 0.36W; (c) results of aortic valve opening duration with P
E = 0.6W; (d) results of cardiac output with P
E = 0.12W; (e) results of cardiac output with P
E = 0.36W; (f) results of cardiac output with P
E = 0.6W.
(a)
(b)
(c)
(d)
(e)
(f)As seen in Figure5, it was found that systemic vascular resistance (SVR), i.e.,R
s, is the most significant uncertainty for all the case studies except in Figure 5(a). To validate the results, additional case studies were performed to evaluate the effect of variations inR
s and HR on the aortic valve opening duration with respect to a relatively larger range of pump powers, and the results were summarized in Figure 6.Figure 6
Sensitivity index ofR
s and HR corresponding to pump motor power.As seen in Figure6, it was found that the sensitivity index ofR
s increases as the pump power P
E increases, whereas the sensitivity index of HR decreases slowly when P
E is increased. This indicates that the hemodynamic such as the aortic valve opening is more sensitive to variations in the level of activity of a HF patient, i.e.,R
s, especially when the electric pump power is higher. Since the pump is often operated above the minimum value of the pump power, we will focus on the level of activity R
s in this work for the rest of the study. It is important to note that similar results were previously reported by clinical studies; i.e., a rotary LVAD pump is more sensitive to the afterload of the left ventricle [26, 27].Based on discussion,R
s is identified as the most sensitive factor in the cardiovascular-LVAD system that can affect the cardiac hemodynamic. Following the procedures as discussed in Section 2.2, a gPC model is developed to study the effect of the variations in R
s on the aortic valve opening duration and the cardiac output.


## 3.3. Formulation of the Stochastic gPC Model

It is assumed that systemic vascular resistance (SVR), i.e.,R
s in this work, follows a normal distribution. Three mean values ofR
s are used for algorithm verification, i.e., 0.5 mmHg/ml/s, 1 mmHg/ml/s, and 2 mmHg/ml/s, which represent different levels of physical activity, i.e., very active (0.5 mmHg/ml/s), moderately active (1 mmHg/ml/s), and inactive (2 mmHg/ml/s). For example, the very active state represents that the patient is climbing stairs, while the inactive state means the patient is resting or sleeping. To build a gPC model of SVR, the order of the polynomial chaos expansion is set to 1, i.e.,q= 1 in (21), since the uncertainty is normally distributed. Using (24), the expansion of each state variable in Table 2 would involve 2 terms, i.e.,Q = 1, since there is one uncertainty n
g
=
1 and the highest order of polynomial expansion of uncertainty is 1 (q=1).To introduce perturbations in SVR, a 10% variation around each mean value is used, which can be used to determine the gPC coefficients ofR
s in (21). The gPC coefficients of the hemodynamic variables (x
1-x
6) can be calculated by substituting the gPC models ofR
s and each state variable into (5) and by using a Galerkin projection, from which a stochastic model can be formulated. This model can describe dynamic behaviors of the cardiovascular-pump system in the presence of uncertainty in SVR (R
s). For brevity, the stochastic model is presented in Appendix B. The simulation results of the stochastic cardiovascular-pump model are shown in Figure 7.Figure 7
Hemodynamic waveforms in the presence of uncertainty in SVR. For the first and second column, the solid blue lines represent the mean value of each state variable in (5), while the red dotted line represents the 2
n
d-order gPC coefficients that can be used to calculate the standard deviation (σ in the last column) at each time interval. The last column shows the mean, maximum, and minimum value of each state variable, i.e., left ventricular pressure, aortic pressure, aortic flow, and pump flow.The first and second columns in Figure7 show the gPC coefficients of the state variables in (5), i.e., the left ventricle pressure (x
1), the aortic pressure (x
4), and the flows of the pump and aorta (x
5 andx
6), respectively. For the simulations of 20 cardiac cycles, the first 3 cardiac cycles are given in the first column, whereas the second column shows the last two cardiac cycles. The third column shows the variations in these states resulting from perturbations in SVR (R
s). For the results shown in Figure 7, the HR was set to 75 bpm, and the mean value ofR
s was 1 mmHg/ml/s. The maximum elastane (E
max) was set to 1 mmHg/ml, which represents a native heart with mild heart failure, and the pump motor power P
E used in this case study was 0.12W to avoid permanent closure of the aortic valve. It is important to note that the difference between the waveforms of the pump flow (x
6) and the aortic flow (x
5) in each cardiac cycle can be used to estimate the blood flow ejected by a native heart.As seen in the third column of Figures7(c) and 7(f), the left ventricular pressure (x
1) and the aortic pressure (x
4) can be affected by perturbations in SVR (R
s). In addition, it was found that the resulting variation in the aortic pressure is more significant than the left ventricular pressure. Note that the variations inx
1 increase as the electric power increases; however, this is not shown for brevity. Using the gPC coefficients, it is possible to estimate the upper and lower bounds of all state variables at each time interval. For example, Figure 7(f) shows the results of aortic pressure, where σ is the standard deviation calculated with the higher-order gPC coefficients in (27), i.e., 0
<
k
≤
Q. The range defined by the upper and lower bounds can quantify the dynamic values of aortic pressure within three standard deviations of the mean value of aortic pressure, which shows the 99% confidence interval of the aortic pressure at a particular time interval.As explained in Section3.1, the aortic valve opens when the left ventricular pressure is higher than the aortic pressure. Using the upper and lower bounds of the left ventricular pressure and aortic pressure, the effect of perturbations in R
s on the aortic valve opening duration can be quickly estimated. In addition, the mean and the variance in the pump flow can be rapidly calculated with (26) and (27) as shown in Figures 7(i) and 7(l), from which the variation in the cardiac output can be estimated. The calculation of the variance in the aortic valve opening duration and the cardiac output will be further discussed in the following sections.


## 3.4. Dynamics of Aortic Valve Flow and Pump Flow

The effect of perturbations in the SVR (R
s) on the state variables has been discussed in the previous section, where the electric pump power P
E was fixed. In this case study, the main objective is to investigate the effect of uncertainty on the aortic valve dynamics of the heart supported by a LVAD with a time-varying pump power. For this purpose, the electric pump power P
E was changed from 0.12 to 1.56W in order to better study the effect of the pump speed on the aortic valve over a wide range of operating conditions. The aortic valve opening duration can be determined by examining the difference between the left ventricular pressure and aortic pressure, as explained above [13]. When the aortic valve opens, a certain portion of the blood from the left ventricle flows through the aortic valve. As the pump motor power increases, the portion of blood flowing through the aortic valve will be decreased. It was found that when the pump power reaches a certain level, the LVAD takes over heart function and the aortic valve can be fully bypassed, i.e., permanently closed. Figures 8 and 9 show the changes in aortic valve flow and pump flow for different pump power values (P
E) in the presence of perturbations in SVR (R
s).Figure 8
Hemodynamic waveforms of aortic valve flow in the presence of uncertainty in SVR (R
s).
(a)
(b)Figure 9
Hemodynamic waveforms of pump flow in the presence of uncertainty in SVR (R
s).
(a)
(b)Due to the pulsatility in hemodynamic waveforms, Figures8(a) and 9(a) show the maximum, minimum, and mean values of aortic valve flow and the pump flow corresponding to different pump powers, respectively. It is worth mentioning that the maximum, minimum, and mean values of the state variables, i.e., aortic valve flow and pump flow in Figures 8(a) and 9(a), are calculated from the first gPC coefficient in (22) of each state variable, i.e., x
^
j
,
k
=
0. In contrast, the bar plot, representing the confidence interval, is calculated with other higher-order gPC coefficients, i.e., x
^
j
,
k (0
<
k
≤
Q), using (27).Figures8(b) and 9(b) show the simulation results of the hemodynamic waveforms of the aortic valve flow and the pump flow for a specific cardiac cycle, when different pump powers P
E were used. The variances around the maximum and minimum values in both flows through the aortic valve and pump are calculated by using the high order of gPC coefficients corresponding to the maximum values, while the variances of mean values for both flows are obtained based on the upper and lower bounds quantified with standard deviation calculated by gPC coefficients. It is important to note that the pump power P
E can be automatically selected according to the level of activity in order to meet different physiological demands, but this is not discussed in this work, as our objective is to evaluate the dynamic behaviors of the aortic valve over a wide range of pump power.As seen in Figure8(a), it was found that the maximum and mean values of the aortic valve flow decrease with the increase of the pump motor power P
E. Note that the aortic valve was completely closed when the pump power P
E was increased to 0.84W, which can be defined as abreakpoint. In addition, as seen in Figure 8(a), when the pump power P
E is below thebreakpoint, negative values of the minimum aortic valve flow were observed. This can possibly be attributed to regurgitation of flow through the aortic valve, because of aortic compliance. Physiologically, this phenomenon can be caused by the adverse pressure gradient, which can be developed as the aortic flow starts to decelerate quickly after reaching its maximum. This can eventually affect the low momentum fluid near the wall of the aorta, thus inducing reverse flow in the sinus region, which was previously reported [34]. It is important to note that the maximum and mean values of aortic valve flows remain positive. The reverse flow only exists in the minimum values of aortic valve flow, which is negligible and can be decreased as the pump power is increased. After thebreakpoint, i.e., P
E = 0.84W, the aortic flow is reduced to approximately 0 L/min, which means that the LVAD takes over the left ventricle function for the entire cardiac cycle and there is negligible blood that can flow through the aortic valve. Further, as can be seen in Figure 8(b), the pulsatility of the aortic flow decreases as the pump power is increased. When the pump power P
E reaches thebreakpoint, the aortic valve flow can completely lose the pulsatility. Such information can be useful for the controller design to adjust the pump speed to meet the different physiological demands.As shown in Figure9(a), the maximum, minimum, and mean values of pump flow rate increase when the pump power is increased. The maximum values of the pump flow increase until the pump power reaches thebreakpoint and then converges to a constant after thebreakpoint
P
E = 0.84W, while the minimum and mean values of the pump flow keep increasing as the pump power increases. In addition, it is important to note that the dashed lines (purple and red) in Figure 9(b) represent the variations in the pump flow, resulting from the perturbations in SVR (R
s). As seen, when the pump takes over the native heart, the variation in pump flows becomes larger, as compared to the cases when the aortic valve operates normally (i.e., before reaching thebreakpoint). Further, it is worth mentioning that the variations of the aortic valve flow and pump flow shown in Figures 8 and 9 are correlated to the pump motor power. It was found that the aortic valve flow can be significantly affected by perturbations inR
s before thebreakpoint. In contrast, the variation in pump flow is larger when the LVAD pump begins to take over the left ventricular function after thebreakpoint. This observation can be used as a tuning constraint of pump power in the controller design, since larger variations in the pump flow and the aortic flow may weaken the myocardium, which is detrimental to cardiac recovery and can be fatal to HF patients with LVADs. For example, constraints can be used to confine the allowable tuning range of the pump power to avoid inducing larger variations in pump flow, while taking perturbations in SVR (R
s) into account.


## 3.5. Aortic Valve Opening Duration

In Section3.4, it was found that the aortic valve flow rate varies with respect to the electric pump power. The aortic valve can be closed when the pump power reaches a certain level (e.g.,breakpoint), at which the LVAD pump takes over heart function. Based on this observation, the objective in this case study is to examine the aortic valve opening duration in a cardiac cycle, while considering the perturbations in SVR (R
s).As explained in Section3.1, the aortic valve opening duration is measured by calculating the ratio between the time that the aortic valve opens during a cardiac cycle and the duration of the cardiac cycle. For clarity, the mean values of the aortic valve opening duration for a mild HF patient (E
max = 1) and a severe HF patient (E
max = 0.5) with different electric pump powers are shown in Figures 10 and 11, respectively. In addition, the variations around the mean values of aortic valve opening duration are summarized in Tables 3 and 4 for brevity. In this study, HR is set to 75 bpm for both cases, and R
s follows the same probabilistic description as done in Section 3.3, i.e., three different levels of physical activity of a patient.Table 3
Aortic valve opening duration calculated from gPC coefficients for a mild heart failure patient.
Pump power,P
E (watts)
0.12
0.36
0.6
0.84
1.08
1.32
1.56
R
s = 0.5
Opening duration (%)
29.70±0.80
26.47±0.94
23.51±1.08
20.40±1.23
16.88±1.41
12.55±1.70
0±9.20
R
s = 1.0
Opening duration (%)
23.39±0.74
18.23±0.91
13.15±1.24
-
-
-
-
R
s = 2.0
Opening duration (%)
16.04±0.66
0±4.66
-
-
-
-
-Table 4
Aortic valve opening duration calculated from gPC coefficients for a severe heart failure patient.
Pump power,P
E (watts)
0.12
0.36
0.6
0.84
1.08
1.32
1.56
R
s = 0.5
Opening duration (%)
26.27±0.85
19.89±0.74
12.06±1.03
-
-
-
-
R
s = 1.0
Opening duration (%)
18.54±0.57
0±1.62
-
-
-
-
-
R
s = 2.0
Opening duration (%)
8.65±0.57
-
-
-
-
-Figure 10
Aortic valve opening duration with respect to the level of physical activity and pump power for a mild heart failure patient (E
max = 1).Figure 11
Aortic valve opening duration with respect to the level of physical activity and pump power for a severe heart failure patient (E
max = 0.5).As seen in Figures10 and 11, aortic valve opening duration highly depends on the level of activity for a patient and the severity level of heart failure. For example, as seen in Figure 10, when the patient is very active, the aortic valve can remain open in a cardiac cycle over a wide range of pump motor power. In contrast, the aortic valve can be easily taken over by the pump when a patient is in an inactive state and when the pump power is higher than 0.36W. In addition, it was found that the operating range of the pump is much smaller, when the level of heart failure is severe. For example, as seen in Figure 11, for a very active patient with severe HF, the aortic valve will be closed when the pump power is larger than about 0.7W. However, the aortic valve can remain open for the whole range of the pump power for a very active mild HF patient (see Figure 10). Thus, thebreakpoint, where the aortic valve closes permanently, can be determined with respect to the level of physical activity as well as the severity of HF. These observations can provide useful information when the pump motor power is used as a control variable. For example, it is difficult to control the pump power to maintain the aortic valve open for a severe heart failure patient, since the allowable operating range of the pump power to keep the aortic valve open is small. This may reduce the cardiac perfusion required to meet the physiological demands. For severe HF patients, an appropriate controller can be developed by considering the severity of HF and by taking into account possible levels of activity, which finds a trade-off between the desired cardiac output and the aortic valve opening duration.To validate the results, the aortic valve opening duration was also compared with previously reported work. As mentioned in [35], the average opening duration of the aortic valve is 30.5% for a native healthy heart without a LVAD, 27% for a mild HF patient, and 25% for a severe HF patient. This clearly shows that our results are in good agreement with these clinical observations. In addition, it was found that in this work the aortic valve opening duration is highly related to the level of activity and the severity of the native heart.In addition, Tables3 and 4 show the confidence interval of the aortic valve opening duration. It was found that the variance in the aortic valve opening duration increases, as the mean value ofR
s decreases and as the pump power increases. For example, as seen in Table 3, for an active and mild heart failure patient (i.e., R
s = 0.5 mmHg/ml/s), the variations in the aortic valve opening duration increase as the pump power increases. It was found that the variation in aortic valve opening duration is about 2 percent point of the mean values on average for these possible pump powers shown in Table 3. However, for an active and severe heart failure patient as seen Table 4 (R
s = 0.5 mmHg/ml/s), the average of the variation in the aortic valve opening duration is approximately 0.8%.


## 3.6. Cardiac Output

Cardiac output (CO) is conventionally calculated as the product of the stroke volume (V) and the HR asC
O
=
V
×
H
R. However, for HF patients with implanted LVADs, the calculation of cardiac output needs to consider the blood ejected by both the native heart and the LVAD pump. To accomplish this, the cardiac output can be calculated by integrating the aortic flow (x
5) [8]. For a cardiac cycle, (30)~(32) can be used to estimate the cardiac output as (30)
V
T
=
V
P
+
V
h
(31)
∫
t
t
+
t
c
x
5
ξ
=
∫
t
t
+
t
c
x
6
ξ
+
V
h
(32)
C
O
T
=
C
O
P
+
C
O
h
,where V
T in (30) is the total blood volume pumped into aorta and V
P and V
h are the blood volume ejected by the LVAD pump and the native heart, respectively. The blood volume can be calculated with (31) for a cardiac cyclet
c, from which the total cardiac output of the aorta can be approximated with the summation of two cardiac outputs generated by the pump C
O
P and the native heart C
O
h, by multiplying both sides of (31) with the HR. Note that since the perturbations in SVR (R
s) are considered in this work, both the pump flow and the aortic flow are functions of a random event ξ that can be used to approximate the variation in SVR. Thus, the integration in (31) is calculated over the domain defined by ξ. To solve (31), the trapezoidal rule is used in this work. Using the gPC model of state variable such asx
5 andx
6, it is possible to provide a measure of confidence interval in the cardiac output prediction.Note thatV
h in (31) is the blood volume of the native heart, which can be estimated from a PV-loop using the gPC models. Figure 12 shows the PV-loops of a mild HF patient with different levels of activity. In this case, the pump power is set to 0.12Wand the HR is 75 bpm. As seen in Figure 12, any perturbations inR
s can affect the stroke work, thus affecting the cardiac output.Figure 12
PV-loops generated by the developed stochastic model using gPC theory (E
max = 1.0, P
E = 0.12W, and HR = 75 bpm).In addition, the results of the total cardiac output with respect to different mean values ofR
s and pump motor powers are shown in Figures 13 and 14 for a mild HF patient and a severe HF patient, respectively. For simplicity, only three pump powers were investigated, i.e., 0.12, 0.84, and 1.56W, respectively. These values were chosen according to the dynamic behavior of the aortic valve. For example, the smallest value of the pump power, i.e., 0.12W, is used to ensure the aortic valve can operate normally. The medium value of the pump power, i.e., 0.84W, is used to study the cardiac output when the pump is operated at thebreakpoint as discussed before. In contrast, the largest value of pump power, i.e., 1.56W, is used to study the effect of perturbations on cardiac output, when the pump takes over heart function.Figure 13
Simulated cardiac output with respect to mean values ofR
s and pump motor power for mild heart failure (E
max = 1.0).Figure 14
Simulated cardiac output with respect to mean values ofR
s and pump motor power for severe heart failure (E
max = 0.5).In Figures13 and 14, the verticalbar represents the mean value of the cardiac output, while theerror-plot in each vertical bar represents the variation around the mean value. As seen, the cardiac output and the variation in the cardiac output decrease as the mean value ofR
s increases, since the HF patient is less active, and less blood is required. In addition, it was found that the cardiac output increases as the pump power increases. Further, it was found that the variation in the cardiac output decreases as heart function gets weaker. For example, as seen in Figure 13, the variation in the cardiac output for an active and mild heart failure patient is about 6.1 L/min, when the pump power is set to 0.12W. In contrast, the variation in cardiac output for an active and severe heart failure patient is approximately 4.2 L/min with the same pump power as seen in Figure 14.


## 4. Conclusions

In this paper, a stochastic cardiovascular-LVAD model is developed to predict the dynamic behaviors of the aortic valve and cardiac output in the presence of uncertainty in systemic ventricular resistance (SVR). First, the effect of uncertainty in different model parameters on the aortic valve opening duration and the cardiac outputs was evaluated through a half-normal based sensitivity analysis. For the most sensitive model parameter, i.e., SVR representing the level of activity of a heart failure patient, a generalized polynomial chaos (gPC) expansion was used to approximate the perturbations around a set of mean values of SVR. Further, the effect of uncertainty on the cardiac hemodynamic such as the left ventricular pressure was evaluated using a Galerkin projection. It was found that perturbations in the level of activity (SVR) can significantly affect the aortic valve flow and the pump flow, which can consequently affect the aortic valve opening duration and the total cardiac output in each cardiac cycle. To ensure proper operation and avoid permeant closure of the aortic valve, an upper limit of the pump power (i.e.,breakpoint) is defined. As discussed in Results, thebreakpoint of the pump power should be specified with respect to the different levels of activity and different severities of heart failure patients. The more severe the heart failure, the lower the value of the pump power. In addition, the variation in the aortic valve opening duration and the cardiac output can be quickly calculated with the gPC model. This is useful for the control design to automatically adjust the pump power to meet the different physiological demands of the human body, since larger variation in the aortic valve opening and cardiac output may weaken the myocardium, which is detrimental to cardiac recovery and can be fatal to HF patients with LVADs. The understanding of the contribution of the pump in the overall task of the ejected blood can be very useful for developing a reliable and adaptive controller for the LVAD pump.



---

*Source: 7179317-2019-01-10.xml*