Evaluation of a Model-Based Hemodynamic Monitoring Method in a Porcine Study of Septic Shock

Abstract

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 𝑅²>0.76. 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 [1–5].

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 [9–16]. 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 (𝑅_sys) [24–26], an increase in pulmonary vascular resistance (𝑅_pul) [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 m² 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 (𝑃_ao,𝑃_pa) and the left and right ventricular pressure and volume waveforms (𝑃_lv, 𝑉_lv, 𝑃_rv, and 𝑉_rv). Continuous waveforms of 𝑃_ao, 𝑃_pa, 𝑃_lv, 𝑉_lv, 𝑃_rv, and 𝑉_rv (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, 30–32] 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].

Figure 1 

Overview of the six-chamber cardiovascular system model which relates the pressures (𝑃), volumes (𝑉), and flows (𝑄) in the cardiovascular system using parameters of elastance (𝐸), resistance (𝑅), and inertance (𝐿); driL and driR represent cardiac muscle activation and act as the driver functions for the CVS model.

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 𝑅_sys 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.

Table 1 

List of identified parameters and corresponding measurements used to identify them for the simplified systemic, simplified pulmonary, and six-chamber models. Each line on the table represents a measurement/parameter pair, where the measurement on the line is used to directly identify the corresponding parameter. All the parameters in this table, except for the valve resistances, are uniquely identified for each set of measurements. The valve resistances are identified for each set of measurements and then averaged for a subject.

The parameter identification process was performed in sequential stages to increase convergence stability. In this process the inertances in the CVS model were ignored (𝐿_mt=𝐿_av=𝐿_tc=𝐿_pv=0), 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 (A.4) in the appendix) are held as constants (see Table 5) during the identification process.

(a)

(b)

Figure 2 

Simplified submodels of (a) the systemic and (b) pulmonary circulations with inertia and ventricular interaction removed. Note that for comparison the orientation of the pulmonary circulation has been reversed with respect to Figure 1.

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: 𝑃_new=Measurement𝑃ModelOutput^old. (1) Parameters that were inversely related to their chosen model output were identified using 𝑃_new=ModelOutput𝑃Measurement^old. (2) 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 𝑅_sys, 𝐸_ao, and 𝑅_av, respectively (see (A.21), (A.22), and (A.24) 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.

Figure 3 

Example of the features from the aortic pressure waveform which are matched during the model identification process.

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 (𝑡_mt,𝑡_tc) 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 𝑉_lv and 𝑉_rv. However, clinically, GEDV and SV can be measured using minimally invasive thermodilution techniques [38–41]. The aortic pressure waveform can be estimated from radial artery pressure using one of the several transfer function methods [42–47]. Pulmonary artery indices can be determined using a pulmonary artery catheter (PAC). Vena cava pressure, 𝑃_vc, can be identified from central venous pressure which is normally measured in the ICU, removing the need to measure 𝑡_tc. Mitral valve closure time, 𝑡_mt, 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. 𝑃<0.05 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 [49–51].

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 (𝑃<0.05) 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.

(a)

(b)

(c)

(d)

(e)

Figure 4 

Evolution of the averaged hemodynamic measurements recorded during the trials. *indicates 𝑃<0.05 for expected temporal changes over T0–T30, T0–T60, or T30–T60 due to the induction of septic shock. The data is presented as mean ± one standard deviation.

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 (𝑃_lv,max,𝑃_rv,max) 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 (𝑅²≥0.76). Modelled left ventricular outputs (LVEDV, 𝑃_lv,max) had a small negative bias (−3.4 mL, −1.6 mmHg), whereas right ventricular outputs (RVEDV, 𝑃_rv,max) 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 𝑃_lv,max and 𝑃_rv,max. All percentage errors were within and acceptable range (<20%), except for the 95th percentile error in the modelled maximum right ventricular pressure.

Table 2 

The median and range of the measured left and right end diastolic volumes (LVEDV, RVEDV) and maximum left and right ventricular pressures (𝑃_lv,max, 𝑃_rv,max). Bias and precision metrics (2SD), median percentage errors with 5th and 95th percentile bounds, and correlation coefficients of the modelled LVEDV, RVEDV, 𝑃_lv,max, and 𝑃_rv,max.

(a)

(b)

Figure 5 

Regression (a) and Bland-Altman analysis showing 2 standard deviation limits (b) of the modelled and measured left ventricular end diastolic volume (LVEDV).

(a)

(b)

Figure 6 

Regression (a) and Bland-Altman analysis showing 2 standard deviation limits (b) of the modelled and measured right ventricular end diastolic volume (RVEDV).

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 (𝐸_es,rvf, 𝐸_𝑎, and 𝐸_es,rvf/𝐸_𝑎). These metrics were determined experimentally in [23] using right ventricular pressure-volume analysis, for 𝐸_es,rvf, 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 (𝐸_𝑎,model=[𝑅_pv+𝑅_pul]/𝑇). From Figure 7 it is seen that the subject-specific models in most cases track the experimental trends in 𝐸_es,rvf and 𝐸_𝑎. Weaker relationships are seen between the modelled and experimentally RVAC metrics (𝐸_es,rvf /𝐸_𝑎). However, the model did capture all the averaged temporal trends for the cohort with 𝑅²=0.62, 94, and 0.71 for 𝐸_es,rvf, 𝐸_𝑎, and 𝐸_es,rvf /𝐸_𝑎. 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.

Figure 7 

Comparison of subject-specific and averaged modelled right ventricle end systolic elastance (𝐸es,rvf), pulmonary elastance (𝐸𝑎), and right ventricular arterial coupling (𝐸es,rvf/𝐸𝑎), represented by the black symbols, against experimentally derived metrics from [23], represented by the white symbols.

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.

Figure 8 

Modelled subject-specific indices of preload (end diastolic volume), afterload (vascular resistance), and contractility (end systolic elastance). The grey symbols represent the left ventricle or systemic metrics. The black symbols represent right ventricle or pulmonary metrics.

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 [24–26]. In the model, the parameter systemic vascular resistance (𝑅_sys) represents the main component of the afterload on the left ventricle. The average trends for 𝑅_sys, along with its pulmonary counterpart, pulmonary vascular resistance (𝑅_pul), can be seen in Figure 8. These identified trends show a significant decrease in 𝑅_sys after T60, symptomatic of sepsis. Moreover, a steady increase in 𝑅_pul 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 [24–26].

3.5. Contractility

In the subject-specific CVS models contractility is represented by parameters of left and right ventricular end systolic elastance (𝐸_es,lvf,𝐸_es,rvf), 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 𝐸_es,lvf and 𝐸_es,rvf 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 𝑅_sys 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 𝑅_sys 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 𝑅_sys 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 (𝑃_lv,max and 𝑃_rv,max) 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 𝑅²>0.90, except for the maximum right ventricular pressure where 𝑅²=0.76, 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 𝐸_es,rvf and 𝐸_𝑎 from [23]. The modelled results tended to underestimate 𝐸_es,rvf. 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 (𝑉_𝑑,rvf) is constant and equal to 23 mL, as reported in [36]. In reality, 𝑉_𝑑,rvf 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 𝐸_es,rvf. 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 (𝐸_es,rvf/𝐸_𝑎). 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 𝑅²=0.71 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 𝐸_es,rvf 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 [24–26]. 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 (𝑅_sys,𝑅_pul), preload (LVEDV, RVEDV), and contractility (𝐸_es,lvf,𝐸_es,rvf). 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 𝐸_es,lvf and 𝐸_es,rvf. 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, 𝑅_sys decreased in all the pigs, but in later stages of the trials 𝑅_sys increased in Pigs 1 and 2, increasing aortic blood pressure in these pigs. These results tentatively suggest that the LPHF was responsible for increasing 𝑅_sys 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 [58–61]. 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.

Appendix

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: 𝐿_av𝑄̇_av=𝑃_lv−𝑃_ao−𝑄_av𝑅_av,𝐿_mt𝑄̇_mt=𝑃_pu−𝑃_lv−𝑄_mt𝑅_mt,𝐿_pv𝑄̇_pv=𝑃_rv−𝑃_pa−𝑄_pv𝑅_pv,𝐿_tc𝑄̇_tc=𝑃_vc−𝑃_rv−𝑄_tc𝑅_tc. (A.1) 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: 𝑄_sys=𝑃_ao−𝑃_vc𝑅_sys,𝑄_pul=𝑃_pa−𝑃_pu𝑅_pul. (A.2) 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 (𝑃_th) in the thoracic cavity as follows: 𝑃_pu=𝐸_pu𝑉_pu−𝑉_𝑑,pu+𝑃_th,𝑃_pa=𝐸_pa𝑉_pa−𝑉_𝑑,pa+𝑃_th,𝑃_vc=𝐸_vc𝑉_vc−𝑉_𝑑,vc+𝑃_th,𝑃_ao=𝐸_ao𝑉_ao−𝑉_𝑑,ao. (A.3) The free wall pressure of the left and right ventricles (𝑃_lvf, 𝑃_rvf) 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𝑃_lvf=driL⋅𝐸_es,lvf⋅𝑉_lvf−𝑉_𝑑,lvf+(1−driL)⋅𝑃_0,lvf⋅𝑒^{𝜆lvf(𝑉lvf−𝑉0,lvf)},𝑃−1_rvf=driR⋅𝐸_es,rvf⋅𝑉_rvf−𝑉_𝑑,rvf+(1−driR)⋅𝑃_0,rvf⋅𝑒^{𝜆rvf(𝑉rvf−𝑉0,rvf)}−1. (A.4) 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: 𝑉̇_pv=𝑄_pul−𝑄_mt,𝑉̇_pa=𝑄_pv−𝑄_pul,𝑉̇_vc=𝑄_sys−𝑄_tc,𝑉̇_ao=𝑄_av−𝑄_sys,𝑉̇_rv=Heaviside𝑄_tc𝑄_tc−Heaviside𝑄_pv𝑄_pv𝑉̇_lv,=Heaviside𝑄_mt𝑄_mt−Heaviside𝑄_av𝑄_av. (A.5)(A.6)(A.7)(A.8)(A.9)(A.10) 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, 𝑉_spt, can be evaluated and used to calculate the real left and right ventricular volumes (𝑉_lv,𝑉_rv), as shown by (A.13) and (A.14) as follows: driS⋅𝐸_es,spt⋅𝑉_spt−𝑉_𝑑,spt+(1−driS)⋅𝑃_0,spt⋅𝑒^{𝜆spt(𝑉spt−𝑉0,spt)}−1=driL⋅𝐸_es,lvf⋅𝑉_lv−𝑉_spt+(1−driL)⋅𝑃_0,lvf⋅𝑒^{𝜆lvf(𝑉lv−𝑉spt)}−1−driR⋅𝐸_es,rvf⋅𝑉_rv+𝑉_spt−(1−driR)⋅𝑃_0,rvf⋅𝑒^{𝜆rvf(𝑉rv+𝑉spt)}−1,driS=driL+driR2,𝑉_lvf=𝑉_lv−𝑉_spt,𝑉_rvf=𝑉_rv+𝑉_spt. (A.11)(A.12)(A.13)(A.14) The pressure generated by the pericardium (𝑃_pcd), which encases the heart, is proportional to the sum of the left and right ventricular volumes (𝑉_pcd) as described by (A.15) and (A.16). The total external pressure acting on the cardiac chamber (𝑃_peri), 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, 30–32], for a more detailed description of the derivation of the CVS model equations. One has𝑃_pcd=𝑃_0,pcd𝑒^{𝜆pcd(𝑉pcd−𝑉0,pcd)},𝑉−1_pcd=𝑉_lv+𝑉_rv,𝑃_peri=𝑃_pcd+𝑃_th,𝑃_lv=𝑃_lvf+𝑃_peri,𝑃_rv=𝑃_rvf+𝑃_peri. (A.15)(A.16)(A.17)(A.18)(A.19) 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.

Table 3 

Identified CVS model parameters and the measurements (convergence set points) used to identify them.

Table 4 

Identified cardiac driver functions (normalised time varying elastances) and the measurements used to calculate them.

Table 5 

CVS model constants (from [28, 29]).

Table 6 

Full list of the CVS model outputs. Check marks indicate outputs that are also experimentally measured.

Table 7 

Initial systemic model parameter inputs.

Table 8 

Initial pulmonary model parameter inputs.

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 𝑅_mt, 𝐸_ao, and 𝑅_sys in the systemic submodel by matching the modelled SV, PP_ao and 𝑃_ao,mean to their measured value as follows: 𝑅_mt,new=SV_lv,approxSV_true𝑅_mt,old, (A.20) where SV_lv,approx=max(𝑉_lv)−min(𝑉_lv), 𝐸_ao,new=PP_ao,truePP_ao,approx𝐸_ao,old, (A.21) where PP_ao=max(𝑃_ao)−min(𝑃_ao), 𝑅_sys,new=𝑃_ao,mean,true𝑃_{ao,mean,approx}𝑅_sys,old, (A.22) where 𝑃_ao,mean=(max(𝑃_ao)+min(𝑃_ao))/2.

Once the modelled SV, PP_ao, and 𝑃_ao,mean have converged, new estimates of 𝑃_pu and 𝑅_av are identified. 𝑅_mt, 𝐸_ao, and 𝑅_sys are then reidentified with the new estimates of 𝑃_pu and 𝑅_av, and the whole process is repeated until 𝑃_pu does not change between iterations, and 𝑑𝑃_ao,max/𝑑𝑡 converges to its measured value. In this process, 𝑃_pu and 𝑅_av are identified independently of 𝑅_mt, 𝐸_ao, and 𝑅_sys, as they interact and are dependent on the correct convergence of these parameters as follows: 𝑃_pu=𝑃_lv𝑡_mt,𝑅_av,new=𝑑𝑃_{ao,max,approx}/𝑑𝑡𝑑𝑃_ao,max,true/𝑑𝑡𝑅_av,old. (A.23)(A.24) During this process 𝐸_es,lvf 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. 𝑅_tc, 𝐸_pa, and 𝑅_pul are identified first. Then 𝑃_vc and 𝑅_pv are updated, these values are used to reidentify 𝑅_tc, 𝐸_pa, and 𝑅_pul, and the process is repeated until SV, PP, 𝑃_ao,mean, 𝑃_vc, and 𝑑𝑃_ao,max/𝑑𝑡 have converged. During this process 𝐸_es,rvf is held as its original value to be identified later. Consider 𝑅_tc,new=SV_rv,approxSV_true𝑅_tc,old, (A.25) where SV_rv,approx=max(𝑉_rv)−min(𝑉_rv), 𝐸_pa,new=PP_pa,truePP_pa,approx𝐸_pa,old, (A.26) where PP_pa=max(𝑃_pa)−min(𝑃_pa), 𝑅_pul,new=𝑃_pa,mean,true𝑃_{pa,mean,approx}𝑅_pul,old. (A.27) where 𝑃_pa,mean=(max(𝑃_pa)+min(𝑃_pa))/2𝑃_vc=𝑃_rv𝑡_tc𝑅_pv,new=𝑑𝑃_{pa,max,approx}/𝑑𝑡𝑑𝑃_pa,max,true/𝑑𝑡𝑅_pv,old (A.28)(A.29) Once, the systemic and pulmonary submodelss have converged, 𝐸_es,lvf and 𝐸_es,rvf are updated. First, the sum of the end systolic elastances (𝐸_es,sum) is updated by matching GEDV as follows: 𝐸_es,sum,new=GEDV_approxGEDV_true𝐸_es,sum,old, (A.30) where GEDV_approx=max(V_lv)+max(V_rv).

Then a relationship, derived from the afterloads on the left and right ventricles, is used to calculate the relative contribution of 𝐸_es,lvf from this sum. The derivation of this relationship is shown in [22] as follows: 𝐶_𝐸=𝑃_ao,mean,true−𝑃_vc𝑃_ao,mean,true+𝑃_pa,mean,true≈𝐸_es,lvf𝐸_es,lvf+𝐸_es,rvf𝐸=const,_es,lvf=𝐶_𝐸,mean𝐸_es,sum (A.31)(A.32) Finally, 𝐸_es,rvf can also be found using 𝐸_es,rvf=𝐸_es,sum,new−𝐸_es,lvf. (A.33) During this step in the identification process, ventricular interaction, in the form of septal and pericardial dynamics (𝑉_spt and 𝑃_pcd) are calculated and added to the systemic and pulmonary submodelss to couple the systems. Once 𝐸_es,lvf and 𝐸_es,rvf have been updated, and 𝑉_spt and 𝑃_pcd have been calculated, the systemic and pulmonary models are reidentified using the new values for these parameters and interactions. This overall process, for identifying 𝐸_es,lvf and 𝐸_es,rvf, 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 𝑃_pa and 𝑃_vc as follows: 𝑃_pu,new=SV_lv,approxSV_true𝑃_pu,old,𝑃_vc,new=SV_rv,approxSV_true𝑃_vc,new. (A.34)(A.35)Once the systemic and pulmonary models have been identified with the averaged valve resistances, the final parameters of the six-chamber model, 𝐸_vc and 𝐸_pu, 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 𝐸_vc and 𝐸_pu. The parameter identification process is concluded when the mean of the six-chamber 𝑃_vc,6 and 𝑃_pu,6 matches their corresponding identified values from the simplified submodels (𝑃_vc,simple, 𝑃_pu,simple), producing a full six-chamber subject-specific model of the CVS as follows: 𝐸_vc,new=𝑃_vc,simplemean𝑃_vc,6𝐸_vc,old,𝐸_pu,new=𝑃_pu,simplemean𝑉_pu. (A.36)(A.37)

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

Figure 9 

Step-by-step overview of the parameter identification used to create subject-specific CVS models.