Modeling the potential impact of rectal microbicides to reduce HIV transmission in bathhouses

1.   Introduction

The cardiovascular system is a complex dynamic system, which contains rich dynamics information to predict the physiological or pathological state. The traditional method for studying cardiovascular hemodynamics is to perform surgery on animals or humans, which is usually limited by ethics and experimental conditions ^[1,2]. As the method and tool in vitro, the numerical model and mock circulatory model can provide technical support for clinical noninvasive diagnosis ^[3]. The numerical modeling is the most widespread research technique due to its proven capability, efficiency and minimized cost in comparison to in-vitro studies in mock flow loop, or in-vivo studies in animals and humans ^[4].

The cardiovascular lumped parameter (CLP) model is a simplified numerical method to study hemodynamics. Researchers have done a lot of work in the CLP modeling in the past ten years. Some lumped parameter models for a single circulation loop were modeled ^[5,6,7]. For example, Figliola et al. ^[5] modeled a lumped parameter model including only right heart and pulmonary circulation. Simaan et al. ^[7] modeled a fifth-order lumped parameter circuit model for left heart and systemic circulation. These simplified Windkessel models can be used to simulate the hemodynamics of a single circulation loop. Some complete CLP models have been developed ^{[8,9,10,11,12,13]}. For example, a CLP model including systemic circulation, pulmonary circulation, and the heart with four chambers was proposed ^[8]. To simulate the hemodynamics of patients with heart failure, Gu et al. ^[9] established the CLP for patients with heart failure, and verified the effectiveness of physiological mechanism by the model. Abdi ^[10] used a lumped method to institute a cardiovascular model with 42 chambers, including arteries, veins, capillaries, and heart chambers. This model was used to evaluate the effect of the abnormal heartbeat on cardiovascular system performance. Although these models are capable of simulating the hemodynamics of the intact cardiovascular system, they are not combined with ventricular assist devices.

Since their first surgical implantation in the 1960s ^[14], ventricular assist devices play an active role in the treatment of heart failure. The new generation of LVADs use axial or centrifugal rotary pumps ^[15]. Compared with the first generation of pulsating pumps, the rotary pumps have the advantages of small size, high reliability, and durability. It is helpful for clinical management to comprehend the coupling relationship between the native heart and the rotary pump. Coupling modeling of the cardiovascular system and ventricular assist device (VAD) is a feasible way to explain the hemodynamic relationship between a native heart and a rotary pump. Shi et al. ^[4,14] proposed a coupled model for the cardiovascular system and VAD to simulate the hemodynamics under left ventricular (LV) failure with VAD support. CARDIOSIM^© is a developed simulator based on the CLP model by Lazzari ^{[15,16,17,18]}. It not only can be used to simulate the physiological and pathological hemodynamics, but also can be combined with various types of ventricular assist devices (pulsating pumps, rotary pumps, etc.). CARDIOSIM^© has also been applied in medical education ^[19]. Lim et al. ^[1] developed a coupled model between the cardiovascular system and implantable rotary blood pump, and it includes the left heart, right heart, systemic circulation, pulmonary circulation, and implantable rotary blood pump. Di et al. ^[20,21,22] established a complete CLP model, which could be connected with pulsating pumps or rotary pumps. In the case of cardiac dysplasia or heart failure, it is helpful to understand the effects of single or two VADs on hemodynamics. However, the autonomic regulation of the cardiovascular system is not considered in the above models.

Baroreflex is a kind of short-term autonomic nerve regulation, which plays a great role in short-term pressure control ^[23]. The work of Ursino ^[24] and Ottesen ^[25] indicates that the baroreflex system can maintain the stability of the hemodynamic environment by regulating neural effectors. Consequently, it is more autonomous and accurate for simulating hemodynamics that a combination of the CLP model and the baroreflex model to form a closed-loop control system. Certainly, baroreflex regulation was also considered in some coupled models for the cardiovascular system and VADs ^{[26,27,28,29,30,31,32]}. In these models, ventricular characteristics or systemic vascular resistance are usually altered to simulate the hemodynamics of the cardiovascular system with impaired heart function. Bozkurt et al. ^[27] reduced left ventricular systolic elastance to 0.5 mmHg/mL from 2.5 mmHg/mL to simulate a dilated cardiomyopathy (DCM) condition in the cardiovascular system model. Cox et al. ^[28] changed left ventricular contractility, initial left ventricular volume and systemic vascular resistance to simulate end-stage DCM. Lim et al. ^[32] decreased the end-systolic elastance value to 40% of normal to implement a heart-failure model. Heart failure is a complex syndrome, it is not only related to the cardiovascular system but also the baroreflex system. These models only consider ventricular characteristics or systemic vascular resistance to simulate heart failure, which is different from clinical observations. However, other hemodynamic parameters in the cardiovascular system, such as systemic vascular compliance and pulmonary vascular resistance, as well as some parameters in the baroreflex system, should also be considered. Otherwise, some changes in hemodynamic parameters during LVAD support and the prediction results of the LVAD operating condition will be affected. The application scope of the above model is limited. Therefore, it is necessary to develop a comprehensive numerical model and consider the mechanism of changes in the cardiovascular system and the baroreflex system to improve the simulation and prediction effect.

In this study, we aim to establish a CLP model coupled with the baroreflex model and the LVAD model to form a comprehensive numerical model, which provides a platform for the hemodynamics simulation and the noninvasive evaluation of LVAD operating condition.

This paper is organized in the following manner: The modeling of the coupling cardiovascular system is introduced in Section 2. Section 3 presents some results including the hemodynamics under different conditions and the hemodynamic variables at different pump speeds. Section 4 discusses the operating condition of HeartMate Ⅲ and its effect on hemodynamics. In addition, the limitations of this study are also described in Section 4. Lastly, conclusions are drawn on the outcomes of this study in Section 5.

2.   Materials and method

The comprehensive numerical model is illustrated in Figure 1. The numerical model was made up of three parts: a developed CLP model for describing the human cardiovascular circulation system, a baroreflex model for short-term neuroregulation, and an LVAD model for assisting diseased heart. The rotary pump was connected with the diseased heart in parallel. Thus, the pump inlet cannulae was connected to the left ventricular apex and the outlet cannulae was connected to the aorta.

2.1.   CLP Model

The lumped parameter method is to concentrate on the human cardiovascular circulation system into several chambers with organ and tissue functions. Each chamber is composed of Windkessel elements representing hemodynamic parameters. This method can not only reduce the number of parameters to improve computational efficiency, but also obtain global hemodynamics. In this paper, the cardiovascular circulation system was divided into four parts: left heart, right heart, systemic circulation, and pulmonary circulation. The heart is a four-chamber pump with time-varying elasticity. The direction of blood flow is controlled by the four valves of the heart. Each chamber of heart was modeled as a Windkessel model.

2.1.1.   Ventricles

Suga et al. ^[33] show that the pressure-volume relationship of the ventricle can be expressed by a time-varying ventricle elastance. The ventricle elastance characterizes the elastic of ventricular myocardium and describes the dynamics of ventricular systole and diastole. Numerically, the ventricle elastance is generally defined as the ratio of pressure and volume in the ventricle. In the research field of the cardiovascular model, the ventricle elastance has been recognized and used ^{[1,4,7,11,13,18,21]}.

The time-varying elastances of the left and right ventricles:

where ${E_{lv}}(t)$ and ${E_{rv}}(t)$ are the time-varying elastances of left and right ventricles, ${C_{lv}}(t)$ and ${C_{rv}}(t)$ are the time-varying compliances of left and right ventricles, which are the reciprocals of ${E_{lv}}(t)$ and ${E_{rv}}(t)$, ${P_{lv}}(t)$ and ${P_{rv}}(t)$ are pressures of left and right ventricles, ${V_{lv}}(t)$ and ${V_{rv}}(t)$ are volumes of left and right ventricles, ${V_{lv, 0}}$ and ${V_{rv, 0}}$ are theoretical volumes in the left and right ventricles at zero pressure.

However, when the pressure and volume of the ventricle are not be provided, the time-varying elastance cannot be determined. Accordingly, the calculation formulas commonly used in engineering are as follows ^[7]:

where ${E_{n, lv}}({t_n})$ and ${E_{n, rv}}({t_n})$ are the normalized elastic functions. The actions of the left and right ventricles are synchronous, therefore, the same normalized elastic function is used for the left ventricle and right ventricle.

where $t$ and $T$ are time variable and heart period, ${\rm{HR}}$ is heart rate, ${E_{\max, lv}}$ and ${E_{\min, lv}}$ are left ventricular end-systolic elastance and left ventricular end-diastolic elastance, ${E_{\max, rv}}$ and ${E_{\min, rv}}$ are right ventricular end-systolic elastance and right ventricular end-diastolic elastance.

2.1.2.   Atria

The atrium also contracts and relaxes periodically similar to the ventricle, but its contractility is much smaller than that of the ventricle. The main function of the atrium is to store venous blood. Similarly, the time-varying elastance was adopted to describe the pressure-volume relationship of the atrium. The states of the left and right atria are modeled by the time-varying compliances ${C_{la}}(t)$ and ${C_{ra}}(t)$. ${E_{la}}(t)$ and ${E_{ra}}(t)$ are the time-varying elastances of left and right atria. ${E_{la}}(t)$ and ${E_{ra}}(t)$ can be expressed by Eqs (8) and (9). The meanings of other parameters have resemblance with those for ventricles.

Although the activities of the atrium and ventricle are asynchronous, the formulas for calculating atrial elastances are similar to those for the ventricles. The expressions of ${E_{la}}(t)$ and ${E_{ra}}(t)$ are as follows ^[14]:

where ${E_{\max, la}}$ and ${E_{\min, la}}$ are left atrial end-systolic elastance and left atrial end-diastolic elastance, ${E_{\max, ra}}$ and ${E_{\min, ra}}$ are right atrial end-systolic elastance and right atrial end-diastolic elastance, ${e_{la}}(t)$ and ${e_{ra}}(t)$ are the activation functions of the left and right atria.

The synchronous property of left and right atria gives rise to the same activation functions for left and right atria.

where ${T_{asb}}$ is the beginning time of atrium contraction and ${T_{asb}} = 0.9T$, ${T_{ase}}$ is the end time of atrium contraction and ${T_{ase}} = 0.09T + {T_{asb}}$.

2.1.3.   Valves

The heart valves prevent blood from flowing back, it comprises the mitral and tricuspid valves between the atria and ventricles, as well as the aortic and pulmonary valves between the ventricles and arteries. The opening and closing of the valve are similar to the short circuit and open circuit of an ideal diode. Consequently, the circuit model for the valve can be represented by resistor in series ideal diode. Take the aortic valve in Figure 1 for example. The circuit model of the aortic valve is composed of aortic valve resistance ${R_{{\rm{AV}}}}$ and ideal diode ${D_{{\rm{AV}}}}$, and the opening or closing of the valve is controlled by the difference between the left ventricular pressure ${P_{lv}}$ and the aortic pressure ${P_{ao}}$. If ${P_{lv}}{\rm{ - }}{P_{ao}} > 0$, the aortic valve is opened and the flow rate through the valve is ${{({P_{lv}} - {P_{ao}})} / {{R_{{\rm{AV}}}}}}$, otherwise the aortic valve is closed and the flow rate is 0.

The states of the valve are determined by the transvalvular pressure difference. The function of transvalvular pressure difference $\gamma (\xi)$ is introduced:

The valve flow rate ${Q_V}$ is as follow:

where $\xi$ is the transvalvular pressure difference, ${Q_V}$ is the flow rate through the valve and the $V$ is the abbreviation of the valve.

2.1.4.   Systemic circulation and pulmonary circulation

The systemic circulation is separated into aorta, artery, arteriole, capillary, and vein segments. The aorta was regarded as a vascular chamber with resistance, inertance, and compliance. The artery segment was considered as a vascular chamber with resistance and compliance. The primary role of the arteriole and capillary is to provide resistance to blood flow. Arteriole and capillary were modeled as a vascular chamber with resistance. The vascular resistance is provided by the artery, arteriole, and capillary. All levels of veins were modeled as a vascular chamber. Collecting and storing blood are primary roles played by veins that are possessed of greater compliance. The model for pulmonary circulation is similar to that for systemic circulation.

There is an analogy between fluid network and electric network. An electrical analogy is typically used where pressure is analogous to voltage, flow rate is analogous to electric current, resistance is analogous to resistor, inertance is analogous to inductor, and compliance is analogous to capacitor ^[12,34]. On account of an electrical network is applied to the blood circulation system, this paper used the Kirchhoff's law to expound the coupling relationship between the branches of the blood circulation system. Using the notation given in Figure 1, the equations of the instantaneous pressure and flow rate for systemic circulation are given as follows.

The flow rate balance equations are:

The pressure balance equation is:

2.2.   Baroreflex model

According to Ursino et al. ^[24], the baroreflex model was conducted in the order of afferent neural pathway, efferent neural pathways and effectors. In this study, the baroreflex model was modified to control the corresponding effectors in the cardiovascular system: systemic vascular resistance, left ventricular end-systolic elastance, right ventricular end-systolic elastance, and heart period. The open-loop control system for baroreflex regulation is shown in Figure 2.

2.2.1.   Afferent neural pathway

Ursino et al. ^[24] described the relationship between arterial pressure and the activity of the afferent nerves by considering two main aspects arising from experiments. First, the relationship between arterial pressure and afferent nerve activity in static conditions exhibits a sigmoidal shape with upper saturation and lower threshold. Second, the arterial baroreceptors are sensitive not only to the arterial pressure but also to its rate of change. These two properties have been reproduced by representing the arterial baroreceptors plus the afferent neural pathway as a first-order linear differential equation, with a static gain and a rate-dependent gain, connected in series with a sigmoidal static function.

where ${\tau _p}$ and ${\tau _z}$ are the time constants, aortic pressure ${p_{ao}}$ is used as the input pressure to arterial baroreceptors, $\tilde p$ is the output variable of a first-order linear differentiator, ${f_{as}}$ is the frequency of spikes in the afferent fibers, ${f_{\min }}$ and ${f_{\max }}$ are the lower and upper saturation of the frequency discharge, ${p_n}$ is the pressure value at the central point of the sigmoidal functional, ${k_a}$ is a parameter with the dimension of pressure, related to the slope of the static function at the central point, and ${1 /{{k_a}}}$ denotes the baroreceptor sensitivity.

2.2.2.   Efferent neural pathways

The efferent pathway includes the efferent sympathetic pathway and efferent vagal pathway. Increasing the frequency of spikes in the afferent fibers results in a decrease in the frequency of the sympathetic fibers and in an increase in vagal activity. These relationships have been reproduced by using a monotonically decreasing monoexponential static curve for the sympathetic efferent activity and a monotonically increasing exponential curve, with upper saturation, for the vagal activity.

where ${f_{es}}$ is the frequency of spikes in the efferent sympathetic nerves, ${f_{es, \infty }}$, ${f_{es, 0}}$ and ${k_{es}}$ are constants with ${f_{es, 0}} > {f_{es, \infty }}$, ${f_{ev}}$ is the frequency of spikes in the efferent vagal fibers, ${f_{ev, \infty }}$, ${f_{ev, 0}}$ and ${k_{ev}}$ are constant parameters with ${f_{ev, 0}} > {f_{ev, \infty }}$, and ${f_{as, 0}}$ is the central value in Eq 20.

2.2.3.   Effectors

All the effectors in the cardiovascular model change in response to stimulation of sympathetic nerves. The response of a generic effector to sympathetic stimulation includes a pure delay, a monotonic logarithmic static function, and a linear first-order dynamics with a real time constant. The use of logarithmic relationships is justified because the effector static response exhibits a rapid slope at low sympathetic frequencies and then progressively flattens. Of course, the static function is monotonically increasing when the systemic vascular resistance and the end-systolic elastances (left and right ventricles) are considered but monotonically decreasing when heart period is considered. Finally, the first-order dynamics simulates the time required for the effectors to progressively complete its action.

where $\theta$ denotes the generic controlled parameter (${R_{svr}}$, ${E_{\max, lv}}$, ${E_{\max, rv}}$, or $T$), $\Delta \theta$ is the parameter change caused by sympathetic stimulation, ${\sigma _\theta }$ is the output of the static characteristic, ${\tau _\theta }$ and ${D_\theta }$ are the time constant and pure latency, ${f_{es, \min }}$ is a threshold for sympathetic stimulation, ${G_\theta }$ is a constant gain factor, and ${\theta _0}$ denotes generic controlled parameter in the absence of cardiac innervation.

Heart period is affected not only by the sympathetic but also by the vagal activity. However, the static relationship is assumed to be linear, i.e., heart period increases proportionally to vagal stimulation. The final heart period level is then computed by summing the positive changes induced by the vagal stimulation, the negative changes induced by the sympathetic stimulation, and a constant level representing heart period in the absence of cardiac innervation. Hence, the two mechanisms interact linearly in the heart period control. Of course, this is a crude simplification of the real sympathetic-parasympathetic interaction. However, it provides an acceptable approximation of heart rate changes under a variety of physiological conditions.

where $\Delta {T_s}$ and $\Delta {T_v}$ denote sympathetic and vagal stimulation, ${\sigma _{{T_s}}}$ and ${\sigma _{{T_v}}}$ are the sympathetic and vagal activities, ${\tau _{{T_s}}}$ and ${D_{{T_s}}}$ are the time constant and pure latency for sympathetic, ${\tau _{{T_v}}}$ and ${D_{{T_v}}}$ are the time constant and pure latency for vagus, ${G_{{T_s}}}$ and ${G_{{T_v}}}$ are constant gain factors, and ${T_0}$ denotes heart period in the absence of cardiac innervation.

2.3.   LVAD model

A continuous flow LVAD model was connected to the left ventricular and aortic of the CLP model in Figure 1. It added one state variable to the original model, namely the pressure difference between the pump outlet and pump inlet $H$ as a function of pump flow rate and speed ^[4]:

where ${Q_{irp}}$ is the pump flow rate (mL/s), $\omega$ is pump speed (revolutions per minute, rpm), and ${\beta _a}$, ${\beta _b}$, ${\beta _c}$ are coefficients associated with a rotary pump. In this study, ${\beta _a}$ = −1.8×10⁻³ mmHg·s²·mL⁻², ${\beta _b}$ = −1.2×10⁻⁵ mmHg·s·mL⁻¹·rpm⁻¹ and ${\beta _c}$ = 7.3×10⁻⁶ mmHg·rpm⁻² are derived from a centrifugal pump HeartMate III. Also, corresponding parameters for the resistance and inertance of the outflow cannulae are ${R_i}$ = 0.0677 mmHg·s·mL⁻¹, and ${L_i}$ = 0.0127 mmHg·s²·mL⁻¹. HeartMate III is directly connected to the left ventricular apex and it does not possess inflow cannulae. Therefore, in this study, the parameters for the resistance and inertance of the inflow cannulae are set to 0.

A nonlinear time-varying resistor ${R_k}$ is inserted between the inflow cannulae and the left ventricle to simulate the potential suction event ^[35].

where $\alpha$ is weight parameter and its value is taken as −3.5 s/mL, $\bar P$ is threshold pressure, and its value is taken as 1 mmHg.

2.4.   Simulation protocol

In this study, most parameters under normal conditions in the CLP model and the baroreflex model have values taken from the published literature. The parameters in the CLP model were assigned the values shown in Table A.1 based on the parameter selection in ^[7,14,31]. The parameters in the baroreflex model were set as shown in Table A.2 according to Ref. ^[24]. Table A.1 and Table A.2 were listed in Appendix A.

Compared with normal conditions, heart failure can cause a series of changes, including changes in the cardiovascular system and the baroreflex system. Dilated cardiomyopathy (DCM) is the most common type of nonischemic cardiomyopathy and a major cause of congestive heart failure. The hemodynamic characteristics of patients with DCM usually show some changes, including decreased ventricular end-systolic elastance, increased heart rate, increased systemic vascular resistance, decreased systemic vascular compliance, and increased pulmonary vascular resistance ^[9]. Meanwhile, in the baroreflex system, the baroreflex sensitivity and baroreflex gain for heart failure decreased ^[36,37]. These qualitative and quantitative changes can be mapped into the heart failure model. In the heart failure model, initial changes in left ventricular end-systolic elastance, heart rate, and systemic vascular resistance were simulated by decreasing ${E_{\max, lv0}}$, ${T_0}$, and increasing ${R_{svr0}}$. Besides, increases in ${C_{sar}}$, ${R_{pvr}}$, and ${R_{pv}}$ represented increased systemic arterial compliance and pulmonary circulation resistance. A decrease in ${1 / {{k_a}}}$ indicated a decrease in baroreflex sensitivity. The gain constants (${G_{{R_{svr}}}}$, ${G_{{E_{\max, \, lv}}}}$, ${G_{{E_{\max, \, rv}}}}$, ${G_{{T_s}}}$, and ${G_{{T_v}}}$) were reduced to simulate the change in the baroreflex gain. The parameters were adjusted to simulate the pathological condition of LV failure, and it was set as shown in Table A.3. Table A.3 was listed in Appendix A.

In this study, a series of hemodynamic responses were conducted under normal condition (no LV failure or LVAD support), LV failure condition (no LVAD support) and LV failure condition with different levels of LVAD support (2200 rpm, 3700 rpm, 4300 rpm). In addition, some variables were compared under different conditions, including pressure, flow rate, volume, heart rate, etc. The area of ventricular pressure-volume (P-V) loop can usually be used to measure ventricular stroke work ^[21]. Consequently, the ventricular P-V loop can provide information about the unloading of the ventricle. The P-V loops under different conditions were also compared in this study.

With different levels of pump support, the LVAD operating condition can be evaluated by changes in pump flow rate and physiological variables. Firstly, the HeartMate III pump speed increased linearly from 0 rpm to 5500 rpm within 60 s. The pump speed range, without pump regurgitation and ventricular suction, could be preliminarily determined by the information of the instantaneous pump flow rate during the linear increase in pump speed. The preliminary speed range was obtained by an instantaneous increase in speed, and it may have errors compared to the stable result with constant speed support. Thus, it is necessary to analyze the hemodynamic variables under constant speed support. Then, the constant pump speeds with 100 rpm interval were selected in the preliminary pump speed range, and the stable values of hemodynamic variables corresponding to constant speeds were calculated.

According to the numerical model described above, a program was developed in MATLAB 2010b (The MathWorks Inc., USA) to perform different simulation experiments. The ode45 solver was used to solve time-varying differential equations. In the system of differential equations, the variable values obtained in the last ode45 calculation was used as the initial value in the next ode45 calculation, and the calculation was repeated until the set time was reached. The time step was set to 0.001 second. All simulation times were set to 60 seconds. The system reaches a stable periodic solution after about 15 seconds. The stable result from the 58th to 60th seconds was chosen in all simulated cases for analysis.

3.   Results

3.1.   Hemodynamic response under normal condition

Figure 3 depicts the changes in pressure, flow, and volume under normal condition.

The heart rate was about 68 bpm (beat per minute, bpm).

The left ventricular systolic pressure was 120 mmHg, the aortic pressure changed between 73 and 118 mmHg, and the systemic venous pressure was about 7 mmHg. The right ventricular pressure changed between 3 and 27 mmHg, pulmonary artery changed between 10 and 27 mmHg, and pulmonary venous pressure was about 7 mmHg. The pressures of the left and right atria fluctuated slightly, ranging from 4 to 7 mmHg and 3 to 7 mmHg, respectively.

The left or right ventricle had a stroke volume of about 81 mL and a cardiac output of about 5.5 L/min.

The peak flow rate of the aorta or pulmonary artery was close to 700 mL/s.

These data and waveforms are all consistent with hemodynamic information in the Textbook of Medical Physiology ^[38,39].

3.2.   Hemodynamic response under LV failure condition

The hemodynamic response under the LV failure condition is given in Figure 4.

The heart rate was adjusted to 91 bpm. LV failure caused an increase in heart rate.

The aortic pressure range was reduced from the normal value of 73–118 mmHg to the current 53−77 mmHg, the value of aortic pulse pressure (the difference between the maximum systolic pressure and minimum diastolic pressure) was reduced from the normal value of 45 mmHg to the current value of 24 mmHg, and the left ventricular systolic pressure was reduced from the normal value of 120 mmHg to the current value of 78 mmHg. The pressures of the left atrium and pulmonary vein were elevated to about 10 mmHg, almost increased 43% from the normal value of about 7 mmHg. When the LV failure develops further and the pressure of the pulmonary vein is elevated to above 18 mmHg, pulmonary edema will manifest itself and serious consequence of pulmonary failure will follow ^[39].

The peak flow rate of the aorta or pulmonary artery was reduced about 60% of the normal healthy values.

It is observed that there are prominent changes to all four-chamber volumes in the heart. The current LV volume varied between 179 and 211 mL, greatly elevated from the normal range of 57−138 mL. The stroke volume decreased from the normal value of 81 mL to the current value of 32 mL. The ejection fraction (the ratio of stroke volume to end-diastolic volume) is usually not less than 50%, but the current ejection fraction of the left ventricle was about 15%. The current value of cardiac output in the left ventricle (LVCO) was 2.9 L/min, which was reduced 47% from the normal value of 5.5 L/min. The left atrial volume was also increased from the normal range of 28−54 mL to the current range of 43−65 mL. Accompanying the dilation in the left heart, there was a great volume reduction in the right ventricle and atrium. The right ventricle volume was reduced from the normal range of 39−120 mL to the current range of 26−58 mL. The right atrial volume was also reduced from the normal range of 22−47 mL to the current range of 12−26 mL.

3.3.   Hemodynamic response under LVAD support

Figure 5 shows the changes in pressure, flow rate, and volume under HeartMate III support at low (2200 rpm), medium (3700 rpm), and high (4300 rpm) speeds.

Under LV failure condition with pump support at 4300 rpm, left ventricular pressure approached 0 mmHg and the pressure and flow rate in aorta almost lost its pulsation. Further increased in the support level could cause a collapse in the left ventricle.

At 2200 rpm pump operating speed, the flow rates in the rotary pump and aorta had the intermittent regurgitation (the negative flow rate in numerical model represents regurgitation) and the aortic peak flow rate was larger than that of rotary pump over a heart period. It suggested that the total aortic flow rate was provided by both the left ventricle and the rotary pump.

At 3700 rpm or 4300 rpm, the total aortic flow rate and the pump flow rate completely coincided and the aortic pressure was greater than the left ventricular pressure over a heart period. It indicated that the current aortic valve was closed and the total aortic flow rate was completely provided by the rotary pump.

Figure 6 gives the left and right ventricular P-V loops under normal condition, LV failure condition, and different levels of pump support.

In contrast to the P-V loop in normal condition, the left ventricular P-V loop shifted to the right due to left ventricular end-diastolic volume dilation in the LV failure condition (Figure 6a). At the same time, the area of the left ventricular P-V loop decreased because of the decreased left ventricular end-systolic pressure and stroke volume. In addition, the right ventricular P-V loop shifted to left and its area decreased in LV failure condition relative to that in normal condition (Figure 6b). Compared with the P-V loop in LV failure condition, the left ventricular P-V loop shifted to the right and its area (stroke work) increased due to increased left ventricular end-diastolic volume and end-systolic pressure at 2200 rpm operating speed. Moreover, at 3700 rpm or 4300 rpm pump operating speed, the left ventricular P-V loop shifted to left and its area decreased, even near to 0 (Figure 6a). With the increase of pump support level, the right ventricular P-V loop continuously shifted to the right and its area (stroke work) was also increasing relative to that in LV failure condition (Figure 6b). And the area of the right ventricular P-V loop at 3700 rpm was relatively closer to that under normal condition.

Ranges of pressures and volumes in left and right heart, arterial and venous pressures, cardiac output, ventricular end-systolic elastances, systemic vascular resistance, and heart rate are shown in Table 1 under different conditions. Under LVAD support, the cardiac output in the left heart was provided by both left ventricle and rotary pump.

Compared with the LV failure condition, the ranges of pressures and volumes in the left or right heart increased at 2200 rpm. But the ranges of pressures and volumes in the left heart decreased and them in the right heart increased, at 3700 rpm or 4300 rpm. Furthermore, both arterial pressure and cardiac output increased with the LVAD support. At 3700 rpm, the mean aortic pressure (MAP) was 100 mmHg and cardiac output in the left ventricle (LVCO) or right ventricle (RVCO) was 6 L/min, which was closer to that in normal condition. With the increase of the LVAD support level, the values of left or right ventricular end-systolic elastance, systemic vascular resistance, and heart rate decreased. At 3700 rpm LVAD operating speed, the heart rate was restored from 91 bpm in LV failure condition to the current 72 bpm.

3.4.   Effect of pump speed on hemodynamic variables

Figure 7 shows the instantaneous pump flow rate during a linear increase in pump speed from 0 rpm to 5500 rpm. The pump flow rate had instantaneous or persistent regurgitation when the pump speed was between 0 and 3594 rpm. And the pump regurgitation disappeared when pump speed was higher than 3594 rpm. However, the amplitude of the pump flow rate signal suddenly increased when the speed reached 4712 rpm. It indicated that the left ventricle had an onset of suction. Therefore, the preliminary pump speed range should be between 3594 rpm and 4712 rpm.

For convenience of calculation, the constant speeds were varied from 3600 rpm to 4700 rpm with 100 rpm interval, and the corresponding changes in hemodynamic variables are shown in Figure 8.

With the increase of pump speed, the pressures and volumes in the left heart decreased gradually. On the contrary, the pressures and volumes in the right heart gradually increased. Meanwhile, aortic pressure and pulmonary artery pressure were increasing. Moreover, the pulsatility of aortic pressure gradually weakened, and that of pulmonary artery pressure gradually enhanced. The increase of speed also enhanced the cardiac output and the minimum pump flow rate over a heart period. There was no negative value of the minimum pump flow rate within 3600−4700 rpm, which indicated that pump regurgitation disappeared when the speed was higher than 3600 rpm. However, when the pump operating speed is higher than 4400 rpm, the minimum pump flow rate decreased, the aortic pulse pressure increased, the left ventricular pressure approached 0 mmHg, and the left ventricular volume was less than 50 mL, which implied that the left ventricular suction occurred.

4.   Discussion

In this study, a CLP model for describing the circulatory loop was developed to simulate hemodynamic responses. The CLP model was coupled with the baroreflex model and LVAD model, which provided a platform for the simulation of hemodynamics and the evaluation of LVAD operating condition.

LVADs have become the standard treatment for advanced heart failure. Nevertheless, pump regurgitation may occur under pump support at low speeds. Inadequate support may cause mean aortic pressure and cardiac output to be lower than the normal level, which is difficult to meet the patients' normal physiological perfusion. Pump support at high speed allows the left ventricle to unload, but it also increases right ventricular load (Figure 6). However, long-term overloading for the right ventricle may cause right heart failure ^[40].

The HeartMate III operating speed should be in the range of 3600−4400 rpm with the selected parameter settings to avoid pump regurgitation and ventricular suction. Generally, an optimal speed range is set for clinical continuous-flow LVADs, and it can be determined by monitoring cardiac index (cardiac output) and left ventricular size based on hemodynamic data ^[40]. In this study, cardiac output in the left ventricle and left ventricular end-diastolic volume are 5.5 L/min and 138 mL in normal conditions, respectively. In Figure 8, the cardiac output of 5.6 L/min at 3600 rpm is near to the normal cardiac output. And at 3800 rpm, the left ventricular end-systolic volume is 139 mL that is closer to the normal value of 138 mL. Hence, the low-end of the optimal speed range can be set to 3600 rpm and the high-end can be set to 3800 rpm.

In the physiological or pathological cardiovascular system, the baroreflex model objectively describes the changes in ventricular end-systolic elastances, systemic vascular resistance, and heart rate. Different levels of LVAD support also affect ventricular end-systolic elastances, systemic vascular resistance, and heart rate. However, heart rate is usually maintained as a constant in cardiovascular models without baroreflex regulation, which may not be in line with the actual situation.

A pulse pressure of greater than 15 mmHg is considered to be pulsatile flow while a pulse pressure of less than 15 mmHg is considered to be nonpulsatile flow ^[41]. In this study, the aortic pulse pressure is less than 15 mmHg for the operating speed between 3600 rpm and 4400 rpm in Figure 8. The continuous flow LVADs usually operate at a constant speed to provide the patient with adequate blood flow. However, long-term LVAD constant support diminishes pulse pressure and vascular pulsatility, which can lead to late complications such as aortic insufficiency, gastrointestinal bleeding, and hemorrhagic strokes ^[31,42]. Recently, the new generation of commercial products, such as HVAD (Medtronic, Inc., Minneapolis, MN, USA) and HeartMate III (Abbott Inc., Abbott Park, IL, USA), come with an artificial pulse. The artificial pulse is formed by intermittent speed change, and its main goal is to reduce the rate of pump thrombosis. Importantly, the cycle of variable speed is asynchronous with the native heart period. Moreover, the aortic pressure waveform becomes non-physiological ^[43]. Therefore, certain issues should be studied and solved, including modulation and optimization of pulsating speed with heart rate synchronization, hemodynamic analysis at pulsating speed support, and design of the physiological controller. These will be the subject of further research.

The limitations of this study are as follows.

From a hemodynamic viewpoint, the numerical model can reflect clinical observations under LV failure and LVAD support conditions. The numerical model cannot replace the importance and significance of in vivo research, and cannot reproduce all the expected clinical responses (such as neurohumoral response, myocardial remodeling). But it can prove the feasibility of some concepts to provide a valuable initial step for further in vivo research.

In this study, the valve model is assumed that valves can be opened and closed immediately, and it ignores the process of valve closure. Thus, the ideal valve model cannot simulate the valve flow rate during valve closure. The non-ideal diode modeling for the heart valves is a feasible way to simulate the dynamic process of valves ^[44,45,46].

The experimental data of HeartMate III are from earlier versions of the product design ^[4]. Since the pump data was first published, the HeartMate III design has gone through several versions. However, the method in this study is still valid.

In this study, the regulation of venous volume in the baroreflex model is ignored. This is because the systemic venous compliance is configured to be constant. But this does not affect the change in systemic venous pressure.

5.   Conclusion

In this study, a numerical model had been developed to simulate the hemodynamics and to evaluate the operating condition of continuous-flow left ventricular assist device. The hemodynamic responses under various conditions (normal, left ventricular failure, and different levels of pump support) had been simulated successfully, and the hemodynamic variables under different conditions had been compared. The simulation results had shown that the developed numerical model provides a feasible simulation platform for reproducing hemodynamics. Also, the centrifugal pump HeartMate III had been taken as an example of continuous-flow LVAD, and its operating condition had been evaluated and discussed. The simulation results had shown that in the simulation system, HeartMate III is safe to operate between 3600 rpm and 4400 rpm, and that pump operating speeds are adjusted between 3600 rpm and 3800 rpm to provide optimal physiological perfusion. What makes this paper regret is that some system parameters are quoted or set as constants, it gives rise to some limitations on the final results. Actually, these parameters are variable because of the interpatient and intrapatient variations. The future work will simulate combining the patient data, and obtain the hemodynamic response for variable physiological or pathological parameters.

Acknowledgments

The study is supported by the Key Technology Research and Development Program of Shandong (CN) (No. 2017GSF221009).

Conflict of interest

The authors declare that there is no conflict of interest.

Appendix A