Computational Simulation of Cardiac Function and Blood Flow in the Circulatory System under Continuous Flow Left Ventricular Assist Device Support during Atrial Fibrillation

Prevalence of atrial fibrillation (AF) is high in heart failure patients supported by a continuous flow left ventricular assist device (CF-LVAD); however, the long term effects remain unclear. In this study, a computational model simulating effects of AF on cardiac function and blood flow for heart failure and CF-LVAD support is presented. The computational model describes left and right heart, systemic and pulmonary circulations and cerebral circulation, and utilises patient-derived RR interval series for normal sinus rhythm (SR). Moreover, AF was simulated using patient-derived unimodal and bimodal distributed RR interval series and patient specific left ventricular systolic functions. The cardiovascular system model simulated clinically-observed haemodynamic outcomes under CF-LVAD support during AF, such as reduced right ventricular ejection fraction and elevated systolic pulmonary arterial pressure. Moreover, relatively high aortic peak pressures and middle arterial peak flow rates during AF with bimodal RR interval distribution, reduced to similar levels as during normal SR and AF with unimodal RR interval distribution under CF-LVAD support. The simulation results suggest that factors such as distribution of RR intervals and systolic left ventricular function may influence haemodynamic outcome of CF-LVAD support during AF.


Introduction
Atrial fibrillation (AF) is associated with an increased risk of stroke, heart failure and cognitive impairment, and it also increases the risk of death in heart failure patients [1][2][3][4]. The distribution of RR interval frequencies in AF patients is mainly unimodal or bimodal [5]. Bimodal RR interval distribution can be a sign of dual atrioventricular nodal pathway physiology [6]. AF and heart failure may coexist and show bidirectional interactions, with heart failure promoting AF and AF with an uncontrolled ventricular rate leading to ventricular tachycardiomyopathy [7,8]. Moreover, preceding and pre-preceding RR interval durations seem to play a role in the left ventricular systolic function in patients with chronic AF and dilated cardiomyopathy [9].
The prevalence of AF is high in heart failure patients supported by a continuous flow left ventricular assist device (CF-LVAD) and it reaches up to 46% [10]. However, the effects of AF remain unclear in heart failure patients implanted with a CF-LVAD for long-term support, and there is no consensus on the risks associated with AF during heart pump support [11]. There are studies reporting that AF is not associated with risk of death, stroke or thromboembolism during CF-LVAD support [12][13][14][15]. LVADs may even induce favourable atrial structural and electrical remodelling in the patients with AF; however, they may also be associated with increased mortality rates and risk of thromboembolic events after CF-LVAD implantation [10,16,17], contradictory to the previous studies.

RR Interval Series
In this study, RR intervals were extracted from the patients' Electrocardiogram (ECG) data sets available on PhysioNet (https://physionet.org/). These were used as the duration of consecutive heartbeats to simulate AF and normal SR [25] over a ten minute period. The distribution of RR intervals during AF can be unimodal, bimodal or multimodal [5] depending on the average heartbeat duration over a certain period [26]. Unimodal and bimodal RR interval distributions for AF were taken from a long-term AF database, as given in [27]. The data in this database was digitised at a sampling rate of 128 Hz at 16 bit resolution. The signals were band-pass filtered utilising a cut-off frequency at 1 and 50 Hz to avoid power line interference [27]. RR interval data from three different subjects was used in the simulations. The mean of the RR intervals simulating normal sinus rhythm (SR) was 0.828 s with a 0.051 s standard deviation and 0.006 coefficient of variation. The mean of the RR intervals simulating AF with a unimodal distribution was 0.512 s with a 0.106 s standard deviation and 0.207 coefficient of variation, whereas the mean RR interval duration simulating AF with a bimodal distribution was 0.884 s with a 0.260 s standard deviation and 0.294 coefficient of variation in the cardiovascular system model. The RR intervals and frequency histograms of the normal SR and AF data used in the simulations is given in Figure 1.

Cardiac Function
Ventricular functions were described using the relationship between ventricular pressure and volume signals by utilising both active and passive contraction components. The left ventricular pressure signal (plv) during the active contraction of the left ventricle is driven by an activation function (fact,lv) and it utilises the end-systolic elastance (Ees,lv), left ventricular volume (Vlv) and left ventricular zero pressure volume (Vlv,0) [28]. The left ventricular passive pressure component (plv,p) was modelled using an exponential relationship between pressure and volume signals [28].
( ) lv,a es,lv lv lv,0 act,lv Left ventricular volume (Vlv) was described using the left ventricular radius (rlv), long axis length (llv) and an additional coefficient (Klv), which includes effects of the contraction in the long axis and scales the proportion between the left ventricular radius and volume over a cardiac cycle [28]. Change of the left ventricular radius (rlv) was described utilising the flow rates through the aortic and mitral valves (Qav, Qmv), left ventricular volume, long axis length and the coefficient Klv [28].
The activation function (fact,lv) was described using instantaneous time (t), the durations of active contraction (T1) and relaxation (T2) phases and the duration of an RR interval (RR) [29].

Cardiac Function
Ventricular functions were described using the relationship between ventricular pressure and volume signals by utilising both active and passive contraction components. The left ventricular pressure signal (p lv ) during the active contraction of the left ventricle is driven by an activation function (f act,lv ) and it utilises the end-systolic elastance (E es,lv ), left ventricular volume (V lv ) and left ventricular zero pressure volume (V lv,0 ) [28]. The left ventricular passive pressure component (p lv,p ) was modelled using an exponential relationship between pressure and volume signals [28].
Left ventricular volume (V lv ) was described using the left ventricular radius (r lv ), long axis length (l lv ) and an additional coefficient (K lv ), which includes effects of the contraction in the long axis and scales the proportion between the left ventricular radius and volume over a cardiac cycle [28]. Change of the left ventricular radius (r lv ) was described utilising the flow rates through the aortic and mitral valves (Q av , Q mv ), left ventricular volume, long axis length and the coefficient K lv [28].
The activation function (f act,lv ) was described using instantaneous time (t), the durations of active contraction (T 1 ) and relaxation (T 2 ) phases and the duration of an RR interval (RR) [29].
The durations of the active contraction (T 1 ) and relaxation (T 2 ) phases in the cardiovascular system model with normal SR were described by utilising linear relationships considering the physiological data [30].
The left atrial pressure (p la ) and volume (V la ) relationship was described using elastance (E la ) [28,29].
Left atrial volume (V la ) was described using the left atrial radius (r la ), long axis length (l la ) and an additional coefficient (K la ) [28]. Change of the left atrial radius (r la ) was described utilising the flow rates through the mitral valve and pulmonary vein (Q mv , Q vp ), left atrial volume, long axis length and the coefficient K la [28].
dr la dt = 3 Q pv − Q mv 4πK la l la 3V la 2πK la l la E la (t) = E min,la + 0.5(E max,la − E min,la )f act.la (t − D), Here, E min,la and E max,la are the minimal and maximal atrial elastances over a cardiac cycle, respectively, and T a is the atrial activation time instant. The right ventricular pressure and volume signals were modelled in a similar way to the left ventricular pressure and volume signals using different parameter values. The right atrial pressure and volume signals were modelled in a similar way to the left atrial pressure and volume signals using different parameter values. Detailed information about the model used to describe atrial ventricular functions can be found in [28].
The atrial contraction contributes to ventricular filling around 10%-15% over each cardiac cycle for a healthy condition [31]. However, the absence of P-wave indicates an impaired function of atrial contraction due to continuous atrial activity resulting from both ectopic and re-entrant sources [32]. The impaired atrial contraction was described in the cardiovascular system models simulating AF using only the minimal atrial elastance (E min,la ) as given below [33].
The systolic elastance in the cardiovascular system models, simulating heart failure with normal SR and AF with unimodal and bimodal RR interval distributions, was adjusted to provide the similar mean blood flow from the left ventricle over a 10 min period. Typical value of systolic elastance in healthy subjects is around 2.5 mmHg/mL; however, it reduces significantly during heart failure [34]. Therefore, systolic elastance in the cardiovascular system model, simulating the heart failure with normal SR, was set to 0.865 mmHg/mL to simulate reduced ejection fraction. A positive correlation exists between the systolic function of left ventricle and the ratio of the preceding RR intervals in heart failure patients with AF [9]. The left ventricular end-systolic elastances in the cardiovascular system models, simulating AF with unimodal (E es,lv,u ) and bimodal (E es,lv,b ) RR interval distributions, were adjusted by adopting patient-specific correlations from heart failure patients, as given in [9], to provide the similar mean blood flows in these models as in the cardiovascular system model simulating heart failure with normal SR. The left ventricular end-systolic elastances in the cardiovascular system models simulating AF were described as below.

Circulatory System
The circulatory system model includes the systemic and pulmonary circulation and cerebral circulation. The model also simulates blood flow in the circle of Willis, as presented in [35]. Blood flow in the circulatory system model was described using a lumped parameter model, which included electrical analogues for resistance (R), compliance (C) and inertia (L). The heart valves were modelled as ideal diodes allowing for one-way blood flow. The aortic blood pressure (p ao ) and flow rate signals (Q ao ), as well as the aortic valve flow signal (Q av ) in the circulatory system model, have been given below. dp ao Here, p aa represents pressure in the aortic arch, C ao , R ao and L ao are the compliance, resistance and inertance in the aorta, respectively, and R av is the aortic valve resistance. The other compartments in the circulatory system were modelled in the same way using different parameter values.

Continuous Flow Left Ventricular Assist Device (CF-LVAD) Support
To simulate CF-LVAD support, a model that estimates the pressure difference across an axial CF-LVAD (∆p CF-LVAD ), considering operating speed of the pump (ω CF-LVAD ), flow rate (Q CF-LVAD ) and change of the flow rate through the pump (dQ CF-LVAD /dt) [36], was integrated into the cardiovascular system model.
In the equations above, L CF-LVAD and R CF-LVAD are the inertance and resistance effects in the pump, respectively. K, k 1 and k 2 are the estimated parameters [36]. The CF-LVAD was operated at a 10,500 rpm rotation speed in the cardiovascular system models simulating the heart pump support. The CF-LVAD model was implemented into the cardiovascular system model by modifying the equations which describe change of the left ventricular radius (Equation (5)) and aortic pressure (Equation (19)) as below.
Systemic arteriolar resistance was 1.25 mmHg/mL for the cardiovascular system model simulating heart failure as it increases to maintain perfusion levels in the body for a failing heart [37]. It was reduced to a healthy value, 0.75 mmHg/mL, for CF-LVAD support in the simulations, because mechanical circulatory support restores the blood flow, and pressure levels in circulatory system baroreflex control is expected to reduce the peripheral systemic resistance with increasing arterial pressure [38].
The electric-analogue of the heart chambers and the circulatory loop used to simulate normal SR and AF with unimodal and bimodal RR interval distributions for heart failure and CF-LVAD support is given in Figure 2. The list of abbreviations and the parameter values used in the circulatory loop and cardiac functions is given in the Appendix A section.
Systemic arteriolar resistance was 1.25 mmHg/mL for the cardiovascular system model simulating heart failure as it increases to maintain perfusion levels in the body for a failing heart [37]. It was reduced to a healthy value, 0.75 mmHg/mL, for CF-LVAD support in the simulations, because mechanical circulatory support restores the blood flow, and pressure levels in circulatory system baroreflex control is expected to reduce the peripheral systemic resistance with increasing arterial pressure [38].
The electric-analogue of the heart chambers and the circulatory loop used to simulate normal SR and AF with unimodal and bimodal RR interval distributions for heart failure and CF-LVAD support is given in Figure 2. The list of abbreviations and the parameter values used in the circulatory loop and cardiac functions is given in the Appendix section. The simulations were performed using MATLAB Simulink R2017a. The set of equations was solved using the ode15 s solver. The maximum step size was 1e-3 s and the relative tolerance was set to 1e-3. All of the cardiovascular system models were run at a 40 s simulation time at a 75 bpm constant heart rate, to obtain a periodic steady state solution, before the RR intervals were used as temporal ranges for the heartbeats. The simulations were performed using MATLAB Simulink R2017a. The set of equations was solved using the ode15 s solver. The maximum step size was 1e-3 s and the relative tolerance was set to 1e-3. All of the cardiovascular system models were run at a 40 s simulation time at a 75 bpm constant heart rate, to obtain a periodic steady state solution, before the RR intervals were used as temporal ranges for the heartbeats.

Results
The mean values and standard deviations of maximal and minimal volumes in the heart chambers, the left and right ventricular ejection fractions, the mean aortic and pulmonary arterial pressures and middle cerebral arterial flow for normal SR and AF with unimodal and bimodal RR interval distributions over each cardiac cycle, and cardiac and mean pump outputs for heart failure and CF-LVAD support over a ten min period are given in Table 1. Table 1. The mean values and standard deviations of maximal and minimal volumes in the heart chambers (V la,max , V la,min , V ra,max , V ra,min , V lv,ed , V lv,es , V rv,ed , V rv,es ), left and right ventricular ejection fractions (EF Vlv , EF Vrv ), mean aortic and pulmonary arterial pressures and middle cerebral arterial flow rate (p ao,mean , p ap,mean , Q mca,mean ) for normal sinus rhythm (NSR), and atrial fibrillation with unimodal and bimodal RR interval distributions (UAF, BAF) over each cardiac cycle, cardiac and mean pump outputs (CO, MPOs) over a ten minute period for heart failure and CF-LVAD support. The waveforms of the left and right atrial and ventricular volumes for the simulated normal SR and AF with unimodal and bimodal RR interval distributions during heart failure and CF-LVAD support over a 6 s period are given in Figure 3. Amplitude of the left and right atrial volume signals was relatively high during AF with bimodal RR interval distribution. Similarly, amplitude of the left and right ventricular volume signals during AF with bimodal RR interval distribution was higher, with respect to the amplitude of left and right ventricular volume signals during normal SR and AF with unimodal RR interval distribution. CF-LVAD support reduced the left atrial and ventricular volumes, whereas right atrial volumes and right ventricular end-diastolic volumes increased under CF-LVAD support for all the simulated cases. The waveforms for the left and right atrial and ventricular pressure signals for the simulated normal SR and AF with unimodal and bimodal RR interval distributions for heart failure and CF-LVAD support over a 6 s period are given in Figure 4. Amplitude of the left and right atrial volume signals was relatively high during AF with bimodal RR interval distribution. Similarly, amplitude of the left and right ventricular volume signals during AF with bimodal RR interval distribution was higher, with respect to the amplitude of left and right ventricular volume signals during normal SR and AF with unimodal RR interval distribution. CF-LVAD support reduced the left atrial and ventricular volumes, whereas right atrial volumes and right ventricular end-diastolic volumes increased under CF-LVAD support for all the simulated cases. The waveforms for the left and right atrial and ventricular pressure signals for the simulated normal SR and AF with unimodal and bimodal RR interval distributions for heart failure and CF-LVAD support over a 6 s period are given in Figure 4. Left atrial pressure during AF with bimodal RR interval distribution remained relatively low for heart failure. CF-LVAD support reduced the overall left atrial pressures for all the simulated cases, whereas right atrial pressures increased under CF-LVAD support. Left ventricular peak pressure was relatively high over some of the cardiac cycles during AF with bimodal RR interval distribution for heart failure, because of the simulated relation between left ventricular systolic elastance and preceding and pre-preceding cardiac cycle duration ratio. CF-LVAD support reduced the left ventricular pressures for all the simulated cases. Right ventricular pressures slightly decreased under CF-LVAD support. The waveforms for the aortic and pulmonary arterial pressure signals and left middle cerebral arterial blood flow rate signal for the simulated normal SR and AF with unimodal Left atrial pressure during AF with bimodal RR interval distribution remained relatively low for heart failure. CF-LVAD support reduced the overall left atrial pressures for all the simulated cases, whereas right atrial pressures increased under CF-LVAD support. Left ventricular peak pressure was relatively high over some of the cardiac cycles during AF with bimodal RR interval distribution for heart failure, because of the simulated relation between left ventricular systolic elastance and preceding and pre-preceding cardiac cycle duration ratio. CF-LVAD support reduced the left ventricular pressures for all the simulated cases. Right ventricular pressures slightly decreased under CF-LVAD support. The waveforms for the aortic and pulmonary arterial pressure signals and left middle cerebral arterial blood flow rate signal for the simulated normal SR and AF with unimodal and bimodal RR interval distributions for heart failure and CF-LVAD support over a 6 s period have been given in Figure 5.

NSR
Appl. Sci. 2020, 10, 876 10 of 16 and bimodal RR interval distributions for heart failure and CF-LVAD support over a 6 s period have been given in Figure 5. Amplitude of the aortic pressure signal during AF with bimodal RR interval distribution was higher with respect to amplitudes of the aortic pressure signals for the other cases. CF-LVAD support reduced the systolic aortic pressure for all the cases while increasing the diastolic pressure. Moreover, amplitudes of the aortic pressure signals were similar under CF-LVAD support for all the cases. Pulmonary arterial pressure reduced under CF-LVAD support for all the simulated cases. The amplitude of the middle cerebral arterial flow rate signal during AF with bimodal RR interval distribution was higher with respect to the amplitudes of the middle cerebral arterial flow rate during normal SR and AF with unimodal RR interval distribution. The amplitudes of the middle cerebral arterial flow rate signals were similar for all cases where there was CF-LVAD support.

Discussion
This study utilises a computational model simulating cardiovascular system dynamics and an integrated CF-LVAD model. The cardiovascular system model allows to simulate different physiological scenarios and mechanisms. For instance, the interaction between the preceding and pre-preceding RR intervals and left ventricular systolic function in patients with chronic AF and dilated cardiomyopathy was simulated in this study. Therefore, the effect of RR interval distributions Amplitude of the aortic pressure signal during AF with bimodal RR interval distribution was higher with respect to amplitudes of the aortic pressure signals for the other cases. CF-LVAD support reduced the systolic aortic pressure for all the cases while increasing the diastolic pressure. Moreover, amplitudes of the aortic pressure signals were similar under CF-LVAD support for all the cases. Pulmonary arterial pressure reduced under CF-LVAD support for all the simulated cases. The amplitude of the middle cerebral arterial flow rate signal during AF with bimodal RR interval distribution was higher with respect to the amplitudes of the middle cerebral arterial flow rate during normal SR and AF with unimodal RR interval distribution. The amplitudes of the middle cerebral arterial flow rate signals were similar for all cases where there was CF-LVAD support.

Discussion
This study utilises a computational model simulating cardiovascular system dynamics and an integrated CF-LVAD model. The cardiovascular system model allows to simulate different physiological scenarios and mechanisms. For instance, the interaction between the preceding and pre-preceding RR intervals and left ventricular systolic function in patients with chronic AF and dilated cardiomyopathy was simulated in this study. Therefore, the effect of RR interval distributions on the haemodynamic outcome CF-LVAD support can be simulated using the presented cardiovascular system model.
The mean RR interval duration in the cardiovascular system model simulating AF with unimodal RR interval distribution was 0.512 s corresponding to a 117 bpm heart rate. Although such a high heart rate is common in untreated AF patients, the current guidelines suggest to reduce the heart rate below 100 bpm in AF patients [39]. The RR interval duration in the cardiovascular system model simulating AF with bimodal RR interval distribution was 0.884 s corresponding to 68 bpm. This is below the suggested target heart rate. It should also be noted that the RR interval durations used in this study are patient specific and the AF treatment is not included in the utilised numerical models.
The left ventricular ejection fraction was relatively low in the cardiovascular system simulating AF with unimodal RR interval distribution. It has already been shown that there is a strong correlation between the ejection fraction and the ratio of preceding and pre-preceding RR intervals during AF [40]. A higher preceding and pre-preceding RR interval ratio results in high ejection fraction as in the cardiovascular system model simulating AF with bimodal RR interval distribution [40].
The low right ventricular ejection fraction during normal SR and AF with bimodal RR interval distribution increased above 50% under CF-LVAD support similar to a healthy cardiovascular system [41]. Although the right ventricular ejection fraction during AF with a unimodal RR interval distribution increased as well under CF-LVAD support, it was around 40% which indicates an impaired function. It is already known that the right ventricular function may be impaired in heart failure patients with AF under CF-LVAD support [42]. The cardiovascular system model produced impaired right ventricular function for relatively short cardiac cycles and unimodal RR interval distribution. Systolic pulmonary arterial pressures also reported to be relatively high in LVAD implanted heart failure patients with AF with respect to patients without AF [20]. Again, the cardiovascular system model produced relatively high systolic pulmonary arterial pressures for unimodal RR interval distribution under CF-LVAD support (Table 1). Higher aortic peak pressures and middle cerebral arterial peak flow rates may result in increased cerebral hypertensive events [24]. Relatively high aortic peak pressures and middle cerebral arterial peak flow rates were relatively high produced for bimodal RR interval distribution without CF-LVAD support ( Figure 5). The peak values in the aortic pressure and middle cerebral arterial flow rate signals were reduced to similar values under CF-LVAD support for all the cases ( Figure 5). Such a result may show the beneficial effect of CF-LVAD therapy on the patients with AF. However, it should be noted that the mean values of the aortic pressures and middle cerebral arterial flow rates were different in each case under heart pump support. Although cardiovascular system model generated similar cardiac outputs for each case, the mean pump outputs were remarkably different. The different mean pump outputs were obtained for each case depending on the duration of the cardiac cycles utilised in the simulations and varying CF-LVAD resistance used to describe CF-LVAD support (Equation (21)).
The left and right atrial volumes are similar in healthy cardiovascular system [43][44][45]. However, relatively high left atrial volumes were simulated in the cardiovascular system models for all of the cases because the simulations were performed considering impaired systolic left ventricular function which may cause dilation in the left atrium [46]. It should be noted that the left atrial volumes were reduced to values comparable to healthy conditions for all the simulated cases under CF-LVAD support.
The simulation results in this study depend on the selected patient-specific relations and parameter settings. However, it should be noted that utilised patient-specific relations and parameter settings were able simulate conditions such as reduced right ventricular ejection fraction or high systolic pulmonary arterial pressure, which are also observed in the clinical studies [20,42] for different cases.
The simulations were performed using data from two different AF patients and one healthy individual. The mean values and variability of the haemodynamic parameters for each cardiac cycle are given in the Results section. A statistical comparison between the groups was not performed because the data from one patient were utilised for each case. Nonetheless, the results provide useful insights about the effect of preceding and pre-preceding RR intervals of AF and CF-LVAD support. The cardiovascular system model utilises patient-derived series to simulate AF. Incorporating dynamical models that simulate electrocardiogram signals or dual pathway atrioventricular nodal conduction [47][48][49][50] with the existing cardiovascular system model may allow to simulate the dependence of the cardiovascular response on the rate and regularity of ventricular activation in AF [51]. Elaborating the cardiovascular system model to simulate atrioventricular nodal dynamics and the evaluation of CF-LVAD support and validation of the results will be future work.

Conclusions
The presented cardiovascular system model simulates AF for unimodal and bimodal RR interval distributions and Normal SR. Clinically-observed haemodynamic outcomes under CF-LVAD support during AF such as reduced right ventricular ejection fraction or high systolic pulmonary arterial pressure were simulated utilising patient-derived RR interval series and the interaction between the preceding and pre-preceding RR interval and left ventricular systolic function. The results in the submitted study show that patient specific conditions such as duration of preceding and pre preceding RR intervals may play a role on the outcome of the CF-LVAD therapy during atrial fibrillation. Moreover, the simulation results provide insights and possible future research directions for AF under CF-LVAD support.
Funding: This research received no external funding.

Conflicts of Interest:
The author declares no conflicts of interest.