Framework for patient-specific simulation of hemodynamics in heart failure with counterpulsation support

Despite being responsible for half of heart failure-related hospitalizations, heart failure with preserved ejection fraction (HFpEF) has limited evidence-based treatment options. Currently, a substantial clinical issue is that the disease etiology is very heterogenous with no patient-specific treatment options. Modeling can provide a framework for evaluating alternative treatment strategies. Counterpulsation strategies have the capacity to improve left ventricular diastolic filling by reducing systolic blood pressure and augmenting the diastolic pressure that drives coronary perfusion. Here, we propose a framework for testing the effectiveness of a soft robotic extra-aortic counterpulsation strategy using a patient-specific closed-loop hemodynamic lumped parameter model of a patient with HFpEF. The soft robotic device prototype was characterized experimentally in a physiologically pressurized (50–150 mmHg) soft silicone vessel and modeled as a combination of a pressure source and a capacitance. The patient-specific model was created using open-source software and validated against hemodynamics obtained by imaging of a patient (male, 87 years, HR = 60 bpm) with HFpEF. The impact of actuation timing on the flows and pressures as well as systolic function was analyzed. Good agreement between the patient-specific model and patient data was achieved with relative errors below 5% in all categories except for the diastolic aortic root pressure and the end systolic volume. The most effective reduction in systolic pressure compared to baseline (147 vs. 141 mmHg) was achieved when actuating 350 ms before systole. In this case, flow splits were preserved, and cardiac output was increased (5.17 vs. 5.34 L/min), resulting in increased blood flow to the coronaries (0.15 vs. 0.16 L/min). Both arterial elastance (0.77 vs. 0.74 mmHg/mL) and stroke work (11.8 vs. 10.6 kJ) were decreased compared to baseline, however left atrial pressure increased (11.2 vs. 11.5 mmHg). A higher actuation pressure is associated with higher systolic pressure reduction and slightly higher coronary flow. The soft robotic device prototype achieves reduced systolic pressure, reduced stroke work, slightly increased coronary perfusion, but increased left atrial pressures in HFpEF patients. In future work, the framework could include additional physiological mechanisms, a larger patient cohort with HFpEF, and testing against clinically used devices.

Despite being responsible for half of heart failure-related hospitalizations, heart failure with preserved ejection fraction (HFpEF) has limited evidence-based treatment options. Currently, a substantial clinical issue is that the disease etiology is very heterogenous with no patient-specific treatment options. Modeling can provide a framework for evaluating alternative treatment strategies. Counterpulsation strategies have the capacity to improve left ventricular diastolic filling by reducing systolic blood pressure and augmenting the diastolic pressure that drives coronary perfusion. Here, we propose a framework for testing the e ectiveness of a soft robotic extra-aortic counterpulsation strategy using a patient-specific closed-loop hemodynamic lumped parameter model of a patient with HFpEF. The soft robotic device prototype was characterized experimentally in a physiologically pressurized ( -mmHg) soft silicone vessel and modeled as a combination of a pressure source and a capacitance. The patient-specific model was created using open-source software and validated against hemodynamics obtained by imaging of a patient (male, years, HR = bpm) with HFpEF. The impact of actuation timing on the flows and pressures as well as systolic function was analyzed. Good agreement between the patient-specific model and patient data was achieved with relative errors below % in all categories except for the diastolic aortic root pressure and the end systolic volume. The most e ective reduction in systolic pressure compared to baseline ( vs. mmHg) was achieved when actuating ms before systole. In this case, flow splits were preserved, and cardiac output was increased ( . vs. . L/min), resulting in increased blood flow to the coronaries ( . vs. . L/min). Both arterial elastance ( . vs. . mmHg/mL) and stroke work ( . vs. . kJ) were decreased compared to baseline, however left atrial pressure increased ( . vs. . mmHg). A higher actuation pressure is associated with higher systolic pressure reduction and slightly higher coronary flow. The soft robotic device prototype achieves reduced systolic pressure, reduced stroke work, slightly

. Introduction
In heart failure (HF) the heart is unable to appropriately supply the receiving organs with blood, leading to over five million hospitalizations yearly (1)(2)(3)(4). To date, we have limited evidence-based treatment options for patients with heart failure with preserved ejection fraction (HFpEF) even though they represent nearly half of the HF population (5). HFpEF is characterized by suboptimal filling of the LV during diastole due to lower compliance of the ventricular walls (6). The syndrome leads to higher LV pressure levels to obtain a certain LV volume than in healthy subjects (7). Consequently, patients experience dyspnea on exertion, reduced quality of life, and reduced life expectancy (6,8).
Recent efforts have proposed multiple device-based treatment options such as the intra-atrial shunt device (9, 10), use of conventional left ventricular assist devices (11,12) for atrial decompression, a valveless pulsatile left ventricular assist device (13), and a left atrial assist device (14). However, the development of treatment strategies for HFpEF is complicated by the heterogeneity in pathophysiologies and associated comorbidities (3,15,16). A recent study compared the above therapies for multiple cohorts of patients with HFpEF (17) and found that each treatment may be effective only in a sub-cohort of HFpEF. As an alternative to grouping patients, we aimed to present a framework for patient-specific testing of treatment strategies and use it to test the usefulness of aortic counterpulsation in HFpEF.
Counterpulsation is a treatment strategy often applied in HF with reduced ejection fraction (HFrEF), but understudied for its usefulness in HFpEF. Counterpulsation reduces afterload and systolic pressure, and increases coronary perfusion (18). Intra-aortic balloon pumps (IABPs) are the most commonly used form of counterpulsation; the device is placed in the descending aorta and directs diastolic blood flow toward the periphery (18). The risk factors of IABPs include vascular injury, embolic complications during deployment or removal, and bacterial infections (19), motivating a need for non-blood contacting extra-aortic actuation. The C-Pulse System (Nuwellis Inc., Eden Prairie, Mn) is an extra-aortic counterpulsation device for HF which has shown to reduce systolic pressure and afterload in clinical trials (20,21). The hemodynamic benefits of counterpulsation, such as decreased systolic pressures and increased coronary filling, have previously been suggested in small cohorts to improve diastolic filling (22,23). Potential benefits of counterpulsation in patients with HFpEF have not been previously assessed.
Patient-specific lumped parameter computational models provide a low-cost option for proof-of concept testing of innovative treatment options in interaction with the cardiovascular system. Lumped parameter computational models are efficient (24) and are particularly suitable for testing in heterogeneous patient populations as each component can be tuned to match anatomical and hemodynamic data in a patient-specific manner (25,26). Given the heterogeneity of the HFpEF population, a patient-specific approach is particularly important to allow an extension to a patient cohort in the future.
We demonstrated and evaluated a framework (Figure 1) for prototyping a soft robotic device prototype and modeling its interaction with patient-specific HFpEF hemodynamics. An extra-aortic soft-robotic counterpulsation device prototype was characterized experimentally, as a part of this work. We assess the impact of this device prototype on lowering peak systolic pressures and increasing coronary filling using a closed-loop patient-specific lumped parameter model at different timings and actuation levels.
. Materials and methods . . Patient-specific lumped parameter model of HFpEF We created a patient-specific lumped parameter model of the cardiovascular system tuned to the hemodynamics of a patient with HFpEF. The study was approved by the local IRB board. Aortic blood flows measurements were obtained by 2D phase contrast cardiac magnetic resonance imaging and blood pressures as upper arm cuff pressures within 1 h of acquisition (male, 87 years, HR = 60 bpm). The patient was selected based on high diastolic filling pressures assessed by echocardiography (Lat E/e' = 14.5, Med E/e' = 18.4 at/above 15), a criteria defined by the American Society of Echocardiography .
/fcvm. . Initial conditions of pressure, flow, and volume for the first cardiac cycle were estimated using a minimum squared error function. The function starts solving for the initial conditions based on values found in literature. A time-periodic state is reached after five cardiac cycles.

. . . The vessels
The lumped parameter model of the main arterial vessels (aortic arch) was created with the following steps. We first constructed an anatomic model using SimVascular (http:// simvascular.github.io) (Point 1 in Figure 2) based on a computational tomography (CT) scan of the HFpEF patient. First, we drew pathlines along the centerlines of the arterial vessels of interest. We then segmented the CT scan by outlining each vessel lumen manually to produce 10 to 20 evenly spaced cross sections normal to the vessel centerlines (29). We then ran an automated script in SimVascular to calculate the resistances of each aortic branch vessel based on the average of its varying cross sectional area assuming Poiseuille flow (30). Each vessel element was reduced to three resistor subelements placed in series. At the end of each branch a resistor capacitor resistor (RCR) element tuned to achieve the patient's RCRRSub, resistance capacitance resistance element at the right subclavian artery; RCRRCar, resistance capacitance resistance element at the right carotid artery; RCRLCar, resistance capacitance resistance element at the left carotid artery; RCRLSub, resistance capacitance resistance element at the left subclavian artery; RCRCoro, resistance capacitance resistance element at the coronary arteries; RCRSyst, resistance capacitance resistance element of the systemic circulation.

. . . The valves
The mitral valve (MV) and the aortic valve (AV) (Point 4 in Figure 2) were modeled as pressure-actuated valves (33). The modeled heart valves only allow the flow to go through when the pressure on the inlet side is higher than the pressure on the outlet side. The state s(t) of the valve changes between 0 and 1 with 0 being a fully closed valve.

. . Lumped parameter model of the soft robot
The prototype of the soft robot is based on the McKibben actuator principle and consists of a single non-elastic pneumatically actuated bladder in a mesh which is wrapped around the ascending aorta ( Figure 3A). When injecting a user-determined pressure p act , the McKibben actuator will shorten circumferentially while thickening circumferentially and thus contract the aorta. The maximum width of the device prototype upon inflation is 12 mm. The inflation of the actuator induces a pressure on the ascending aorta, which in turn induces a change in volume inside the vessel. We used the capacitance-based approach proposed by Schampaert et al. to model the soft robot (34). The capacitor element receives an input by an experimentally defined pressure source p sr composed of a contraction part followed by a relaxation part (Equation 1) mimicking the soft robot's effect on a vessel (Section 2.2.2).

. . . Analytical definition
The parameters defining the soft robot actuation were the duration of the contraction t act , the delay with respect to the beginning of the cardiac cycle δt act (measured as R-wave of the electrocardiogram), the pressure inside the actuator p act , and the pressure inside the aorta p aorta . The duration of the actuator contraction was always 300 ms, such that for delays above 700 ms the actuator only deflated at the beginning of the next cycle. The duration of the contraction of the soft robot was t on while t off was the duration of a full actuation cycle including the relaxation. The times t on and t off varied for different p act and were determined experimentally (Supplementary material). Our experimental results confirmed that a higher p act leads to higher peak pressures while the peak pressure is reduced for higher p aorta .
HR − δt act is the time t normalized with respect to the cycle and actuation delay. g act (p aorta , p act ) is an experimental coefficient. Time coefficients k 1 and k 2 were determined experimentally and are equal to 13.04 and 34.66 respectively. δt act is the delay of the soft robot actuation. Coefficients c and d ensure continuity between the contraction and relaxations functions at moments t on and t off .
g act (p act , p aorta ) = a 1 * p 2 aorta + a 2 * p aorta + a 3 * p act + a 4 + a 5 * p aorta * p act (2) where the coefficients a 1 , a 2 , a 3 , a 4 and a 5 were obtained by minimizing the root mean square error (RMSE = (g act − g exp ) 2 ) with respect to experimental data (g exp ) using the Curve Fitting Tool from MATLAB (Version R2020b, MathWorks, Natick, MA, USA).

. . . Experimental characterization
To determine the coefficients of p sr , a closed-, non-actuated loop filled with a glycerine-water (40-60%) mix mimicking the viscosity of blood ( Figure 3A) was subjected to the contraction of the soft robot at 60 bpm with p act at 6, 7, 8, 9, 10, 11, and 12 psi. The duration of actuation t act was set to 300 ms which was previously found to be most effective in generating pressure between different parts of the cardiac cycle (35). We used 2.5 cm diameter silicone tubing (high-temperature silicone rubber, durometer 35A) with 1 mm wall thickness and total length of 2 m. A one-directional flow valve was used to inhibit back flow in the loop. The pressure p aorta was tuned to 50, 60, 70, 80, 90, 100, 125, and 150 mmHg by changing the glycerine-water mixture volume in the loop via a syringe (Figure 3). No significant kinking of the silicone vessel wall was observed during actuation of the soft robot at p aorta = 75mmHg as measured my magnetic resonance imaging (Supplementary material). Pressure and flow were measured using a pressure transducer (Mikro-Tip Catheter Transducer SPR-350s, Millar Inc., Houston, Texas, USA) and flow probe (ME14PXL294, Transonic Systems Inc., Ithaca, New York, USA). The measured data was recorded using the LabChart software (ADInstruments, Dunedin, New Zealand) and the PowerLab 8/35 combined with the Power Bridge Amp hardware (ADInstruments, Dunedin, New Zealand).

. . Coupling the patient-specific and soft robot lumped parameter models
The formulation of p sr was validated through RMSE analysis (Supplementary material). After validating the patientspecific lumped parameter model through comparison with hemodynamic data, the soft robot lumped parameter element was added between the first two ascending aorta vessel elements of the patient-specific cardiovascular model. A difference in capacitance and resistance between the experimental set-up and the patient leading to a different peak p sr . We solved this by multiplying p sr by a correcting factor. The factor was determined for p aorta = 100 mmHg. The resulting coupled system was solved numerically with a time step of 0.1 ms over 10 cardiac cycles. . Results
The waveform as well as the peak flow rate of the model, presented in Figure 4A, match well with the flow measured in the HFpEF patient. The RMSE between the in vivo flow data and the patient-specific model calculated over a full cardiac cycle is 36.1 mL/s. The RMSE can be partly explained through the small time shift between both flow peaks. The modeled PV loop of the LV ( Figure 4B) is composed of an isovolumetric pressure increase and decrease phases during LV contraction and relaxation as well as a small pressure increase during diastole. The flow splits into the main branches of the cardiovascular system are in line with the in vivo flow splits presented in literature (36). Moreover, the PV loop is closed, proving that the computational model has stabilized and is running identical cycles as it has reached a time periodic state.

. . The impact of the soft robotic device prototype on the patient
We present the influence of the actuator delay δt act between 50 and 950 ms with a step of 50 ms and the different p sr between 6 and 12 psi with a step of 2 psi on the patient where peak systolic pressure reduction, lower stroke work, and higher  An actuation delay of 650 ms results in maximal peak systolic pressure reductions (141 vs. 147 mmHg), allows for the highest cardiac output (5.34 vs. 5.17 L/min) and the lowest stroke work (10.6 kJ compared to 11.8 kJ) compared to baseline ( Figure 5). In this optimal case, the arterial elastance EA (EA = ESP/CO, where ESP is the end systolic pressure) is also reduced compared to baseline. The delay of 650 ms together with 300 ms actuation time mean that the soft robot relaxes shortly after the start of the ventricular contraction. However, the left atrial pressure increases upon actuation (11.5 mmHg compared to 11.2 mmHg). On the contrary, the worst case actuation was found to be at a delay of 0 ms and thus occurring simultaneous to the ventricular contraction creating higher peak systolic pressure (164 vs. 147 mmHg) and stroke work (11.8 vs. 11.1 kJ).
When comparing the PV loops of the optimal and baseline settings we see that the diastolic filling volume is increased while the end systolic volume is decreased resulting in a larger stroke volume. Moreover, the peak systolic pressure is decreased ( Figure 6). A second pressure and flow peak during diastole can be seen in all branches of the aorta in the optimal case (Figure 7). The pressure peak decays slowly over 0.3 s while the flow peak reaches its maximum after 0.1 s before disappearing after 0.2 s. The highest device-induced diastolic pressure increase at p act = 12 psi is seen in the coronary arteries (25 mmHg), followed by the carotid arteries (23 mmHg) and the descending aorta (11 mmHg). While, the overall cardiac output increased, this diastolic support did only lead to a small increase in coronary flow. In the worst case, the pressure and flow created by the soft robot adds up on top of the systolic activity of the heart, further increasing the afterload.
The flow splits toward the separate parts of the upper body are not known for the patient but their modeled values are presented in Table 2. The flow splits to the peripheral vessels remain similar to baseline in the optimal case. The soft robot increased absolute flows into each branch compared to baseline. The largest increase in flow is directed toward the descending aorta increasing the percentage of received flow. The absolute flow to the coronaries increased by six percent when the patient is supported by the extra-aortic support device.

. . The e ect of actuator pressure
We note a higher diastolic pressure increase during diastole with a higher p sr (Figure 8). In turn, this also leads to stronger reduction of peak systolic pressures. Similarly, the cardiac output increases with the higher p sr (5.22 L/min at 6 psi, 5.26 L/min at 8 psi, 5.29 L/min at 10 psi, and 5.34 L/min at 12 psi).

. Discussion
We created a patient-specific lumped parameter model of a HFpEF patient and a soft robotic extra-aortic counterpulsation device prototype to study the physiological implications of varying operating conditions. Full tuning to the patient-specific flow dynamics was achieved. We found that if operated in the optimal setting, we could achieve lower systolic pressure, a lower stroke work, and slightly improved perfusion of the coronaries at the cost of higher left atrial pressures.
This study evaluated a treatment strategy for HFpEF in a patient-specific model tuned to an individual subject. Previous work on modeling device-based treatment strategies in HFpEF .
/fcvm. . have focused on the use of averaged population characteristics and grouping based on hemodynamic phenotypes (11)(12)(13)(14)37). Burkhoff et al. extensively modeled four distinct phenotypes of HFpEF. The paper acknowledges that the group of patients with HFpEF induced by hypertension is highly variable in its hemodynamic characterization (11). We extended the prior population based approaches using an open-source computational framework to re-create a patient-specific hemodynamic model using the open-source framework Simvascular.
The device modeling was inspired by previous work on modeling counterpulsation devices in lumped parameter and 1D models (38,39). Our approach informed the device model directly from an experimental analysis of the studied .
/fcvm. . counterpulsation device prototype. Combining the model of the device prototype and patient-specific model allowed for patientspecific simulations of their interaction, which is critical given the number of different pathophysiologies found in HFpEF. Our simulation allows one to tune the settings of any potential device to the needs of a specific patient. In the future, the computational model could also allow to predict treatment effects prior to implantation.

FIGURE
Left ventricular PV loop for the best case ( ms delay), the worst case ( ms delay) and baseline (p act = 0 psi).
Using our simulation, we evaluated the effect of counterpulsation on the hemodynamics in HFpEF at different timings. When actuated during optimal setting actuation (p act = 12 psi, delay = 650 ms), an additional diastolic pressure and flow waveform is observed in all cardiovascular branches, which is inline with previous studies (34). The effect on systolic aortic pressure is similar to previous in vivo studies about IABPs and the C-Pulse System, however the cardiac output is increased to a much lesser extent (18,21,39,40). While the effect on decreasing aortic pressure is similar, the soft robotic device prototype covers a significantly smaller aortic arch length then   the C-pulse system. The small width of the actuator could allow for serial integration of multiple actuators in the future. Finally, we observed increased diastolic filling volume allowing for a larger stroke volume, however this comes at the cost of a higher left atrial pressure. In our simulation, the device is able to reduce the load on the heart, evidenced by the reduced stroke work.
The worst case actuation points toward the risks of the counterpulsation device with increased peak systolic pressure. A systolic actuation of the device can lead to an increase of peak systolic pressure by 17 mmHg, while the stroke work of the native heart increases by 0.7 kJ compared to baseline. The soft robotic actuator relies on accurate real-time sensing of the ventricular contraction to avoid harming the cardiovascular system, as has been shown in other pulsating soft robotic devices. While, the existing approaches to trigger IABPs might be useful for triggering of a soft robot, our results point toward the risks of pulsating technologies. In contrast to IABPs, the extraaortic approach mitigates risks associated with infection and thromboembolic considerations, while several other risk factors will need to be addressed such as long-term functioning of the device as well as implications for the aortic valve.
We were also able to show the importance of the chosen p act in tuning the effect of the soft robotic device prototype on the patient's physiology. A higher p act led to lower systolic peak pressures as well as higher cardiac output. The ability to tune the impact allows for scenario-specific adaptation such as rest and exercise but also for gradual increase in support to the heart. Furthermore, a bridge to recovery application can be envisioned through support modulation. Based on the similarity of the flow waveforms between the soft robot and IABPs, no impact on cerebral function is expected (41).
Clinical considerations: The usefulness of counterpulsation strategies is highly debated (35). While IABP are used to bridge patients with HFrEF effectively, the underlying physiological mechanism remain unclear. The current recommendation made by the American College of Cardiology, the American Heart Association, and the Heart Failure Society of America is to lower the systolic blood pressure (SBP) under 130 mmHg for HFpEF patients with hypertension (42). However, Zouein et al. have shown antihypertensive medication to be ineffective in diminishing mortality in patients with HFpEF (4). If lowering systolic pressure (afterload) is beneficial in HFpEF, our study suggest that extra-aortic counterpulsation could be beneficial. However, the design of the current device prototype requires further optimization in order to displace sufficient volume.
Limitations: The study constitutes a framework for patientspecific hemodynamic modeling of device based treatment. The cohort of patients will need to be increased to cover a larger variety of patients as well as the heterogeneity in pathophysiology found in HFpEF. Moreover, the simulation currently lacks integration of the autoregulatory control mechanisms such as the baroreflex as well as coronary blood .
/fcvm. . flow regulations. Furthermore, the effect of the silicon vessel on the evaluation of the soft robots was not specifically addressed in this work. Ongoing experimental work and in-vivo trials with the soft robotic actuator will enable further validation of the findings of the lumped parameter model presented.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary material, further inquiries can be directed to the corresponding author/s.

Ethics statement
The studies involving human participants were reviewed and approved by IRB Veteran Association Palo Alto. The patients/participants provided their written informed consent to participate in this study.

Author contributions
MA, DE, SD, and IC contributed to conception and design of the study. IC organized the data acquisition and provided clinical context. MA and SD performed the experiments and wrote the first draft of the manuscript. JP and AM provided access to the original code. All authors contributed to manuscript revision, read, and approved the submitted version.
Funding SD received funding from the Swiss National Research Foundation P2EZP2_188964. SD, IC, and DE received funding from the NIH grant RO1 HL131823. AM and JP received funding from NIH grants R01EB029362 and R01EB018302.