Finite state machine implementation for left ventricle modeling and control

Simulation of a left ventricle has become a critical facet of evaluating therapies and operations that interact with cardiac performance. The ability to simulate a wide range of possible conditions, changes in cardiac performance, and production of nuisances at transition points enables evaluation of precision medicine concepts that are designed to function through this spectrum. Ventricle models have historically been based on biomechanical analysis, with model architectures constituted of continuous states and not conducive to deterministic processing. Producing a finite-state machine governance of a left ventricle model would enable a broad range of applications: physiological controller development, experimental left ventricle control, and high throughput simulations of left ventricle function. A method for simulating left ventricular pressure-volume control utilizing a preload, afterload, and contractility sensitive computational model is shown. This approach uses a logic-based conditional finite state machine based on the four pressure-volume phases that describe left ventricular function. This was executed with a physical system hydraulic model using MathWorks’ Simulink® and Stateflow tools. The approach developed is capable of simulating changes in preload, afterload, and contractility in time based on a patient’s preload analysis. Six pressure–volume loop simulations are presented to include a base-line, preload change only, afterload change only, contractility change only, a clinical control, and heart failure with normal ejection fraction. All simulations produced an error of less than 1 mmHg and 1 mL of the absolute difference between the desired and simulated pressure and volume set points. The acceptable performance of the fixed-timestep architecture in the finite state machine allows for deployment to deterministic systems, such as experimental systems for validation. The proposed approach allows for personalized data, revealed through an individualized clinical pressure–volume analysis, to be simulated in silico. The computational model architecture enables this control structure to be executed on deterministic systems that govern experimental left ventricles. This provides a mock circulatory system with the ability to investigate the pathophysiology for a specific individual by replicating the exact pressure–volume relationship defined by their left ventricular functionality; as well as perform predictive analysis regarding changes in preload, afterload, and contractility in time.


Introduction
Every year since 1919, cardiovascular disease (CVD) accounted for more deaths than any other major cause of death in the United States [1]. Based on data collected by the National Health and Nutrition Examination Survey (NHANES), CVD was listed as the underlying cause of death in 30.8% of all deaths in 2014, accounting for approximately 1 of every 3 deaths in the U.S., while CVD attributed to 53.8% of all deaths in that year. Additionally, data accumulated from 2011 to 2014 revealed that approximately 92.1 million American adults currently have one or more types of CVD and that by 2030, projections estimate that 43.9% of the U.S. population will have some form of this disease.
Research has revealed that CVD is a leading contributor to Congestive Heart Failure (CHF) [2]. CHF is a medical condition that occurs when the heart is incapable of meeting the demands necessary for maintaining an adequate amount of blood flow to the body, resulting in ankle swelling, breathlessness, fatigue, and potentially death [2]. In 2012, the total cost for CHF alone was estimated to be $30.7 billion with 68% attributed to direct medical costs. Furthermore, predictions indicate that by 2030, the total cost of CHF will increase almost 127% to an estimated $69.7 billion [1]. This prediction is based on data that revealed that one-third of the U.S. adult population has the predisposing conditions for CHF. With research revealing that 50% of people who develop CHF will die within 5 years of being diagnosed [1,3], the need to evaluate treatments for this widening patient population is of growing importance.
One treatment alternative for patients with late-stage CHF is the use of a ventricular assist device (VAD) to directly assist with the blood flow demands of the circulatory system [2]. Implantable VADs have proven their potential as a quickly implemented solution for bridge to recovery, bridge to transplant, and destination therapy [4]. Given the severity of CHF, and the impending need for supplemental support from these cardiac assist devices, effective methods of identifying the recipient cardiovascular profile and matching that to the operation of the VAD is critical to the success of the intervention.
The effectiveness of CHF diagnosis and treatment therapy depends on an accurate and early assessment of the underlying pathophysiology attributed to a specific type of CVD, typically by means of analyzing ventricular functionality [2,5,6]. Clinical application of non-invasive cardiac imaging in the management of CHF patients with systolic and/or diastolic dysfunction has become the standard with the use of procedures such as echocardiography [7][8][9][10]. Echocardiography is a non-invasive ultrasound procedure used to assess the heart's structures and functionality, to include the left ventricular ejection fraction (LV EF ), left ventricular end-diastolic volume (LV EDV ), and left ventricular end-systolic volume (LV ESV ). Three-dimensional echocardiography of adequate quality has been shown to improve the quantification of left ventricular (LV) volumes and LV EF , as well as provide data with better accuracy when compared with values obtained from cardiac magnetic resonance imaging [2,11]. At present, echocardiography has been shown to be the most accessible technology capable of diagnosing diastolic dysfunction; therefore, a comprehensive echocardiography examination incorporating all relevant two-dimensional and Doppler data is recommended [2]. Doppler techniques allow for the calculation of hemodynamic variations, such as stroke volume (SV) and cardiac output (CO), based on the velocity time integral through the LV outflow tract area.
A left ventricular pressure-volume (LV-PV) analysis, employing hemodynamic principles, has effectively performed as a basis for understanding cardiac physiology and pathophysiology for decades [12,13]. A LV-PV analysis has been primarily restricted to clinical investigations in a research environment; therefore, it has not been extensively used due to the invasive nature of the procedure [14,15]. A broader predictive application for detecting and simulating CHF is more easily attainable with the development of single-beat methodologies that only rely on data collected through non-invasive techniques. These techniques include echocardiographic measurements of the left ventricular volume (LVV), Doppler data, the peripheral estimates of left ventricular pressure (LVP), and the timing of the cardiac cycle [16][17][18][19][20][21].
Utilizing data obtained non-invasively, population and patient-specific investigations can be conducted by simulating the LV-PV relationship obtained through the PV analysis by means of a mock circulatory system (MCS) [22,23]. An MCS is a mechanical representation of the human circulatory system, essential for in vitro evaluation of VADs, as well as other cardiac assist technologies [24][25][26][27][28][29]. An MCS effectively simulates the circulatory system by replicating specific cardiovascular conditions, primarily pressure [mmHg] and flow rate [mL/s], in an integrated bench-top hydraulic circuit. Utilizing these hydraulic cardiovascular simulators and data obtained through a clinical PV analysis, the controls that govern the LV portion of the MCS could be driven to produce the PV relationship of: a CVD profile, specific population, or patient [30]. With research revealing the increasing need for these medical devices [31], a comprehensive in vitro analysis could be completed to assure a particular cardiac assist device treatment will be effective beforehand. The ability of an MCS to be able to replicate the exact PV relationship that defines the pathophysiology for a specific individual allows for a robust in vitro analysis to be completed, and a "patient specific diagnosis" created, ensuring a higher standard of patient care [32,33].
The following is how this manuscript is presented. "Background" section summaries the principal theories governing the modeling of the PV relationship, its background in simulating cardiovascular hemodynamics within an MCS, and how a PV loop controller should perform for subsequent in vitro testing. "Method" section presents the proposed methodology for developing LV-PV control functionality is presented and utilizes a logic-based conditional finite state machine (FSM) and a physical system modeling approach, then the experimental results are presented in "Results" section. "Discussion" section concludes with a discussion regarding the results of this investigation, followed by "Conclusion" section which outlines the limitations of the approach and future investigations.

Pressure-volume relationship
The efficacy of the PV relationship, often referred to as a PV loop, to describe and quantify the fundamental mechanical properties of the LV was first demonstrated in 1895 by Otto Frank [34]. Frank represented the cardiac cycle of ventricular contraction as a loop on a plane defined by ventricular pressure on the vertical axis and ventricular volume on the horizontal. By late twentieth century, the PV analysis was considered the gold standard for assessing ventricular properties, primarily due to the researched conducted by Suga and Sagawa [35][36][37]. Yet, this approach has failed to become the clinical standard for evaluating LV functionality due to the invasive nature of the procedure [14,15]. However, due to recent advances single-beat methodologies, the practical application for PV analysis is expanding [18][19][20]. Most recently are the efforts published in 2018 by Davidson et al. with regard to the development of a beat-by-beat method for estimating the left ventricular PV relationship using inputs that are clinically accessible in an intensive care unit (ICU) setting and are supported by a brief echocardiograph evaluation [20].
There has been extensive clinical and computational research into understanding the PV relationship, which is presented in Fig. 1 [12,21,30,38]. However, for the purpose of repeatability within a MCS, the culmination of this knowledge can be summarized by simplifying the performance of the LV through three principal factors: preload, afterload, and contractility [24,25]. These have significant implications on VAD performance [39].
A schematic of the LV pressure-volume loop in a normal heart is presented in Fig. 1a. In Phase I, ventricular filling occurs with only a small increase in pressure and a large increase in volume, guided along the EDPVR curve. Phase I can additionally be divided in two sub-phases, rapid filling governed by elastance of the ventricle and atrial systole that brings the ventricle into optimal preload for contraction. Phase II constitutes the first segment of systole called isovolumetric contraction. Phase III begins with the opening of the aortic valve; ejection initiates and LV volume falls as LV pressure continues to increase. Phase III can be divided into two sub-phases: rapid ejection and reduced ejection. Isovolumetric relaxation begins after the closure of the aortic valve constituting Phase IV.
Ventricular preload refers to the amount of passive tension or stretch exerted on the ventricular walls (i.e. intraventricular pressure) just prior to the systolic contraction [14,29]. This load determines the end-diastolic sarcomere length and thus the force of contraction. Because the true sarcomere length is not easily measured clinically, preload is typically measured by ventricular pressure and volume at the point immediately preceding isometric ventricular contraction. This correlation is described through the end-systolic pressure-volume relationship (ESPVR); as well as through the end-diastolic pressure-volume relationship (EDPVR). The effects of increasing preload on the PV relationship is displayed in Fig. 1b; reduced isovolumetric contraction period and increased stroke volume.
Afterload is defined as the forces opposing ventricular ejection [14]. Effective arterial elastance (E a ) is a lumped measure of total arterial load that incorporates the mean resistance with the pulsatile factors that vary directly with heart rate, systemic vascular resistance, and relates inversely with total arterial compliance. E a is directly defined as the ratio of left ventricular end-systolic pressure (LV ESP ) to SV. In practice, another measure of afterload is the LV ESP at the moment ventricular pressure begins to decrease to less than systemic arterial pressure. The effects of increasing afterload are presented in Fig. 1c; increase in peak systolic pressure and decrease in stroke volume.
A acceptable clinical index of contractility that is independent of preload and afterload has not been completely defined [29]. In non-pathological conditions, contractility is best described by the pressure-volume point when the aortic valve closes. Contractility is typically measured by the slope of the ESPVR line, known as E es , which is calculated as P V [38]. An additional index of contractility is dP/dt max which is the derivative of the maximum rate of ventricular pressure rise during the isovolumetric period. The effects of increasing contractility on the PV relationship is revealed in Fig. 1d; revealing the ability for the stroke volume to accommodate with increasing peak systolic pressure.
For a given ventricular state, there is not just a single Frank-Starling curve, rather there is a set or family of curves [29]. Each curve is determined by the driving conditions of preload, afterload, and inotropic state (contractility) of the heart. While deviations in venous return can cause a ventricle to move along a single Frank-Starling curve, changes  [30]). a Schematic of LV pressure-volume loop in a normal heart. In Phase I, preceding the opening of the mitral valve, ventricular filling occurs with only a small increase in pressure and a large increase in volume, guided along the EDPVR curve. Phase II constitutes the first segment of systole called isovolumetric contraction. Phase III begins with the opening of the aortic valve; ejection initiates and LV volume falls as LV pressure continues to increase. Isovolumetric relaxation begins after the closure of the aortic valve constituting Phase IV. b Effects of increasing preload on a LV-PV loop with afterload and contractility remaining constant. Loop 2 has an increased preload compared to loop 1 by rolling the arterial elastance (E a ) line parallel while keeping the slope (E a ) constant, resulting in an increase in SV. c Effects of increasing afterload on a LV-PV loop with preload and contractility held constant. This consists of increasing the slope of the E a line. d Effects of increasing contractility on a LV-PV loop with preload and afterload remaining constant. This consists of increasing the slope (E es ) of the ESPVR line. Note that in b, c, and d, loop 2 represents the increase in the respective principle factor, i.e. preload, afterload, and contractility, when compared to loop 1 in the driving conditions can cause the PV relationship of the heart to shift to a different Frank-Starling curve. This allows clinicians to diagnose the pathophysiological state of a dysfunctional heart by analyzing the PV relationship of a patient. Additionally, it provides the ability to simulate diseased states: heart failure [14], valvular disease [29], or specific cardiovascular dysfunction seen in pediatric heart failure [40].

Pressure-volume loop computational modeling
Comprehensive computationally modeling of the LV-PV relationship has been effectively reported since the mid-1980s, following the extensive work completed by Suga and Sagawa [34][35][36]. In 1986, Burkhoff and Sagawa first developed a comprehensive analytical model for predicting ventricular efficiency utilizing Windkessel modeling techniques and an understanding of the PV relationship principles previously developed by Suga and Sagawa. With the advancement and routine use of innovative technologies in the early twenty-first century (e.g. conductance catheter, echocardiography), there was a significant increase in research efforts to determine the potential clinical applications [12][13][14][15], improving predictive strategies [16][17][18][19], and refining computational models [41][42][43].
An elastance-based control of an electrical circuit analogue of a closed circulatory system with VAD assistance was developed in 2009 by Yu et al. [42]. Their state-feedback controller was designed to drive a voice coil actuator to track a reference volume, and consequently generate the desired ventricular pressure by means of position and velocity feedbacks. The controller was tested in silico by modifying the load conditions as well as contractility to produce an accurate preload response of the system. The MCS analogue and controller architecture was able to reproduce human circulatory functionality ranging from healthy to unhealthy conditions. Additionally, the MCS control system developed was able to simulate the cardiac functionality during VAD support.
In 2007, Colacino et al. developed a pneumatically-driven mock left ventricle as well as a native left ventricle model and connected each model to a numerical analogue of a closed circulatory system comprised of systemic circulation, a left atrium, and inlet/ outlet ventricular valves [43]. The purpose of their research was to investigate the difference between preload and afterload sensitivity of a pneumatic ventricle, when used as a fluid actuator in a MCS, when compared to elastance-based ventricle computational model. Their research concluded that the elastance-based model performed more realistically when reproducing specific cardiovascular scenarios and that many MCS designs could be considered inadequate, if careful consideration is not made to the pumping action of the ventricle. Subsequent in vitro testing utilizing this control approach successfully reproduced an elastance mechanism of a natural ventricle by mimicking preload and afterload sensitivity [25]. Preload was modified by means of manually changing the fluid content of the closed loop hydraulic circuit, while afterload was varied by increasing or decreasing the systemic arterial resistance within a modified Windkessel model.

Recent advancements in contractility-based control
An MCS simulates the circulatory system by accurately and precisely replicating specific cardiovascular hemodynamic variables, mainly the respective pressure (mmHg) and flow rate (mL/s) for key circulatory constituents, in an integrated bench-top hydraulic circuit [23]. While this human circulatory system model is not an all-inclusive replacement for an in vivo analysis of a cardiac assist device's design, it is an effective method of evaluating fundamental design decisions beforehand by determining its influence on a patient's circulatory hemodynamics in a safe and controlled environment. Published research efforts typically either involve the development of the system [22,25,26,[44][45][46] or the dissemination of the results of a particular in vitro investigation [27,28].
In 2017, Wang et al. was able to replicate the PV relationship with controllable ESPVR and EDPRV curves on a personalized MCS based on an elastance function for use in the evaluation of VADs [21]. The numerical elastance models were scaled to change the slopes of the ESPVR and EDPVR curves to simulate systolic and diastolic dysfunction. The results of their investigation produced experimental PV loops that are consistent with the respective theoretical loop; however, their model only includes a means of controlling preload and contractility with no afterload control. Their model assumes afterload remains constant regardless of preload changes; due to the Frank-Starling mechanism, the ventricle reached the same LV ESV despite an increase in LV EDV and preload.
Jansen-Park et al., 2015, determined the interactive effects between a simulated patient with VAD assistance on an auto-regulated MCS which includes a means of producing the Frank-Starling response and baroreflex [24]. In their study, a preload sensitive MCS was developed to investigate the interaction between the left ventricle and a VAD. Their design was able to simulate the physiological PV relationship for different conditions of preload, afterload, ventricular contractility, and heart rate. The Frank-Starling mechanism (preload sensitivity) was modeled by regulating the stroke volume based on the measured mean diastolic left atrial pressure, afterload was controlled by modifying systemic vascular resistance by means of an electrically controlled proportional valve, and contractility was changed depending on the end diastolic volume. The effects of contractility, afterload, and heart rate on stroke volume were implemented by means of two interpolating three-dimensional look-up tables based on experimental data for each state of the system. The structure of their MCS was based on the design developed by Timms et al. [27]. The results of their investigation revealed a high correlation to published clinical literature.
In 2011, Gregory et al. was able to replicate a non-linear Frank-Starling response in a MCS by modifying preload by means of opening a hydraulic valve attached to the systemic venous chamber [44]. Their research was able to successfully alter left and right ventricular contractility by changing preload to simulate the conditions of mild and severe biventricular heart failure. The EDV offset and a sensitivity gain were manually adjusted through trial and error to produce an appropriate degree of contractility with a fixed ventricular preload. The shape of the ESPVR curve was then modified by decreasing MCS volume until the ventricular volumes approached zero. These efforts, validated using published literature, improved a previously established MCS design developed by Timms et al. [28].
These control architectures were primarily hardware determined, rather than software-driven. In some cases, reproducibility is inhibited due to the tuning of hemodynamic conditions by manually adjusting parameters until a desired response is achieved. Utilizing a conditional logic-based conditional finite state machine (FSM) and physical system modeling control approach, a software-driven controller could be developed to respond to explicitly-defined preload, afterload, and contractility events. This would enable the regulation of the PV relationship within an MCS's LV section, without the limitation of dedicated hardware.

Logic-based finite state machine (FSM) and physical system modeling tools
MathWorks' Simulink ® is a model based design tool utilized for multi-domain physical system simulation and model-based design [47]. Simulink ® provides a graphical user interface, an assortment of solver options, and an extensive block library for accurately modeling dynamic system performance. Stateflow ® is a toolbox found within Simulink ® for constructing combinatorial and sequential decision-based control logic represented in state machine and flow chart structure. Stateflow ® offers the ability to create graphical and tabular representations, such as state transition diagrams and truth tables, which can be used to model how a system reacts to time-based conditions and events, as well as an external signal. The Simscape ™ toolbox, utilized within the Simulink ® environment, provides the ability to create models of physical systems that integrate block diagrams acknowledged by real-world physical connections. Dynamic models of complex systems, such as those with hydraulic and pneumatic actuation, can be generated and controlled by assembling fundamental components into a schematic-based modeling diagram. An additional toolbox that was utilized in this approach was the Simscape Fluids ™ toolbox which provides component libraries for modeling and simulating fluid systems. The block library for this toolbox includes all the necessary modules to create systems with a variety of domain elements, such as hydraulic pumps, fluid reservoirs, valves, and pipes. The advantage of using these toolbox libraries is that the blocks are version controlled and conformal to regulatory processes that mandate tractable computational modeling tools.

Overview of methodology and model architecture
A method for simulating LV-PV control functionality utilizing explicitly defined preload, afterload, and contractility is needed for cardiovascular intervention assessment. The resulting solution must be capable of being compiled for hardware control of an MCS; deterministic processing compatible logic and architecture that would enable runtime setpoint changes. The approach used was a logic-based conditional finite state machine (FSM) based on the four PV phases that describe left ventricular functionality developed with a physical system hydraulic plant model using Simulink ® . The proposed aggregate model consists of three subsystems to include: a preload/afterload/contractility-based setpoint calculator ("PV loop critical point determination" section), a FSM controller ("PV loop modeling utilizing a state machine control architecture approach" section), and a hydraulic testing system ("Hydraulic testing model utilizing MathWorks' Simulink ® and SimscapeTM toolbox" section). The last subsystem acts as the simulated plant to evaluate the control architecture that is formed by the first two subsystems. The proposed method allows for multiple uses which include the simulation of parameter effects in time and the simulation of personalized data, revealed through an individualized clinical PV analysis. This method provides the means to be simulated in silico and can be subsequently compiled for control of in vitro investigations. This provides an MCS with the ability to investigate the pathophysiology for a specific individual by replicating the exact PV relationship defined by their left ventricular functionality; as well as perform predictive analysis regarding changes in preload, afterload, and contractility with time. Of critical importance were the non-isovolumetric state behavior: non-linear EDPVR curve, rate-limited ejection, and energy-driven model of contraction. This investigation was developed utilizing Matlab R2017b and a Dell T7500 Precision workstation with 8.0 gigabytes of RAM, a Dual Core Xeon E5606 processor, and a Windows 7 64-bit operating system.

PV loop critical point determination
A preload, afterload, and contractility sensitive computational model was developed utilizing Simulink ® for determining critical points for switching between PV loop states; the four phases described in . These can be resolved by the three equations that describe ESPVR, EDPVR, and E a . ESPVR is typically described as a linear equation with a positive slope (E es ) and a negative or positive y-intercept, EDPVR can be defined with a third-order polynomial, while E a is also linear and has a negative slope with a positive y-intercept [13]. Eqs. 1, 2, and 3 define the system of equations used to produce the critical points, where ESPVR, EDPVR, and E a are Eqs. 1, 2, and 3 respectively.
The point where Eqs. 1 and 3 intercept is LV ESV and LV ESP and solving produces Eqs. 4 and 5.
(1) Lastly, due to isovolumetric contraction, LV EICV equals LV EDV . The final unknown variable value to complete the four-phase cycle is LV EICP . This is resolved by utilizing an offset value based on LV ESP . , an open-source program used to extract data from images, these coefficients can be obtained from a plot of a patient's left ventricular pressure-volume analysis of preload change.

PV loop modeling utilizing a state machine control architecture approach
Utilizing Simulink ™ Stateflow ® , a sequential decision-based control logic represented in Mealy machine structure form was developed to control the transition between LV-PV phases. A Mealy machine is appropriate because this application requires that the output values are determined by both its current state and the current input values. A state transition diagram is presented in Fig. 3. The Variables in the block are parameters that are held constant: Piston cross-sectional area (A), b3, b2, b1, b0, Isovolumetric Rate, Isovolumetric Contraction Offset, Systolic Ejection Rate, and Systolic Ejection Offset. The Inputs are parameters that can change with time and are (7) Isovolumetric Rate is defined as the rate of change in the output variable, F, during isovolumetric relaxation and contraction. For isovolumetric relaxation, this rate is one-third the magnitude when compared to isovolumetric contraction. The Isovolumetric Contraction Offset is defined as the value subtracted from the LV EDV to start  the initialization of the Phase 2 state to compensate for the radius of curvature created due to transitioning from fill to eject, as well as the means by which end-diastolic pressure and volume are clinically quantified. The Systolic Ejection Rate is defined as  Upstream to the mitral valve is a Hydraulic Reference source block governed by the EDPVR curve function with respect to simulated volume, LVV, and increased by an offset of 2 mmHg to ensure proper flow through the mitral check valve. This establishes a dynamic LAP, the initial pressure condition of the left heart. LAP is outputted from the model here for plotting purposes. Downstream to the aortic valve is a Spring-Loaded Accumulator block. This block element consists of a preloaded spring and a fluid chamber. As the fluid pressure at the inlet of the accumulator becomes greater than the prescribed preload pressure, fluid enters the accumulator and compresses the spring, creating stored hydraulic energy. A decrease in the fluid pressure causes the spring to decompress and eject the stored fluid into the system. The spring motion is restricted by a hard stop when the fluid volume becomes zero, as well as when the fluid volume is at the prescribed capacity of the fluid chamber. These settings are utilized to regulate the compliance, V P , of the aorta. Immediately following is Hydraulic Pressure Sensor measuring AoP.
Additionally, a needle valve was positioned downstream to the aortic valve to simulate the resistance to flow contributed to the branching arteries of the aortic arch, as well as provide the capability to simulate the effects of increasing and decreasing resistance with time. As was previously stated, all block element parameter values that were modified are noted in the diagram presented in Fig. 4, while any parameters not noted were left standard to the block's original parameter values. For any element parameter denoted as 'Variable' , these values were not left constant for all simulations presented. For each simulation, these values are displayed in Table 1.

Results
The computational model effectively executed the trials assessing the performance of the FSM architecture. Solver settings and simulated fluid type were held constant through the analysis. The results presented were produced with the MathWorks' ode14x (fixedstep, extrapolation) using a fundamental sampling time of 1 1024 s. This solver was chosen to accelerate the simulations and ensure the resultant model is compatible with deterministic hardware systems. Validation of this solver was performed against a variablestep variable-order solver (ODE15 s) to ensure accuracy. The fluid selected is a glycerol/ water mixture with a fluid density of 1107.1 kg/m 3 and a kinematic viscosity of 3.3 centistoke [49]. These characteristics equate to a fluid temperature of 25 °C or 77 °F.
The input variables utilized for each presented simulation are displayed in Table 1, while the results of each simulation are displayed in Table 2. All simulations were performed utilizing discrete changes, evenly incremented between the designated initial and final LV ESP , LV ESV , LV EDP , and LV EDV over a 10 s total simulation time. Each discrete variable is controlled by means of a Lookup Table element block that outputs the modified variable value, dependent on the specific cycle count number. Note, any variable presented as a vector, changes with each cycle count, i.e. [1, 2, 3, · · · , n] where the nth value represents the input variable value for the entirety of the corresponding cycle. If a simulation has more cycles than input vector elements, then the system continues with a zero-order hold of the last value.
The parameters for the Spring-Loaded Accumulator block were developed based on a desired LVP response due to aortic compliance. The desired response consisted of a physiological correct AoP waveform and a peak-to-peak AoP amplitude of approximately 40 mmHg, corresponding to a normal range of 120/80. The baseline of this response was created at a heart rate of 60 bpm and a compliance of 1. This corresponded to an Isovolumetric Rate of 225 N*sample/s, a Resistance value of 0.03, a Fluid Chamber Capacity of 517.15 mmHg, a Preload Pressure of 0.01 psi, and a Pressure at Full Capacity of 10.01 psi. Given the relationship 1 R * C = I , where R is resistance, C is compliance, and I is the Impedance, I was held constant for all simulations using I = 33.333. For the simulations that required a heart rate beyond 60 bpm, the Isovolumetric Rate had to be consequently increased. Utilizing this relationship to sustain a peak-to-peak AoP amplitude of 40 mmHg, the Fluid Chamber Capacity and the Preload Pressure was held constant, while Resistance and Pressure at Full Capacity was modified to produce the desired heart rate while sustaining aortic performance. Lastly, the Initial Volume of Fluid for each simulation was calculated to create an initial LVP corresponding to LV ESP . This was done to decrease the amount of initial cycles necessary to achieve simulation stability to 1. All values utilized for these parameters are presented in Table 1. Error was calculated as the absolute value of the difference between the desired and simulated LV ESP , LV ESV , LV EDP , and LV EDV .
A LV-PV loop; LVP, LAP, and AoP versus time; and volume versus time graphs for the 10 s total simulation time was presented for each simulation. Note, the driving force

Computational model verification
The LV-PV loop critical point computational model and FSM approach was effective at driving the hydraulic testing model to produce the characteristic LV-PV relationship as presented in Fig. 5. The computational model parameters are the same as those presented in Fig. 2. As can be shown from the graph, with known ESPVR, EDPVR, and E a curves, the computational model successfully provided the correct LV ESP , LV ESV , LV EDP , LV EDV , LV EIRP , and LV EIRV transition points within the state transition logic to produce the prescribed LV-PV relationship. Table 1 contains all input parameters and Table 2 presents the results of all simulations performed. For each LV-PV loop graph, the initial LV end-systolic and end-diastolic datasets are denoted with circle points. Figure 5a displays the LV-PV loop based on data collected using DataThief on loop 1 of Fig. 1b. The results presented reveal an error between the desired and simulated end-systolic and end-diastolic transition points in the datasets of less than 1 mmHg and 1 mL, respectively.
The system takes one cycle to initialize from a rest state before the control topology functions consistently for the remainder of the simulation. Additionally, the isovolumetric and systolic offsets and rates, necessary to achieve this response are noted in Table 1. Figure 5a also presents the LVP, LAP, and AoP versus time and volume versus time graphs for the 10 s total simulation time. The reference force [N] produced by the FSM as well as the piston position [cm] can be derived from these time graphs. The outlined approach was effective at simulating the characteristic LV-PV relationship. Preload, afterload, and contractility changes in time were simulated by means of manipulating the input variables of the computational model via evenly-spaced discrete increments that change per cycle count. The LV-PV loop, pressure versus time, and volume versus time graphs are presented for each simulation. Displayed in a is the derived LV-PV loop, based on the computational model parameters determined using DataThief on loop 1 of Fig. 1b and presented in Fig. 2. The parameters for this LV-PV loop constitutes the initial conditions for the subsequent simulations. b presents the system correctly responding to a discrete change in preload. c reveals the correct afterload change response to the PV relationship. d displays the correct system response to contractility change. Each simulation was run for a total simulation time of 10 s and the system takes one cycle before it settles. The system functions consistently for every preceding cycle. The heart rate begins at approximately 60 bpm for each simulation. The reference force Preload, afterload, and contractility changes in time As presented in Fig. 5b-d, the outlined approach was effective at simulating preload, afterload, and contractility changes in time by means of discretely manipulating the computational model over time. The initial parameters of the computational model are the same as those presented in Fig. 5a and presented in Table 1. Presented for each simulation is the LV-PV loop; LVP, LAP, and AoP versus time; and volume versus time graphs for the 10 s total simulation time.
As shown in Fig. 5b, the system displays the correct preload change response to the PV relationship as displayed in Fig. 1b. The E a was initially defined by the equation P = −1.7504(V) + 185.02 . The y-axis intercept was increased from 185.02 mmHg at a rate of 5 mmHg per cycle, ending with a y-axis intercept of 215.02 mmHg for the last completed cycle. The results report an error of less than 1 mmHg and 1 mL for all targeted pressures and volumes.
Presented in Fig. 5c, the system reveals the correct afterload change response to the PV relationship as shown in Fig. 1c. E a is initially defined by the equation P = −1.7504(V) + 185.02 . The y-axis intercept was decreased from 185.02 mmHg at a rate of 15 mmHg per cycle, ending with a y-axis intercept of 110.02 mmHg for the last completed cycle. The slope of the E a was decreased from − 1.7504 mmHg/mL concluding with a slope of − 1.0408 mmHg/mL. This rate of change for the E a slope was derived from the 15 mmHg per cycle y-axis rate of increase to achieve a consistent x-intercept, as shown in Fig. 1c. The results indicate an error of less than 1 mmHg and 1 mL for all targeted datasets.
As presented in Fig. 5d, the system displays the correct contractility change response to the PV relationship as revealed in Fig. 1d. The ESPVR curve is initially defined by the equation P = 2.9745(V) − 17.133 . The slope of the ESPVR curve was decreased from 2.9745 mmHg/mL, concluding with a slope of 1.2245 mmHg/mL for the last completed cycle. The results report an error of less than 1 mmHg and 1 mL for all targeted pressures and volumes. Figure 6 displays the results of simulating Heart Failure with Normal Ejection Fraction (HFNEF) and the Control developed by means of a preload reduction analysis conducted in 2008 by Westermann et al. [50] and presented in Fig. 1 of their investigation. The ESPVR, E a , and EDPVR curve coefficients were developed utilizing DataThief to find the associated LVESP, LV ESV , LV EDP , and LV EDV for the initial and final loops, as well as evaluate the EDPVR curve. These datasets were analyzed over a 10 s total simulation time and for each simulation are the LV-PV loop; LVP, LAP, and AoP versus time; and volume versus time graphs. Both simulations reflect a mean heart rate [bpm] within the range of mean values noted in the reference material. All parameter values are presented in Table 1 and the results are in Table 2.

Clinical assessment of outlined approach
The Control is presented in Fig. 6a. The ESPVR curve was found to be defined by the equation P = 1.2407(V) + 33.857 and the EDPVR curve was found to be P = 2.6928E − 7(V ) 3 + −9.3013E − 6(V ) 2 + 0.026968(V ) + 2.9515 . E a is initially defined by the equation P = −1.1365(V) + 211.17 and defined by the equation P = −1.4501(V) + 160.11 for the final cycle. The slope and y-intercept of E a was divided into evenly-spaced increments to constitute 4 intermediate discrete steps between the initial and final cycle parameters. The results indicate an error of less than 1 mmHg and 1 mL for all targeted datasets. HFNEF is presented in Fig. 6b. The ESPVR curve was found to be P = 0.99741(V) + 72.586 and the EDPVR curve was found to be 010234 . E a is initially defined by the equation P = −1.4054(V) + 235.76 and defined by the equation P = −1.3754(V) + 160.43 for the final cycle. The slope and y-intercept of E a was divided into evenly-spaced increments to constitute 4 intermediate discrete steps between the initial and final cycle parameters. The results produced an error of less than 1 mmHg and 1 mL for all targeted datasets. Fig. 6 The outlined approach was effective at simulating Heart Failure with Normal Ejection Fraction (HFNEF) and the Control developed by means of a preload reduction analysis conducted in 2008 by Westermann et al. [50] and presented in Fig. 1 of their investigation. The ESPVR, E a , and EDPVR curve coefficients were developed utilizing DataThief to find the associated LV ESP , LV ESV , LV EDP , and LV EDV for the initial and final loops, as well as evaluate the EDPVR curve. These datasets were analyzed over a 10 s total simulation time and for each simulation is the LV-PV loop; LVP, LAP, and AoP versus time; and volume versus time graphs. a presents the Control where the slope and y-intercept of E a was divided into evenly-spaced increments to constitute 4 intermediate discrete steps between the initial and final cycle parameters. HFNEF is presented in b. The slope and y-intercept of E a was also divided into evenly-spaced increments to constitute 4 intermediate discrete steps between the initial and final cycle parameters. For both simulations, the results produced an error of less than 1 mmHg and 1 mL for all targeted datasets and reflect a mean heart rate [bpm] within the range of mean values noted in the reference material.

Discussion
A novel method for simulating LV-PV control functionality utilizing explicitly defined preload, afterload, and contractility was delivered for cardiovascular intervention assessment. The proposed aggregate model consists of three subsystems which include a preload, afterload, and contractility sensitive computational setpoint calculator ("PV loop critical point determination" section), a FSM controller ("PV loop modeling utilizing a state machine control architecture approach" section), and a hydraulic testing system ("Hydraulic testing model utilizing MathWorks' Simulink ® and SimscapeTM toolbox" section). The computation model provides pressure and volume setpoints based on the coefficients revealed by best fit equations for ESPVR, EDPVR, and E a . The acquired setpoints drive the FSM controller to perform the prescribed PV relationship. Then the hydraulic testing system, which reproduces conditions comparable to those found in a left heart MCS with cardiac piston actuation, simulates the PV relationship defined by the inputs to the computational model. The resulting solution was capable of being compiled for hardware control in a MCS through the architecture and solver type employed; deterministic processing is achievable and runtime setpoint changes can be made. Simulink ® and its supplemental product library was effective at developing reproducible clinical conditions, which would be determined through an individualized clinical PV analysis, simulated in silico for this work with ability to translate to future in vitro investigations. This provides an MCS with the capabilities to investigate the pathophysiology for a specific individual, with or without VAD support, by reproducing the precise PV relationship defined by their left ventricular functionality.
In silico verification of the LV-PV loop critical point computational model, FSM control architecture, and hydraulic testing system support this modeling approach as an effective means of simulating the LV-PV relationship. In this work, a novel method for simulating the characteristic EDPVR curve and LAP during diastolic filling was presented. This approach proved to be an effective means of capturing the nuisances in those sections of the PV curve that are critical for diastolic operation of mechanical circulatory support systems and not found in prior computational models [15,41].
As shown in Fig. 5a and Table 2, the computational model was able to create specific points that the FSM was able to utilize as features governing the transition between LV-PV states, given a clinical preload analysis, similar to Fig. 1b. Additionally, the hydraulic testing model was able to produce a suitable degree of realism to be able to evaluate the feasibility of this methodology, producing realistic conditions to include LAP and AoP. The delivered capabilities enable control of the PV relationship beyond that presented in prior work on elastance based control with respect to dynamic afterload response [21,24] and software-oriented control [44].
A key result of this investigation is a novel in silico method for simulating LV-PV relationships based on an analysis of a patient's ESPVR, EDPVR, and E a curves. Displayed in Fig. 6 is the characteristic LV-PV loop of two individuals presented in the research conducted by Westermann et al. [50]. Simulated is Heart Failure with Normal Ejection Fraction (HFNEF) and the Control developed by means of a preload reduction analysis and quantified by means of data capture tools. Both simulations

Conclusion
The limitations of this approach are mainly the ideal hydraulic testing system and use of anticipatory limits in transition points of the PV loop. If a force is exerted into this computational model of the hydraulic system, the system responds with the corresponding pressure instantaneously within that sample period. There was no modeled delay or rise time in the actuation components. This consideration is made in the FSM by increasing force incrementally instead of applying a constant desired force. Some parameters which define the hydraulic system, such as the parameters within the Spring-Loaded Accumulator are ideal assumptions based off a desired performance of the system. The focus of this work was on the control architecture that can be adjusted to a variety of hardware platforms through manipulation of the output signal magnitude and response characteristics. Additionally, pressure sensor feedback is ideal using this modeling approach. The sensor sampling rate was set to 512 Hz and assumed an ideal sensor with low noise. Additionally, a manual offset was made to the transition from diastolic filling to isovolumetric contraction of the system; enabling a ramping from the transition of fill to eject. Moreover, an offset was utilized in the transition from isovolumetric contraction to ejection in order to allow the pressure to slowly increase to the desired LV ESP during ejection. Future work includes a sensitivity analysis regarding resistance, compliance, and force rates. This analysis will be useful in that it will quantify the exact limitations of the hydraulic testing system as well as the range of accuracy of the FSM approach. Isolated in vitro testing of this approach will be conducted on a nested-loop hydraulic system before being incorporated into an MCS for investigating accurate cardiovascular hemodynamic considerations, such as the accuracy of pressure and flowrate sensor feedback. Additionally, what-if scenarios will be conducted on a MCS in order to create feasible scenarios to which a patient may experience.
This research will assist in producing an investigatory method and MCS control logic that will advance the medical community by improving left ventricular in vitro analysis capabilities. The ability of an MCS to be able to replicate the exact PV relationship that define the pathophysiology allows for a robust in vitro analysis to be completed. This ventricular model for ventricular function could also be coupled with aortic and left atrium computational fluid dynamics (CFD) models that require inlet and outlet conditions manifested by the left ventricle. The FSM approach is computationally efficient due to the explicit computation, and simple transition logic, which is preferential when small time steps and high iteration solvers are being employed. It was this efficiency and portability in the outcome that has made this work impactful for a variety of investigative purposes.