Finite state machine implementation for left ventricle modeling and control

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–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–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 \(\frac{{\Delta {\text{P}}}}{{\Delta {\text{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 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–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–15], improving predictive strategies [16–19], and refining computational models [41–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–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

Figure 2 presents the computational model and example developed in Simulink™ to reflect Eq. 4 through 9; utilized to find the critical points which define the initiation of each phase. Figure 2a reflects the system of equations in this example, capable of being solved in real-time. Figure 2b presents a graph of these equations, with critical points noted. For this example, based on data collected using DataThief on loop 1 of Fig. 1b: a1 = 2.9745, a0 = − 17.133, b3 = 2.6435E−5, b2 = − 4.0598E−3, b1 = 0.16687, b0 = 8.5448, c1 = − 1.7504, and c0 = 185.02. The computational system produces LV_EDP = 12.043 mmHg, LV_EDV = 105.71 mL, LV_ESP = 110.13 mmHg, LV_ESV = 42.785 mL, LV_EIRP = 10.323 mmHg, and LV_EIRV = 42.785 mL. Using these parameters, LV Stroke Volume (LV_SV) = 62.93 mL, LV Ejection Fraction (LV_EF) = 0.595, LV Stroke Work (LV_SW) = 6929.9 mmHg*mL. These values are presented in Tables 1 and 2. These coefficient values can be interchanged with clinical values for individualized PV assessment, and can be controlled over time for determining the effects of ventricular functional shifts. Utilizing DataThief [48], 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.

	Original	Preload only	Afterload only	Contractility only	Control	HFNEF
Simulation time [s]	10	10	10	10	10	10
Isovolumetric rate [N*sample/s]	225	225	225	225	500	500
Isovolumetric Contraction offset [mL]	1	1	1	1	1	1
Systolic ejection rate [N*sample/s]	2	2	2	2	9	9
Systolic ejection offset [mmHg]	7	7	7	7	20	20
Resistance	0.03	0.03	0.03	0.03	0.05	0.045
Fluid chamber capacity [mL]	517.15	517.15	517.15	517.15	517.15	517.15
Preload pressure [psi]	0.01	0.01	0.01	0.01	0.01	0.01
Pressure full capacity [psi]	10	10	10	10	6	6.6667
Initial volume of fluid [mL]	110.13	110.13	110.13	110.13	210.667	210.48
A1				[2.9745
				2.9745
				2.7245
				2.4745
	2.9745	2.9745	2.9745	2.2245	1.2407	0.99741
				1.9745
				1.7245
				1.4745
				1.2245]
A0	− 17.133	− 17.133	− 17.133	− 17.133	33.857	72.586
B3	2.6435E−05	2.6435E−05	2.6435E−05	2.6435E−05	2.6928E−07	1.4046E−05
B2	− 4.0598E−03	− 4.0598E−03	− 4.0598E−03	− 4.0598E−03	− 9.3013E−06	− 2.5351E−03
B1	0.16687	0.16687	0.16687	0.16687	0.026968	0.15836
B0	8.5448	8.5448	8.5448	8.5448	2.9515	− 0.010234
C1			[− 1.7504		[− 1.1365	[− 1.4054
			− 1.7504		− 1.1365	− 1.4054
			− 1.60848		− 1.1365	− 1.4054
	− 1.7504	− 1.7504	− 1.46656	− 1.7504	− 1.1992	− 1.3994
			− 1.32464		− 1.2619	− 1.3934
			− 1.18272		− 1.3247	− 1.3874
			− 1.0408]		− 1.3874	− 1.3814
					− 1.4501]	− 1.3754]
C0		[185.02	[185.02		[211.17	[235.76
		185.02	185.02		211.17	235.76
		190.02	170.02		211.17	235.76
		195.02	155.02		200.96	220.7
	185.02	200.02	140.02	185.02	190.75	205.63
		205.02	125.0		180.53	190.56
		210.02	110.02]		170.32	175.5
		215.02]			160.11]	160.43]

	Original	Preload only	Afterload only	Contractility only	Control	HFNEF
LVESP,i [mmHg]	110.13	110.13	110.13	110.13	126.4	140.32
LVESV,i [mL]	42.785	42.785	42.785	42.785	74.589	67.91
LVEDP,i [mmHg]	12.043	12.043	12.043	12.043	9.3686	21.522
LVEDV,i [mL]	105.71	105.71	105.71	105.71	185.81	167.75
LVEIRP,i [mmHg]	10.323	10.323	10.323	10.323	5.023	3.4517
LVEIRV,i [mL]	42.785	42.785	42.785	42.785	74.589	67.91
LVSV,i [mL]	62.93	62.93	62.93	62.93	111.22	99.84
LVEF,i	0.595	0.595	0.595	0.595	0.599	0.595
LVSW,i [mmHg*mL]	6929.9	6929.9	6929.9	6929.9	14058.3	14009.5
Simulated LVESP,i [mmHg]	110.05	110.05	110.05	110.05	125.8	140.1
Error LVESP,i [mmHg]	0.08	0.08	0.08	0.08	0.6	0.22
Simulated LVESV,i [mL]	42.82	42.82	42.82	42.82	73.97	67.19
Error LVESV,i [mL]	0.035	0.035	0.035	0.035	0.619	0.72
Simulated LVEDP,i [mmHg]	12.71	12.71	12.71	12.71	9.87	22.23
Error LVEDP,i [mmHg]	0.667	0.667	0.667	0.667	0.5014	0.708
Simulated LVEDV,i [mL]	106	106	106	106	185.8	168
Error LVEDV,i [mL]	0.29	0.29	0.29	0.29	0.01	0.25
LVESP,f [mmHg]	110.13	129.02	77.061	66.075	92.071	109.51
LVESV,f [mL]	42.785	49.134	31.667	67.953	46.92	37.021
LVEDP,f [mmHg]	12.043	16.783	12.044	12.043	6.1782	6.2607
LVEDV,f [mL]	105.71	122.84	105.71	105.702	110.41	116.64
LVEIRP,f [mmHg]	10.323	10.079	10.597	9.432	4.2242	3.0906
LVEIRV,f [mL]	42.785	49.134	31.667	67.953	46.92	37.021
LVSV,f [mL]	62.93	73.71	74.04	37.75	63.49	79.62
LVEF,f	0.595	0.600	0.700	0.357	0.575	0.683
LVSW,f [mmHg*mL]	6929.9	9509.5	5705.8	2494.3	5845.6	8719.1
Simulated LVESP,f [mmHg]	109.9	128.7	77.84	66.74	92.73	109.9
Error LVESP,f [mmHg]	0.23	0.32	0.779	0.665	0.659	0.39
Simulated LVESV,f [mL]	42.95	49.19	31.9	68.18	46.65	37.16
Error LVESV,f [mL]	0.165	0.056	0.233	0.227	0.27	0.139
Simulated LVEDP,f [mmHg]	12.73	17.21	12.75	12.75	6.711	6.624
Error LVEDP,f [mmHg]	0.687	0.427	0.706	0.707	0.5328	0.3633
Simulated LVEDV,f [mL]	106.1	123.3	106.2	106.2	110.6	116.7
Error LVEDV,f [mL]	0.39	0.46	0.49	0.498	0.19	0.06
ΔLVSV [mL]	0.00	10.78	11.12	− 25.18	− 47.73	− 20.22
ΔLVEF	0.000	0.005	0.105	− 0.238	− 0.024	0.087
ΔLVSW [mmHg*mL]	0.0	2579.6	− 1224.1	− 4435.7	− 8212.7	− 5290.5
Initial Heart Rate [bpm]	59.36	59.36	59.36	59.36	59.71	63.54
Final Heart Rate [bpm]	59.94	54.52	43.73	79.69	89.43	70.7
Mean Heart Rate [bpm]	59.9	57.62	49.28	65.04	73.7	65.36

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 LV_ESP, LV_ESV, LV_EDV, LV_EIRP, time (t), simulated pressure (P), and simulated volume (V). The Output is the output variables of the model, which is Force (F) applied to the piston in Newtons, Cycle_Count, and Heart_Rate [bpm]. The organization of the state transition diagram follows FSM convention: the single curved arrow donates the initial time-dependent conditions of the model, the oval shapes are the states of the model, the dotted hoop arrows denote the output of the state until a specific condition is met, and the straight arrows are the transition direction once the condition annotated is satisfied. Time (t) is an input variable that discretely changes at the Fundamental Sampling Time of the simulation, \(\frac{1}{1024}{\text{s}}\). Correspondingly, the FSM operates at a sampling rate of 1024 Hz. After every complete cycle, the output variables Cycle_Count and Heart_Rate are calculated. Heart rate is determined based on the Cycle_Time that is updated with the current time at the initiation of Phase 1 for every cycle. 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 the rate of change in the output variable, F, during systolic ejection. Systolic Ejection Offset is defined as the value subtracted from the LV_ESP to start the initialization of the Phase 3 state, establishing LV_EICP.

Hydraulic testing model utilizing MathWorks’ Simulink^® and Simscape™ toolbox

A hydraulic testing model was developed for simulating hydraulic performance as presented in Fig. 4. This system was designed to replicate the dynamics of a force-based piston pump model that drives the pressure within a chamber between two opposing check valves. This constitutes similar conditions observed within the left ventricular portion of an MCS. The Simulink^® and Simscape™ block library provided all the necessary components needed to create a hydraulic testing platform capable of simulating this application. All modified parameter values are noted in the diagram, while any parameters not noted were left standard to the block’s original parameter values. Additionally, for any element parameter denoted as ‘Variable’, these values were not left constant for all simulations presented. The values utilized in each simulation, not explicitly declared in Fig. 4, are displayed in Table 1.

The hydraulic testing model is a one-input, four-output system. The input is the force [N] applied to the piston and is regulated by means of the Stateflow^® control architecture. The outputs are simulated left ventricular volume (LVV) [mL], simulated left ventricular pressure (LVP) [mmHg], simulated aortic pressure (AoP) [mmHg], and left atrial pressure (LAP) [mmHg]. LVP and LVV are utilized by the Stateflow^® control logic to govern state transitions while AoP and LAP are used for system fidelity and plotting purposes. The input force is applied to the Ideal Force Source block element which is then directed to an Ideal Translational Motion Sensor which converts an across variable measured between two mechanical translational nodes into a control signal proportional to position. The position signal is then converted into volume [mL] based on a piston diameter of 2 inches, thus a cross-sectional area of π × 2.54² = 20.27 cm². The input force [N] is also applied to a Translational Hydro-Mechanical Converter which converts hydraulic energy into mechanical energy in the form of translational motion of the converter output member. Two check valves (aortic and mitral), positioned in opposing directions, regulate the fluid flow direction as seen in the left ventricular section of an MCS. A Constant Volume element is positioned between the two check valves to simulate a constant volume filling chamber. A Hydraulic Pressure Sensor is positioned between the opposing check valves to monitor LVP, then outputs the simulated values to the Stateflow^® control logic.

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,\(\frac{{\Delta {\text{V}}}}{{\Delta {\text{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.

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.

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 \({\text{P}} = - 1.7504\left( {\text{V}} \right) + 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 \({\text{P}} = - 1.7504\left( {\text{V}} \right) + 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 \({\text{P}} = 2.9745\left( {\text{V}} \right) - 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.

Clinical assessment of outlined approach

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.

The Control is presented in Fig. 6a. The ESPVR curve was found to be defined by the equation \({\text{P}} = 1.2407\left( {\text{V}} \right) + 33.857\) and the EDPVR curve was found to be \({\text{P}} = 2.6928{\text{E}} - 7\left( V \right)^{3} + - 9.3013{\text{E}} - 6\left( V \right)^{2} + 0.026968\left( V \right) + 2.9515\). E_a is initially defined by the equation \({\text{P}} = - 1.1365\left( {\text{V}} \right) + 211.17\) and defined by the equation \({\text{P}} = - 1.4501\left( {\text{V}} \right) + 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 \({\text{P}} = 0.99741\left( {\text{V}} \right) + 72.586\) and the EDPVR curve was found to be \({\text{P}} = 1.4046{\text{E}} - 5\left( V \right)^{3} + - 2.5351{\text{E}} - 3\left( V \right)^{2} + 0.15836\left( V \right) + - 0.010234\). E_a is initially defined by the equation \({\text{P}} = - 1.4054\left( {\text{V}} \right) + 235.76\) and defined by the equation \({\text{P}} = - 1.3754\left( {\text{V}} \right) + 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.