Introduction. The accuracy and clinical applicability of an improved model-based system for tracking hemodynamic changes is assessed in an animal study on septic shock. Methods. This study used cardiovascular measurements recorded during a porcine trial studying the efficacy of large-pore hemofiltration for treating septic shock. Four Pietrain pigs were instrumented and induced with septic shock. A subset of the measured data, representing clinically available measurements, was used to identify subject-specific cardiovascular models. These models were then validated against the remaining measurements. Results. The system accurately matched independent measures of left and right ventricle end diastolic volumes and maximum left and right ventricular pressures to percentage errors less than 20% (except for the 95th percentile error in maximum right ventricular pressure) and all . An average decrease of 42% in systemic resistance, a main cardiovascular consequence of septic shock, was observed 120 minutes after the infusion of the endotoxin, consistent with experimentally measured trends. Moreover, modelled temporal trends in right ventricular end systolic elastance and afterload tracked changes in corresponding experimentally derived metrics. Conclusions. These results demonstrate that this model-based method can monitor disease-dependent changes in preload, afterload, and contractility in porcine study of septic shock.

1. Introduction

Cardiovascular management of the critically ill is often difficult due to the array of complex circulatory interactions that occur within hemodynamically compromised patients. The limited measurement set typically available in the intensive care unit (ICU) is sometimes insufficient to give a full picture of cardiovascular status. Hence, diagnosis and therapy are often based on the experience, skill, and intuition of the attending medical staff. However, difficulties in interpreting cardiac and circulatory measurements can result in misdiagnosis and suboptimal treatment, which may lead to inefficient use of resources, increased length of stay, and in some cases death [15].

To improve hemodynamic monitoring, model-based methods can be used to aggregate common ICU measurements into an easily understandable, physiological form. Other model-based methods have been used to estimate surrogate markers of cardiovascular performance from clinically available data. However, these approaches have only been used to identify a small number of hemodynamic parameters [6] or only describe localised behaviour [7, 8], normally to obtain estimates of CO [916]. To truly understand the effects of many common cardiovascular diseases in the ICU, one must consider the global effects these diseases have on a patient. Hence, a broader approach that describes the effects of disease on parameters of cardiac function, as well as parameters of systemic and pulmonary flows, is required. Patient-specific models, identified from clinical data, could be used to help paint a clearer physiological picture of the patient’s global cardiovascular condition and assist medical staff with decision making [17, 18]. The goal of this research is to show that such an in silico system can accurately identify important disease and treatment dependent changes in the cardiovascular system (CVS).

A previously reported model-based method [19, 20] indicated that subject-specific disease induced hemodynamic changes due to pulmonary embolism could be tracked from common cardiovascular measurements. However, the method lacked clinical applicability because invasive measurements of left and right ventricular volumes waveforms were required, which are rarely available in the ICU. Thus, an enhanced method [21, 22] has been developed that only uses clinically available or easily inferred cardiovascular measurements, so the system can be used to monitor hemodynamics in a normal, clinical setting.

In the ICU, some of the most prevalent and deadly cardiovascular disorders are severe sepsis and septic shock [1, 2] which are caused by a whole body inflammatory response to an infection, resulting in systemic vasodilation. To further prove the monitoring ability of this new method, subject-specific CVS models are retrospectively identified and validated using measurements from a porcine study on septic shock and large-pore hemofiltration therapy [23]. Hence, in this study, the ability of subject-specific CVS models to monitor the signs of sepsis, including a decrease in systemic vascular resistance [2426], an increase in pulmonary vascular resistance [26], and an increase right ventricle end diastolic elastances (RVEDV) [26], were analysed.

2. Materials and Methods

This study used cardiovascular measurements recorded during a porcine trial studying the efficacy of large-pore hemofiltration (LPHF) for treating septic shock [23]. Although the efficacy of LPHF is not of particular interest to this study, the induction of sepsis followed by LPHF provides a dynamic and clinically realistic background on which to test this model and methodology. Measurements taken from four pigs during this study provide a suitable dataset for identifying subject-specific CVS model parameters and safely validating the model-outputs. In this work, subject-specific CVS models were only identified from measurements commonly available in the ICU. Changes in the identified parameters in these models were analysed to track the effect of septic shock in the pigs. No prior assumptions were made or implemented regarding the trends in these identified parameters. The more invasive measurements recorded in the porcine study, which are not widely obtainable clinically, were used to validate the accuracy of the CVS models.

2.1. Subjects

Four pure Pietrain pigs weighing 20–30 kg were premedicated, anesthetized, and ventilated as described in [23, 27]. The animals received a 0.5-mg/kg endotoxin infusion (lipopolysaccharide from Escherichia coli serotype 0127:B8; Sigma Chemical, St. Louis, MO, U.S.A) over 30 minutes (from T0 to T30) to initiate sepsis. From T60 onwards the pigs underwent a zero balance continuous venovenous hemofiltration at a rate of 45 mL/kg/h with a 0.7 m2 large-pore (78 ) membrane with a cut off of 80 kDa (Sureflux FH 70, Nipro, Osaka, Japan) and a BaxterBM 25-BM 14 hemofiltration device (Baxter Health Care, Munich, Germany). Each animal acted as its own control in this study with baseline measurements taken prior to the endotoxin infusions (T0), representing the undiseased state of the pig. All procedures and protocols used in the porcine experiments were reviewed and approved by the Ethics Committee of the Medical Faculty at the University of Liege (Belgium).

2.2. Measurements

The pulmonary trunk was accessed via medial sternotomy, and a micromanometer-tipped catheter (Sentron, Cordis, Miami, FL) was inserted into the main pulmonary artery through a stab wound in the RV outflow tract and positioned approximately 2 cm downstream of the pulmonary valve. Aortic pressure was measured using a micromanometer-tipped catheter inserted into the descending thoracic aorta through the right femoral artery. A 7F, 12 electrode conductance micromanometer-tipped catheter (CD Leycom, Zoetermeer, Holland), was positioned in each of the left and right ventricles, so that all electrodes were in their respective cavities. Further details on the trials can be found in [23].

These catheters provided measurements of the aortic and pulmonary pressure waveforms () and the left and right ventricular pressure and volume waveforms (, , , and ). Continuous waveforms of , , , , , and (sampling frequency, = 200 Hz) were recorded from the four catheters at 30 minutes intervals (T0, T30, T60, …, T210, and T240) between initiation of anaesthesia (T0) and 4 hours subsequent (T240). At each 30 minute interval 6–12 heartbeats of these waveforms were recorded. For each set of waveforms one good heartbeat worth of data was selected. From this selected heartbeat discrete convergence set points for the parameter identification process were extracted. In this study, 34 sets of measurements were extracted from the raw data providing 34 sets of convergence set points from the four pigs. Two data sets were unusable as steady-state measurements could not be obtained.

2.3. Cardiovascular System Model

This study uses a previously validated cardiovascular system model [19, 28, 3032] defined by a set of equations describing cardiac and circulatory flow. This model describes the “lumped” or global dynamics of major regions of the CVS, ignoring smaller localised behaviour. The lumped nature allows the model to be relatively light computationally, requiring a smaller number of parameters. However, it is complex enough to capture all the clinically important hemodynamics of the CVS.

The six chamber CVS model represents the left and right ventricles, aorta, pulmonary artery, vena cava, and pulmonary vein, as illustrated schematically in Figure 1. Pressures and volumes in the chambers and flows between the chambers are defined using parameters of resistance to flow , blood inertia , and chamber elastance or stiffness . Cardiac muscle activation of the ventricles is represented using time varying elastances (driL and driR) which act as the driver functions for the CVS model. The ventricular chambers in the model are actively elastic meaning pressure changes nonlinearly with respect to volume, to model the effects of contraction and relaxation in these chambers. The noncardiac chambers are passively elastic. Thus, pressure is linearly related to volume in these chambers. Flow in and out of the ventricles is controlled via pressure-gated heart valves, and the parallel interaction between the ventricles is modelled through pericardial and septal dynamics. The combination of all these model features and parameters provides a means for accurately depicting ventricular, arterial, and venous blood pressures and volumes, and the blood flow rates through the heart valves, body, and lungs. The full model definition can be found in Hann et al. [33].

Please note that this model does not implicitly simulate any autonomic control mechanism. However, such mechanisms can be indirectly simulated through alteration of the model parameters. For example, sympathetic control of the blood pressure could be simulated by increasing the parameter in the CVS model.

2.4. Parameter Identification Method: Pig-Specific Model

The parameter identification method used for this study was described by Revie et al. [21, 22]. The advantage of this method over the previously employed technique [19, 20] is that it only requires data from measurements typically available in an ICU environment. Table 1 shows the small set of measurements necessary to identify the patient-specific parameters required for the model.

The parameter identification process was performed in sequential stages to increase convergence stability. In this process the inertances in the CVS model were ignored , as they were found to have negligible effect on the modelled flow. To simplify parameter identification, the six-chamber CVS model was divided into two submodels of the systemic and pulmonary circulations, as shown in Figure 2. These models were decoupled by removing the ventricular interaction between the circulations. The submodels were independently simulated using left and right ventricular driver functions (driL and driR) and initial estimates for parameter values. The driver functions were identified using the method described in [34, 35] and help define the ventricular systolic and diastolic function in the identified models. These driver functions play an important role in defining pig and time dependent changes in diastolic ventricular function in the subject-specific CVS models, as the parameters of ventricular filling (as shown in in the appendix) are held as constants (see Table 5) during the identification process.

An iterative proportional gain control method [22] was used to match the model outputs of Table 1 to the corresponding measurements, identifying the desired parameter set of Table 1. Each observable model output was used to identify one model parameter. Parameters () that were proportionally related to their corresponding model output were iteratively identified using the ratio of the true measurement to corresponding model output: Parameters that were inversely related to their chosen model output were identified using The parameters of the systemic submodels were identified first, followed by the pulmonary model. Once the both submodels were identified, they were joined together, with coupling between the models added in the form of septum and pericardium dynamics. The systemic and pulmonary models were then reidentified with the knowledge of this coupling, and the whole process was repeated until whole CVS model had converged. This process was repeated for each measurement set for that pig (i.e., T0, T30, T60, …, T240). From these identified models the valve resistances were averaged, and the remaining parameters were reidentified using these averaged valve resistances each measurement set from T0 to T240. Hence, for a given pig, nine unique sets of parameters were identified representing the hemodynamic state of the pig at the point in the experiment. A full description of the model identification process is provided by Revie et al. [22] and in the appendix for the interested reader. Figure 3 shows an example of the systemic model matching the mean, amplitude, and maximum ascending gradient of the aortic pressure waveform, through the identification of , , and , respectively (see , , and in the appendix). It should be noted that the parameter identification method reported in [22] has been modified, so the dead space of volumes of the ventricles is assumed to be equal to 23 mL, as reported in [36], rather than 0 mL.

In this study, global end diastolic volume (GEDV) was assumed to be equal to the sum of the left and right ventricular end diastolic volumes (LVEDV + RVEDV). Stroke volume (SV) was calculated as the average of the amplitudes of the left and right ventricular volume waveforms. The mitral and tricuspid valve closure times () were estimated from the derived subject-specific left and right ventricular driver function, (time varying elastances). These driver functions were calculated using the method described by Hann et al. [37] from SV and features in the aortic and pulmonary artery pressure waveforms.

In this study, GEDV and SV were not directly measured but instead were inferred from the experimental measurements of and . However, clinically, GEDV and SV can be measured using minimally invasive thermodilution techniques [3841]. The aortic pressure waveform can be estimated from radial artery pressure using one of the several transfer function methods [4247]. Pulmonary artery indices can be determined using a pulmonary artery catheter (PAC). Vena cava pressure, , can be identified from central venous pressure which is normally measured in the ICU, removing the need to measure . Mitral valve closure time, , can be estimated from ECG data. Hence, the data set required for model identification represents only measurements that are already measured or inferred in a typical ICU.

2.5. Analysis

Data sets were recorded every 30 minutes from T0 to T240 during the four trials. However, two sets of measurements for Pig 1 at T30 and T90 were unusable, as the pig’s hemodynamics were unstable at these times. Thus, 34 data sets were available for analysis. For each data set, a subject-specific CVS model was retrospectively identified from a subset of available measurements. The model identification method was run on a 2.13 GHz dual core machine with 3 GB of ram. Since this research is still in the development stages it has only been fully tested using the development-orientated but relatively slow MATLAB software (MathWorks, Natwick, MA, USA). Using one processor the identification method took on average 6 minutes and 24 seconds to identify a subject-specific model of the CVS. Preliminary tests using C programming language, which is better suited for real time applications, have suggested that the identification process time can be reduced by a factor of 20–100, to approximately 4–19 seconds per identified model, which is an acceptable run time in a clinical environment.

During parameter identification, all the parameters listed in Table 1 were uniquely identified for each set of measurements (at T0, T30, …, T240), except for the valve resistances which were averaged over the duration of the trials (T0–T240) for each pig. The model outputs of the identified models were compared to the remaining recorded measurements to validate the accuracy of the identified models. Averaged data is presented as mean one standard deviation (1SD). A paired-sample -test was used to check temporal variance over T0–T60 to analyse the effect of the endotoxin intervention. was considered a statistically significant result. The relationship between the modelled and measured maximum left and right ventricular pressures and volumes was analysed statistically using correlation and linear regression analysis, including calculation of absolute percentage errors, and bias and precision analysis (Bland and Altman [48]). A percentage error less than 20% was regarded as an acceptable result as measurement of physiological variables often lack precision, with errors of ±10–20% not uncommon [4951].

3. Results

The temporal variance in the measurements of the cohort was analysed between T0–T30, as well as T0–T60 and T30–T60 using paired sample -tests. Temporal variance was analysed over T0–T60 and T30–T60 because sometimes the symptoms of septic shock do not manifest immediately after the endotoxin infusion at T30. Statistically significant changes () were seen in the measured systolic and diastolic aortic and pulmonary artery pressures and right ventricular end diastolic volume over the first hour, showing the expected influence of the endotoxin infusion on the circulation. Figure 4 shows these main hemodynamic measurements across the four animals.

3.1. Subject-Specific Model Validation

For all the subject-specific models, the model outputs matched the measurements used in the identification process (mean aortic and pulmonary pressures, aortic and pulmonary artery pulse pressure, and stroke volume) to percentage errors of less than 0.5%. This result is expected and indicates the convergence of the identification process.

Validation was achieved by comparing model outputs of left and right ventricular end diastolic volumes (LVEDV, RVEDV) and maximum left and right ventricular pressures () to their corresponding measured value. These measurements were not used directly in the identification process, and thus, represent an independent validation of the method. Bias, precision, and correlation metrics were analysed for validation (Table 2 and Figures 5 and 6 for LVEDV and RVEDV). All four model outputs correlated well with the measured data (). Modelled left ventricular outputs (LVEDV, ) had a small negative bias (−3.4 mL, −1.6 mmHg), whereas right ventricular outputs (RVEDV, ) tended to slightly overestimate (3.8 mL, 4.5 mmHg) the true measurement. The precision (2 standard deviations) of the LVEDV and RVEDV predictions was 10.2 mL and 10.2 mL and 12.8 mmHg and 12.9 mmHg for and . All percentage errors were within and acceptable range (<20%), except for the 95th percentile error in the modelled maximum right ventricular pressure.

For further validation, the temporal trends of the identified right ventricular end systolic elastance, pulmonary afterload, and right ventricular arterial coupling (RVAC), from the model, were compared to corresponding metrics (, , and /). These metrics were determined experimentally in [23] using right ventricular pressure-volume analysis, for , and a windkessel model, for . The afterload on the right ventricle in the model is the sum of the resistance of the pulmonary valve and pulmonary vasculature divided by the period of one heartbeat . From Figure 7 it is seen that the subject-specific models in most cases track the experimental trends in and . Weaker relationships are seen between the modelled and experimentally RVAC metrics ( /). However, the model did capture all the averaged temporal trends for the cohort with , 94, and 0.71 for , , and /. Note that experimental data was not available for Pig 1 at T0, T30, and T90, and Pig 2 at T150, due to difficulties performing the vena cava occlusion manoeuvre.

3.2. Septic Shock Trends

Averaged model metrics were used to test whether the subject-specific CVS models could capture the general trends of septic shock across the cohort of pigs. These metrics include ventricular preload, afterload, and contractility, three of the main determinants of cardiac and circulatory function. The averaged temporal trends for the metrics used in this analysis are shown in Figure 8. Please note that the spike in some of these metrics at T30 due to initial reaction of the pigs to endotoxin insult is not discussed in the following sections. The analysis of preload, afterload, and contractility, as presented in the following, only focuses on overall trends seen in these parameters over the four-hour duration of the trials.

3.3. Preload

In the model, preload on the ventricles is represented by end diastolic volumes in the left and right ventricles (LVEDV, RVEDV). The average modelled LVEDV and RVEDV are shown in Figure 8. These results show the subject-specific models identified an increase in RVEDV during the study indicating the onset of right ventricle distension. In contrast, an initial decrease in LVEDV was observed, suggesting a decrease in venous return to the left heart immediately after the endotoxin infusion. However, LVEDV does start to increase during the second half of the trial indicating the activation of left ventricle compensatory mechanism which combats the effects of low venous return to the heart. Both these observations are known consequences of septic shock [25, 26].

3.4. Afterload

Another main hemodynamic consequence of septic shock is sudden decrease in left ventricle afterload due to a systemic inflammation [2426]. In the model, the parameter systemic vascular resistance () represents the main component of the afterload on the left ventricle. The average trends for , along with its pulmonary counterpart, pulmonary vascular resistance (), can be seen in Figure 8. These identified trends show a significant decrease in after T60, symptomatic of sepsis. Moreover, a steady increase in after T60 was observed by the subject-specific CVS models, indicating an increase in right ventricle afterload due to the effects of the endotoxin. Again, both these findings are indicative of septic shock [2426].

3.5. Contractility

In the subject-specific CVS models contractility is represented by parameters of left and right ventricular end systolic elastance (), as shown in Figure 8. In septic shock a decrease in these parameters is expected due to myocardial depression [25, 26, 52, 53]. However, contrary to this expected trend, the averaged values for and stayed relatively constant over the trial. This unexpected outcome may be the consequence of the hemofiltration therapy, which has been shown to improve cardiac function during endotoxic shock in several experimental studies [23, 54, 55]. On the other hand, these results could be due to autonomous reflex responses, such as baroreflex control, attempting to maintain blood pressure by increasing contractility [56]. As a result, such mechanisms could be counteracting the effects of decreased myocardial state in the pigs. Since autonomous control mechanisms are not directly simulated in the CVS model it is impossible to fully distinguish if changes in the model parameters are due to the effects of these mechanisms or a result of the administered therapy.

3.6. Subject-Specific Response

Although the averaged results, seen previously, show that the model is capable of identifying the hemodynamic trends of the induced disease, it is also important to know how each subject responds to the disorder and corresponding treatment. Hence, it is essential to analyse the individual time varying trends for each pig, not just the averaged trends of the disease. Figure 8 illustrates the subject-specific response of to the endotoxic shock and LPHF therapy. On average, systemic resistance drops by 20.3% across the four pigs over the duration of the trials. However, when analysing the individual response of each of the pigs, it is seen that the returns to baseline for Pig 2, and a substantial improvement is observed in Pig 1 after LPHF treatment (T60 onwards). Whereas, there are large drops in with no apparent improvement with treatment in Pigs 3 and 4, indicating the difference response to septic shock between the pigs.

4. Discussion

4.1. Modelling Errors

For validation the model was compared to independent measurements that were not used in the identification process, as well as pulmonary hemodynamic indices calculated using an independent method [23]. Almost all ventricular volumes (LVEDV and RVEDV) and maximum pressures ( and ) were identified to an error of less than 20% (except the 95th percentile error in the maximum right ventricular pressure), which is acceptable for this type of physiological monitoring. Furthermore, the time varying trends in these four measurements were identified to a high degree of accuracy, with calculated , except for the maximum right ventricular pressure where , which is still a very good result. These results suggest that subject-specific models of the CVS could be used to accurately monitor disease-dependent changes in ventricular preload in an intensive care environment.

In most cases, the temporal trends of right ventricular end systolic elastance and pulmonary afterload from the subject-specific CVS models correlated well to experimentally derived measurements of and from [23]. The modelled results tended to underestimate . The primary reason for this bias is because of the assumption in the model identification process that the dead space volume of the right ventricle () is constant and equal to 23 mL, as reported in [36]. In reality, most likely changes as the contractile state of the pigs changes during the trials, but due to the lack of available measurements it cannot be identified. Hence, in the model the changes in contractile state are solely represented by . Poorer correlations were also noticed for RVAC. The reason for these poor relationships is most likely due to accumulation of errors in calculating this compound metric (/). These include the addition of error due to measurement noise in the experimental metrics and the addition of modelling errors in the identified parameters. However, a strong temporal relationship for RVAC could be seen of when the effects of measurement noise and modelling error were averaged over the cohort (as seen in Figure 7).

Noticeably, there were larger errors in the modelled and for Pig 4, as seen in Figure 7. This pig had a substantially higher heart rate than the other pigs, on average around 150 beats per minute compared to 122, 49, and 113 beats per minute for Pigs 1, 2, and 3. As a result, Pig 4 had a very high cardiac output up to seven times higher than the other pigs. In this extreme case, the CVS model was unable to capture the high flow dynamics of Pig 4. However, work is currently underway to improve the model for these types of cases.

A primary cause of some of the larger errors and difference in polarity of the mean errors between the left and right ventricle measurements, as seen in Table 2, is the inherent difficulty of measuring and calculating the right ventricular volume experimentally with a conductance catheter. Due to the complex shape of the right ventricle, the volume measurements were underestimated, especially for Pigs 3 and 4, where the left ventricular SV appeared to be more than double the right ventricular SV. The CVS is essentially a closed system and the ventricles pump in series. Hence, it is unsustainable for their SV of the ventricles to be a largely different for more than a couple of heartbeats, during a transient period, or else there would be a large buildup of volume somewhere in the circulation. However, the measurements were taken when the pigs were hemodynamically stable and should reflect steady state conditions. Hence, it is obvious that the measured right ventricular volume is underestimating the true volume in the right ventricle.

To overcome these measurement problems, the average value of the measured left and right ventricular stoke volumes was used in the model identification process for both ventricles. Thus, once converged, the identified models generally underestimated the measured left ventricular stoke volume and overestimated the measured right ventricular stroke volume. This issue caused the LVEDV and, consequentially, the left ventricular pressure in the model to be underestimated relatives to the measurements taken in the experiment, with the opposite occurring in the right ventricle. However, more importantly, these errors were generally systematic, so the trends associated with these measurements are still clearly identified and accurate.

4.2. Detecting Septic Shock

The subject-specific models of the CVS captured the main hemodynamic trends of septic shock. A drop in systemic vascular resistance and an increase in pulmonary vascular resistance were identified. As a result, an increase in RVEDV was also seen in the subject-specific CVS models. These changes are well-known consequences of the disease [2426]. Hence, these results indicate that the identified models are able to accurately capture the expected trends of septic shock in the pigs.

4.3. Subject-Specific Modelling

The results show that personalised CVS models can be used to accurately match and predict important hemodynamic markers of afterload (), preload (LVEDV, RVEDV), and contractility (). Clinically, the contractile state and preload on the ventricles can be extremely difficult to continuously measure using traditional techniques. Generally, LVEDV and RVEDV can only be estimated intermittently using echocardiography, and there exists no widely accepted clinical measure of contractility. The model-based approach, used in this paper, provides a way to continuously estimate preload and contractility. These model-based metrics could be used to help guide therapy such as using fluids to increase LVEDV and the administering inotropes to control and . The accurate and continuous monitoring of preload is especially important in septic shock, as it has been shown that early goal-directed fluid resuscitation can increase patient outcomes [57]. Hence, metrics derived from identified subject-specific CVS models could help provide more meaningful targets for similar goal-directed therapies.

The subject-specific models were also able to track the main hemodynamic changes resulting from induced septic shock. These results show that the new model identification method utilised is capable of accurately identifying markers of cardiovascular health, like the previous integral-based method [19, 20] but requires a much smaller, more clinically available set of measurements including features (mean, amplitude, and maximum gradient) of the aortic and pulmonary artery pressures, SV, GEDV, and the mitral and tricuspid valve closure times. This measurement set is significantly smaller than the one used in the integral-based approach, in the sense that only discrete measurements are used instead of full pressure and volume waveforms. However, this new measurement set represents the most important features of these waveforms. For example, instead of fitting the integral of the aortic pressure the new method fits the mean, amplitude, and maximum gradient of aortic pressure in the CVS models, which are of more clinical importance. Hence, although less data is used to identify subject-specific CVS models, the models are actually matched to larger number of features within this data. Therefore, the smaller measurement set does not decrease the number of parameters that can be identified in these models, and its use increases the clinical applicability of the parameter identification method.

From the subject-specific CVS models, identified from this minimal set of measurements, the contrasting reaction of the pigs in response to the induced disease could be seen. Initially, decreased in all the pigs, but in later stages of the trials increased in Pigs 1 and 2, increasing aortic blood pressure in these pigs. These results tentatively suggest that the LPHF was responsible for increasing through the removal inflammatory inducing molecules from the blood stream. However, this treatment appeared to have little effect on Pigs 3 and 4. Larger trials would be required before a more substantiated conclusion on this therapy could be made, which was not the goal of this research. However, it is clear that the identified CVS models were able to accurately segregate those subjects who did (and did not) respond positively to the endotoxin insult, due either to autonomous compensatory reflexes or because of the effects of this therapy.

The parameter identification method used in this study is a significant improvement on previous work [19, 20] especially with respect to clinical feasibility. Although the measurements in this study were only taken intermittently every 30 minutes, the system has the potential to provide continuous real-time information to medical staff. Accurate model identification can be achieved using low fidelity pressure measurements, as only the main features of the aortic and pulmonary artery pressure are required such as the mean, amplitude, and maximum gradient. In practice, the aortic pressure could be inferred from the radial artery pressure using several possible methods [42, 44, 46, 47]. The method could also be implemented in the ICU at little extra cost, as it only uses measurements already taken in the ICU. Thus, the approach could be applied without adding further invasive burden to the patient. In addition, the system can be easily automated, so that a full picture of a patient’s hemodynamic state is available to clinicians with little extra effort required from medical staff.

4.4. Limitations

The number of parameters that can be identified in the CVS model is limited by the paucity of available measurements in the ICU, and thus, the amount of measurements that were used to identify pig-specific CVS models in this research. In this study only 14 out of a possible 37 parameters of the CVS model were identified. However, these 14 parameters control all the important dynamics in the model, including the contractility of the ventricles, vascular stiffness, and resistance to flow through parts of the CVS. In a clinical setting, these parameters can be controlled using common therapies including the administration of inotropes, vasoactive drugs, and/or fluid resuscitation. The other 23 parameters primarily have lesser effects, such as the inertia of blood through the heart valves, or are relatively constant over the population and are thus modelled as such, enabling accurate relative trend identification. The good validation and diagnostic results obtained in this study support the choice of model parameters to be identified.

A further limitation is the use of systemic and vascular resistance as the primary metric for ventricular afterload, although this approach is widely used [5861]. A more complicated description of right ventricular afterload was utilised for validation, as shown earlier in Figure 7, so that the magnitude and units of the metric aligned with the experimentally derived value. However, clinically, systemic and pulmonary vascular resistances are still the most commonly used definitions of afterload [58]. Hence, the systemic and pulmonary vascular resistances were chosen for examination of afterload in this study.

5. Conclusion

This study presented a method for identifying pig-specific CVS models of endotoxic shock. These models were able to identify expected disease-dependent changes in ventricular preload, afterload, and contractility in the pigs, using typically available ICU measurements. These indices of cardiac function can be difficult, impossible, or impractical to monitor clinically but are extremely useful to know when guiding cardiovascular therapy. Hence, this model-based approach could potentially be used to help monitor these determinants of cardiac function and assist clinicians with diagnosis and treatment-based decisions.


The CVS model is defined by (A.1)–(A.16). A full list of the model parameters, constants, and outputs is listed in Tables 3, 4, 5, and 6. Equation (A.1) describe the flow, , through the heart valves. These first-order differential equations state that the pressure difference across a valve is proportional to the sum of the resistive and inertial effects on the flow as follows: Nonvalvular flow is represented using Ohm’s law (A.2). Hence, the flow through the vasculature is proportional to the pressure difference across the section divided by the resistance to flow as follows: Pressure in the noncardiac chambers is proportional to the elastance () and stressed blood volume within the chamber. In these chambers the elastances are assumed to be constant over a heartbeat. The stressed blood volume in (A.3) is represented as the total blood volume () minus the unstressed volume (). Chambers within the thoracic cavity, as shown in Figure 1, are additionally influenced by the intrathoracic pressure () in the thoracic cavity as follows: The free wall pressure of the left and right ventricles (, ) is calculated using (A.4). These equations ignore the effects of the intraventricular septum and pericardium which are taken into account later with (A.11)–(A.19), hence, the name free wall pressures. Myocardial activation of the left and right ventricles are represented using normalised time varying elastance curves (driL and driR) which vary in time between 0 and 1. The first term in both of these equations represent the percentage of the pressure generation in the ventricle due to dynamic contraction of the heart muscles. At maximum contraction, driL and driR equals 1, and the pressure volume relationship in the ventricles is represented by their end systolic pressure volume relationship (ESPVR). The second term represents the percentage of pressure generation due to the passive elastic recoil of the ventricular chambers. When driL and driR equal 0, the pressures in these chambers are bound by their end diastolic pressure volume relationship (EDPVR). In between end systole and end diastole the nonlinear ventricular pressure volume relationship is defined as a weighted sum of the ESPVR and EDPVR. In other words, (A.4) represents the percentage of dynamic and passive pressure generation that is occurring in the ventricles at any point in time, as defined by The volume of blood in the noncardiac chambers varies at a rate proportional to the flow exiting and entering those chambers, as described by (A.5)–(A.8). Furthermore, (A.9) and (A.10), representing the volume in the left and right ventricles, take into consideration the nonlinear effects of valvular flow by using Heaviside functions. The use of Heaviside functions in these equations ensures that there is no backwards flow through the heart valves as follows: Finally, ventricular interaction is modelled through considering the action of the intraventricular septum and the pericardium on cardiac dynamics. Equation (A.11) represents all the pressures acting on the intraventricular septum, where the muscle activation of the septum (driS) is assumed to be equal to the average of the left and right ventricular normalised time varying elastance curves (driL, driR). From (A.11) the blood volume displaced by the septum, , can be evaluated and used to calculate the real left and right ventricular volumes (), as shown by (A.13) and (A.14) as follows: The pressure generated by the pericardium , which encases the heart, is proportional to the sum of the left and right ventricular volumes as described by (A.15) and (A.16). The total external pressure acting on the cardiac chamber , as seen in (A.17)–(A.19), is therefore the sum of the pericardium and intrathoracic pressure, as the heart is located in the thoracic cavity. Please see [28, 3032], for a more detailed description of the derivation of the CVS model equations. One has Equations (A.17)–(A.34) make up the identification process. Table 3 lists all the parameters identified in the identification process and convergence set points used to identify them. Table 4 lists the driver functions, also known as time varying elastances used to simulate the effect of cardiac contraction in the left and right ventricles in the CVS model, and Tables 7 and 8 list the initial parameters estimates used to commence the parameter identification process.

Initially, the driver functions for the left and right ventricles are estimated using GEDV and features from the aortic pressure waveform (for driL) and pulmonary artery pressure waveform (for driR) using the method described in [34, 35]. The systemic model is identified first. The model is initially simulated using the parameters values given in Table 7. A method is then used to iteratively simulate and improve the estimates of , , and in the systemic submodel by matching the modelled SV, PPao and to their measured value as follows: where , where , where .

Once the modelled , PPao, and have converged, new estimates of and are identified. , , and are then reidentified with the new estimates of and , and the whole process is repeated until does not change between iterations, and converges to its measured value. In this process, and are identified independently of , , and , as they interact and are dependent on the correct convergence of these parameters as follows: During this process is left as its constant value to be identified later in the identification process.

The same approach is used to identify the pulmonary sub-model, as the pulmonary model essentially mirrors the systemic model. , , and are identified first. Then and are updated, these values are used to reidentify , , and , and the process is repeated until SV, PP, , , and have converged. During this process is held as its original value to be identified later. Consider where , where , where Once, the systemic and pulmonary submodelss have converged, and are updated. First, the sum of the end systolic elastances () is updated by matching GEDV as follows: where .

Then a relationship, derived from the afterloads on the left and right ventricles, is used to calculate the relative contribution of from this sum. The derivation of this relationship is shown in [22] as follows: Finally, can also be found using During this step in the identification process, ventricular interaction, in the form of septal and pericardial dynamics ( and ) are calculated and added to the systemic and pulmonary submodelss to couple the systems. Once and have been updated, and and have been calculated, the systemic and pulmonary models are reidentified using the new values for these parameters and interactions. This overall process, for identifying and , is repeated until GEDV converges, thus completing the identification of the systemic and pulmonary models.

The systemic and pulmonary models are then identified for the remaining measurement sets, T0–T240 for that pig. The valve resistance are averaged across these models with these averaged valve resistance used to reidentify the systemic and pulmonary models. Hence, (A.20), (A.24), (A.25), and (A.29) are no longer required, and (A.23) and (A.28) are replaced by (A.34) and (A.35) for identifying and as follows: Once the systemic and pulmonary models have been identified with the averaged valve resistances, the final parameters of the six-chamber model, and , can be identified. This is done by using the parameters of systemic and pulmonary to simulate the six-chamber CVS model. The six-chamber model is iteratively simulated, and (A.36) and (A.37) are used to identify and . The parameter identification process is concluded when the mean of the six-chamber and matches their corresponding identified values from the simplified submodels (, ), producing a full six-chamber subject-specific model of the CVS as follows:

Figure 9 outlines the steps used to identify the parameters of the six-chamber CVS model.