Figures
Abstract
Intensive care medicine is complex and resource-demanding. A critical and common challenge lies in inferring the underlying physiological state of a patient from partially observed data. Specifically for the cardiovascular system, clinicians use observables such as heart rate, arterial and venous blood pressures, as well as findings from the physical examination and ancillary tests to formulate a mental model and estimate hidden variables such as cardiac output, vascular resistance, filling pressures and volumes, and autonomic tone. Then, they use this mental model to derive the causes for instability and choose appropriate interventions. Not only this is a very hard problem due to the nature of the signals, but it also requires expertise and a clinician’s ongoing presence at the bedside. Clinical decision support tools based on mechanistic dynamical models offer an appealing solution due to their inherent explainability, corollaries to the clinical mental process, and predictive power. With a translational motivation in mind, we developed iCVS: a simple, with high explanatory power, dynamical mechanistic model to infer hidden cardiovascular states. Full model estimation requires no prior assumptions on physiological parameters except age and weight, and the only inputs are arterial and venous pressure waveforms. iCVS also considers autonomic and non-autonomic modulations. To gain more information without increasing model complexity, both slow and fast timescales of the blood pressure traces are exploited, while the main inference and dynamic evolution are at the longer, clinically relevant, timescale of minutes. iCVS is designed to allow bedside deployment at pediatric and adult intensive care units and for retrospective investigation of cardiovascular mechanisms underlying instability. In this paper, we describe iCVS and inference system in detail, and using a dataset of critically-ill children, we provide initial indications to its ability to identify bleeding, distributive states, and cardiac dysfunction, in isolation and in combination.
Author summary
A common challenge clinicians face across different disciplines is estimating the hidden physiological state of a patient based on partially observed data. Here we describe iCVS (inferring Cardio-Vascular States): a dynamical mechanistic model of the cardiovascular system. We developed iCVS with a translational goal in mind, showing high explanatory power, its inference relies only on routinely available signals, and enables the identification of various clinically important shock states. We demonstrate the use of the model on a dataset that was collected in a pediatric intensive care unit.
Citation: Ravid Tannenbaum N, Gottesman O, Assadi A, Mazwi M, Shalit U, Eytan D (2023) iCVS—Inferring Cardio-Vascular hidden States from physiological signals available at the bedside. PLoS Comput Biol 19(9): e1010835. https://doi.org/10.1371/journal.pcbi.1010835
Editor: Alison L. Marsden, Stanford University, UNITED STATES
Received: December 26, 2022; Accepted: July 20, 2023; Published: September 5, 2023
Copyright: © 2023 Ravid Tannenbaum et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: The authors confirm that all data underlying the findings are fully available without restriction. The code used in the analyses is available at: https://github.com/shalit-lab/ICVS.
Funding: This work was supported by ISF (award number 1950/19 to US) and VATAT - council for higher education (to US and DE). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
Clinical challenge
Intensive care medicine is complex, resource demanding and expensive. High-stakes decisions need to be made rapidly based on the evolving clinical state of the patient and the interpretation by clinicians of continuous physiological and ancillary data. The patient’s cardiovascular state may fluctuate over minutes with rapid transitions due to external perturbations, internal failures of the physiological sub-components, or their control. Moreover, these patients are very diverse in their physiological and disease processes. Inferring the underlying hidden physiological state of a patient from observed data is a critical problem that is encountered by clinicians across different disciplines daily. While the modality of the available data and the frequency of sampling may vary between different clinical scenarios, in almost all cases clinicians must rely on partial observations that are a peripheral reflection of the internal, hidden states that are of importance. Specifically, clinicians use observable measurements of heart rate, arterial and venous blood pressures as well as findings from the physical examination (capillary refill, extremity temperature, etc.) to estimate internal variables that are not readily observable at the bedside such as cardiac output, systemic vascular resistance, filling pressures and volumes, autonomic tone, and more. As the physiological signals are often noisy and ambiguous, creating these mental models is a notoriously hard task that requires clinical expertise. Moreover, since the patient’s state fluctuates frequently, repeat assessments by clinicians are required. The clinicians then use this mental model of the patient’s cardiovascular internal state to derive the causes of instability such as shock type (hypovoloemic, cardiac or distributive for example), and choose appropriate treatments and interventions.
Modeling background
The large amount of data available from ICUs motivated researchers from computational disciplines to develop mathematical tools which aim to support decision-making in this area. An appealing approach to the analysis of continuous physiological times series relies on mechanistic modeling, grounded in a biophysical understanding of systems. Mechanistic modeling of the cardiovascular system often requires an understanding of the cardiac cycle as a pump coupled to the arterial and venous conductance systems with associated autonomic nervous control. This understanding is derived from years of experiments analyzing human and animal cardiovascular physiology. Such models are usually based on sets of ordinary differential equations (ODEs). In a sense, they simulate the thought process of the critical-care physician at the bedside trying to estimate the underlying pathophysiology of the patient to tailor care. However, due to the complexity of even the simplest models and the number of unobserved variables and unknown parameters, such mechanistic models were mostly published in the context of well-curated “toy” datasets and in-silico simulations ([1–5]; barring a few exceptions. See also the thorough review by Chase et al [4]).
Although simplified mechanistic models of the cardiovascular system show promise, so far, to the best of our knowledge, there are none that have been developed specifically for, and tested on data collected from critically-ill patients at the bedside, nor validated in real-life settings. There are several challenges standing in the way of using mechanistic models at the ICU bedside based on real-time data; these have prevented their deployment so far, and also stood in the way of using them to gain insights into the physiology of critical illness. First, as noted above, mechanistic models require estimating multiple hidden physiological parameters and variables: inverse problems are notoriously difficult, and their solutions are often non-unique. To overcome such difficulties, prior works reduced the number of free parameters by assuming constant values for some of the parameters, often taken from the literature to describe a “typical subject” [2, 5–8]. However, physiological properties can vary considerably between subjects, and attempting to use a “one size fits all” parameter often fails to properly characterize most patients. An example where this problem is particularly acute is in pediatric critical care units, where patient weight and normative values can vary significantly [9–11]. Thus, for a model to be useful at the bedside, minimal apriori assumptions regarding the values of hidden parameters must be made. Even so, most published models are still quite complex and require estimation of many hidden variables, and thus rely on sampling of multiple invasive variables that are not readily available outside of a dedicated, specialized laboratory. Many of these models are extremely complex, as they attempt to capture cardiovascular dynamics at second or even sub-second resolution, containing multiple coupled differential equations and hidden parameters rendering them nonidentifiable. On the other hand, limiting model complexity to reduce the number of estimated parameters and variables entails a reduction in the explanatory power of such models. Of note, there are several commercially available tools for bedside clinical use that attempt to estimate hidden cardiovascular variables such as cardiac output or systemic vascular resistance that rely, at least to some extent, on simple and crude mechanistic assumptions and models (as reviewed in [12, 13] or tested in children in [14]). Some of these tools are based on analysis of the arterial pulse wave and either utilize norms derived from published population data, or, calibration based on dilution methods. These tools are mostly used in adult care and have not entered common practice in pediatric critical care, partly due to concerns regarding their accuracy stemming from some of the challenges detailed above. Moreover, the models these tools rely on, if any, are usually extremely simple and do not account for any internal regulatory mechanisms, nor do they incorporate effects over timescales longer than several seconds. Thus, there is a delicate balance as a result of this trade-off between model complexity and estimability using readily available observations.
Second, the measurement frequency of physiological variables in the clinical world varies from continuous monitoring of vital signs to imaging techniques which are conducted once a day or even less. For example, in the ICU, blood pressure and arterial blood oxygen saturation are continuously monitored, while additional indirect estimation of cardiac function or volumes using tools such as echocardiograms are usually performed infrequently, every several hours to days. Models which rely on low-frequency measurements can be good for chronic conditions such as heart failure [2]. However, to estimate the state of unstable patients, mechanistic models should rely only on measurements with high temporal resolution. Of note, up until recently, [15] even these measurements that are commonly captured by the patient’s bedside monitor were not available for offline analysis or real-time ingestion for model inference.
An additional barrier to adapting cardiovascular mechanistic models to bedside use lies in the need to account for various modes of autonomic nervous control: most previously published work focused on a single, simple feedback loop, neglecting other common autonomic modulations.
Our approach: iCVS
Our motivation is translational in nature—develop an inference model that can estimate the hidden cardiovascular state, its control, and any causes for instability, in real-time and using only measurements that are readily available at the bedside for critically-ill patients. Such a model can be used both as a clinical support tool at the bedside, and to foster an understanding of the (patho-)physiology of critical illness. We derive a simple model (iCVS), which retains a high explanatory power of the cardio-vascular system including cardio-vascular control, which can capture different clinical scenarios such as bleeding or evolving hypovolemia, cardiogenic shock and distributive shock, in isolation or in combinations. iCVS takes into account the autonomic control of the cardiovascular system, which couples the dynamics of its different components, mathematically imposing coupling constraints on inferred parameters. This reduces the effective dimensionality of the parameter space over which we look for a solution to an inverse problem, allowing us to infer a larger set of unobserved parameters compared to previous models. The estimated parameters enable us to capture individual patients’ traits.
iCVS is designed to allow deployment at the bedside in the pediatric and adult intensive care unit: no prior assumptions regarding the physiological parameters of the patients are required except age and weight, and the only inputs which are required for the estimation process are arterial and venous blood pressures waveforms, which are commonly measured in the ICU. In order to gain more information about the cardiovascular state, without increasing model complexity, we introduce a novel approach that exploits both slow and fast timescales of the blood pressure trace. From the pulse contour shape of the blood pressure waveform, we extract quantities related to the peripheral resistance, stroke volume, and pulse pressure, while the main mathematical model and its inference run using the mean, smoothed, arterial and venous pressures at the longer, clinically-relevant, time scale of minutes. Thus, the main model and its inference focus on processes that evolve over minutes, and are of interest to clinicians, and further allow for a much simpler model than those that capture all the intricacies of each cardiac cycle [2]. All this is done while still capitalizing on information that can be extracted only from the arterial pressure waveform.
In what follows, we describe the iCVS model and inference system in detail and demonstrate its use and power for real-time estimation using a dataset of 10 patients admitted to a pediatric intensive care unit. While our estimation process entails inferring the full iCVS model including 15 hidden parameters, we focus our demonstration on the identification of three clinical shock states, representing a dire, hard clinical challenge, and a testable task. Shock, due to any cause is a state in which the cardiovascular system fails to deliver an adequate supply of oxygen to the tissues. There are several prototypical physiological shock states, each requiring specific interventions, without which there is a high risk for mortality or morbidity. Identifying shock and its causes is a notoriously hard task due to the nature of the observed physiological signals (demonstrated below), requiring continuous vigilance and clinical expertise. Specifically, using iCVS, we map certain hidden parameter ranges to identify (i) ongoing bleeding, (ii) distributive shock, and (iii) cardiogenic shock due to reduced cardiac contractility.
Results
iCVS: Model overview
The goal of this work is to derive a new model of the cardiovascular system which can be deployed at the patient’s bedside, requires only routinely acquired measurements as input, and provides an inference of the hidden, internal physiological state of the patient in real-time. The model should be flexible enough to capture the different types of shock and account for the inter-subject variability which characterizes the patients in the pediatric and adult critical care unit. The full details of the model and its development are provided in the Materials and Methods section.
Mechanistic models of the cardiovascular system are developed from first principles, tie together changes in blood pressures and volumes, and usually use the latter as the dependent dynamical variables in the set of ODEs [5]. We derive the model in terms of intravascular pressures since arterial and venous blood pressures are directly measurable (in contrast to blood volume). Briefly, our proposed model contains a one-chamber heart that acts like a pump [5], an arterial component and a venous component, and further assumes the peripheral organs can be lumped as a linear resistor (see Fig 1). In addition and importantly, the model contains three regulatory units which are described below. The two dynamical equations which describe the time-dependent evolution of the intravascular pressures are written as follows (the full derivation can be found in Materials and methods): (1) (2) Here Pa and Pv are the mean arterial and venous pressures respectively, Hr is the heart rate, R is the systemic vascular resistance and Pp is the pulse pressure. In addition, Ca and Cv are the arterial and venous compliances respectively, ΔVv0 is the difference between minimal and maximal unstressed venous volume, Iex is the intravascular volume change, and Stot serves as a general marker for the activation of the autonomic nervous system (see below).
The model is a serial circuit and contains one heart chamber which acts like a pump, organs, and arterial and venous compartments. Heart rate, contractility, vascular resistance and unstressed venous volume are affected by the baroreflex (Sb) which is determined by the arterial pressure, and by S which is an independent component that does not depend on the current physiological state. In addition, intra-vascular volume can change through Iex and the vascular resistance is also modulated by MSVR. Created with BioRender.com.
Our proposed model contains fifteen hidden parameters (see Table 1), and five time-dependent observables: mean arterial pressure, mean venous pressure, heart rate, peripheral resistance, and pulse pressure (Table 2). We later show how the mechanistic model can be used to estimate the hidden parameters based on the observables, and identify clinical states based on the observed data (Fig 2).
A. Arterial and venous pressure waveforms are recorded via arterial and central venous catheters that are connected to pressure transducers. Specifically, in this paper waveforms are sampled from critically-ill children hospitalized in an intensive care unit and recorded at 125Hz. B. The raw data is analyzed and five measurements are extracted—the observables. C. The observables are used as an input for the full model (iCVS) estimation and using constrained nonlinear optimization we obtain the set of parameter values that best fit observables. Created with BioRender.com.
Autonomic regulation of the cardiovascular system in iCVS.
In our model, the autonomic nervous system’s activation affects heart rate, cardiac contractility, unstressed venous volume, and peripheral resistance. It consists of two components: (i) The baroreflex (Sb)—the internal feedback loop that modulates heart rate, contractility, unstressed venous volume, and the systemic vascular resistance (R), meant to keep the mean arterial blood pressure as close to a set-point (Pset) as possible. The baroreflex represents the autonomic component that is dependent on the patient’s current cardiovascular state. (ii) An independent component (S), representing the sum of all autonomic modulation which is independent of the cardiovascular state. While the baroreflex Sb describes the activation of different components in response to a reduction in the arterial blood pressure, and thus manifests as negative correlations between the arterial blood pressure and the heart rate (as well as other variables), the independent autonomic component S does not depend on the arterial blood pressure and thus enables the model to capture different interactions between the blood pressure and heart rate, as well as other cardio-vascular components (Fig 1). For example, this component is expected to rise when sympathetic agonists are given, or when a “central arousal state” is present due to patient agitation or in response to a noxious stimulus (a typical fight or flight response).
Non-Autonomic regulation of the cardiovascular system in iCVS.
In iCVS, in addition to the autonomic control, the systemic vascular resistance is independently modulated by a component denoted by MSVR (Fig 1). This component enables the model to describe fluctuations in the peripheral resistance which are uncoupled from the autonomic activation. This can occur for example due to maladaptive vasodilatation that occurs during septic, vasoplegic, or anaphylactic shock, or as a response to intravenous infusion of specific drugs that affect SVR such as Norepinephrine or Vasopressin. Of note, a non-autonomic regulation on the systemic vascular resistance appears also in [7], in the context of metabolic consumption during physical activity. Total blood volume can be changed by an external current Iex which can be either positive, as is the case when intravenous fluids are bolused, or negative, such as during bleeding or loss of intravascular volume due to capillary leak or severe dehydration.
iCVS can simulate different shock states.
The interplay between the fifteen parameters of the model can generate different clinical scenarios. In this work we focus on three different shock states which are described by three hidden parameters. Briefly, shock is a failure of the cardiovascular system to supply adequate blood and oxygen to the body tissues. In order to describe the hypovolemic shock state, we use Iex (intra-vascular volume change), which is negative in case of intravascular fluid or blood loss. A distributive shock state is equivalent to a negative MSVR (non-autonomic vascular resistance modulation). In the next section, we also discuss the cardiogenic shock state which is caused by a reduction in K (heart contractility).
Fig 3 presents examples of different shock conditions of artificial neonate patients which are generated by the mechanistic model. The internal parameters of the subjects are set and two of them (Iex and MSVR) are presented (Fig 3A, 3B, 3I, 3J, 3Q and 3R). The observable parameters (Fig 3C–3G, 3K–3O and 3S–3W) are simulated according to the iCVS model (see Material and methods). The baroreflex, which depends on the current arterial pressure and affects the observables is also presented (Fig 3H, 3P and 3X).
Each column presents a simulation of one artificial patient with a given shock state: A-H—hypovolemic shock state, I-P—distributive shock state and Q-X—combined hypovolemic and distributive shock state. The hidden parameters are fixed, and the observables are simulated accordingly (see Methods). Two of the hidden parameters are presented for each patient: A, I, Q Time dependent intra vascular volume change (Iex). In A, Q liquid withdrawal occurs over a period of 15 minutes. B, J, R Time course of MSVR. In J, R a reduction in MSVR occurs over a period of 15 minutes. C-H, K-P, S-X—the resulting observables: C, K, S Arterial pressure (Pa), D, L, T Venous pressure (Pv), E, M, U Peripheral resistance multiplied by arterial compliance (RC), F, N, V Heart rate (Hr), G, O, W Pulse pressure (Pp). H, P, X The time dependent magnitude of the baro-reflex (Sb).
Panels A-H in Fig 3 present a simulation of a patient in a hypovolemic shock state, as can happen with ongoing bleeding. The intravascular volume change (Iex) is negative (Fig 3A), representing a reduction of intravascular volume over time, while the non-autonomic SVR modulation (MSVR) remains constant (Fig 3B). As expected, we see that the arterial and venous pressures decrease (Fig 3C and 3D). The heart rate and the vascular resistance reflexively rise (Fig 3E and 3F) due to baroreflex activation which acts to maintain blood pressure and compensates for the reduction in the intravascular volume (panel H). Note, that despite the baro-reflex activation that increases the contractility, the pulse pressure (Fig 3G) decreases due to the change in venous pressure and return (see also Eq 22 in Material and methods).
Panels I-P in Fig 3 present a simulation of a patient with a distributive shock state. The parameter Iex is zero (Fig 3I), reflecting constant intravascular volume, while MSVR decreases (Fig 3J), reflecting a reduction in the peripheral resistance which is not caused by the autonomic nervous system. It is possible to see in Fig 3N that the heart rate increases (due to the baroreflex activation, see panel P), the pulse pressure increases as well (due to the baroreflex and rise in venous pressure), and the peripheral resistance decreases due to the reduction in MSVR (Fig 3M).
Panels Q-X in Fig 3 present a simulation of a patient with a combined hypovolemic and distributive shock state. The reduction in the arterial blood pressure (Panel S in Fig 3) is more pronounced than in either shock state separately.
Of note, in addition to simulating shock states, iCVS can emulate common clinical behaviors, such as a fight-or-flight response that is elicited by pain or a noxious stimulus, mediated by autonomic tone modulation (see S1 Fig).
Estimation pipeline of the cardiovascular parameters from bedside-acquired blood pressure waveforms
The estimation pipeline is schematically illustrated in Fig 2 and is described below and in Material and Methods. Briefly, the first step includes the initial processing of the raw data and extraction of five observables. In the next step, the observables serve as an input to an optimization function which is based on iCVS and we estimate the hidden parameters and the full model.
Extracting the observables from the raw data.
As a first step, we use waveforms of venous and arterial blood pressure tracings, acquired at a high frequency that allow extraction of the following five observables: mean arterial pressure, mean venous pressure, heart rate, peripheral resistance, and pulse pressure (see Material and methods). We then use both slow and high-frequency features of the data (Fig 4). From the slow time scale we calculate the mean arterial and venous pressures. From the fast timescale (the scale of a single heartbeat), we extract: heart rate (timing of consecutive beats), peripheral vascular resistance (up to a constant, as detailed in Material and methods), and pulse pressure. These quantities, which cannot be derived without the waveforms recorded at high-frequency (a sampling rate which is much higher than the heart rate and can capture the shape of a single beat), are then smoothed and used for model inference at the slower timescale of minutes. The peripheral resistance is extracted using a non-calibrated pulse contour analysis method relying on assumptions of a simplified Windkessel model and rectangular ejection pattern during systole that has been described in detail elsewhere [16, 17].
A,B. An example of the raw data, sampled at 125 Hz frequency by the bedside patient monitor. A. Arterial blood pressure, B. Venous blood pressure. Note how these real-life data are noisy, fluctuating, and riddled with artifacts related to patient movement and care. C-G. Features that are extracted from the raw data. C,D. Mean arterial and venous pressures (slow time scale features), E-G. Heart rate, peripheral resistance and pulse pressure (fast time scale features).
Estimating the hidden cardiovascular parameters.
The next step in the estimation pipeline is to use the above five observables, as well as information on a patient’s age and weight, to fully estimate the hidden parameters of the iCVS model. Based on the iCVS model, we derive a system of equations that relates the time-dependent measurements to the fifteen hidden physiological parameters (see Material and methods). During the estimation procedure, a set of optimal parameter values that yield the best solution for the time-dependent equations is chosen by minimizing a constrained nonlinear multivariate cost function (see more details in Material and methods). For the illustrations shown below, the estimation procedure is independently conducted on partly overlapping data segments of 300 seconds. We chose a segment size of 300 seconds after an empirical search for a window that averaged out noisy and artifactual fast fluctuations, while still capturing changes over clinically-relevant timescales of several minutes.
Demonstrating iCVS estimation of cardiovascular shock states on real-world data from critically-ill patients
In this section, we illustrate the estimation process for a sample of 10 critically-ill children, hospitalized in an intensive care unit (ICU) at a large, tertiary children’s hospital. These children were admitted to the ICU either for routine observation post-operatively (for the control samples) and were stable hemodynamically, or presented with various shock states, for example, due to massive post-operative bleeding, cardiac dysfunction, infections or other causes, see Table 3 for more information. The dataset which we used to test the model’s performance is detailed in [15]. It contains physiological waveforms recorded from critically-ill children at sampling rates of 125Hz for all invasively-recorded pressures; however, any high-frequency recording that contains the full arterial blood pressure waveform is adequate for the model’s purposes. In addition, the times of all medical interventions, including drug consumption and fluid administration are recorded. Note how these real-life data are noisy, fluctuating and riddled with artifacts related to patient movement and care (see Fig 4A and 4B).
For each patient, age, weight, cause of admission to the intensive care unit and the cardio-vascular state as determined by an expert clinician are given.
As noted above, iCVS is able to deal with patients of varying weights from neonates (a couple of kilograms) to adult-sized patients, as we demonstrate below. The hemodynamic data for all patients in this ICU are collected as routine practice [15, 17]. The diagnostic labels for shock (either hypovolemic/ hemorrhagic, distributive, cardiogenic, or a combination of the above) were derived from a prospective and retrospective review of at least two expert clinicians. These labels were also corroborated using objective measurements from the medical record, when available. For example, from quantification of blood loss postoperatively, assessment of cardiac function or systemic vascular resistance using clinical examination, drugs administered, and ancillary diagnostic tests.
The challenge of estimating the causes for cardiovascular instability and shock from real-world data can be appreciated by examining Figs 5–8. Examination of panels A-E in each figure, which depict the five observables, shows that patients’ data are “messy” with marked fluctuations and trends arising from multiple physiological sources (due to the shock cause and others), unlike the clean simulated data shown before. Moreover, it is evident that simple rules linking changes in these observables to the shock cause may be misleading.
iCVS results for patient number 1 (see Table 3) are presented. A-E.—The observables: mean arterial pressure (A), mean venous pressure (B), heart rate (C.), RC—the peripheral resistance multiplied by the arterial compliance (D) and pulse pressure (E). F-H. Gray circles—optimal parameters which are achieved in each segment of 300 sec. Black line—a smoothed version of the estimated parameters (see Material and methods). F. Relative intra vascular volume change (). G. Non autonomic vascular resistance (MSVR). H. Maximal relative contractility (). Note that in this patient which suffers from a hypovolemic shock state is identified as having negative by the model.
iCVS results for patient number 5 (see Table 3) are presented. A-E. The observables: Mean arterial pressure (A), mean venous pressure (B), heart rate (C), RC—the peripheral resistance multiplied by the arterial compliance (D) and pulse pressure (E). F-H. Gray circles—optimal parameters which are achieved in each segment of 300 sec. Black line—a smoothed version of the estimated parameters (see Material and methods). F. Relative intra vascular volume change (). G. Non autonomic vascular resistance (MSVR). H. Maximal relative contractility (). Note that this patient, suffering from a distributive shock state, is correctly identified as having negative MSVR by the model.
iCVS results for patient number 4 (see Table 3) are presented. Dots mark chest opening. A-E.—The observables: mean arterial pressure (A), mean venous pressure (B), heart rate (C), RC—the peripheral resistance multiplied by the arterial compliance (D) and pulse pressure (E). F-H. Gray circles—optimal parameters which are achieved in each segment of 300 sec. Black line—a smoothed version of the estimated parameters (see Material and methods). F. Relative intra vascular volume change (). G. Non autonomic vascular resistance (MSVR). H. Maximal relative contractility (). Note that in this patient becomes non-negative and the contractility increases after the chest opening procedure, compatible with an improving in cardiac function and bleeding cessation as a result of the chest opening procedure. We note that in panel H a single point with an extremely large value is not plotted, at 9.66 minutes after chest opening ( mmHg). This point correlates with a short transient increase in the venous and arterial pressures, possibly resulting from a bolus of a vasoactive drug given to the patient at the time.
iCVS results for patient number 9 (see Table 3) are presented. A-E. The observables: mean arterial pressure (A), mean venous pressure (B), heart rate (C), RC—the peripheral resistance multiplied by the arterial compliance (D) and pulse pressure (E). F-H. Gray circles—optimal parameters which are achieved in each segment of 300 sec. Black line—a smoothed version of the estimated parameters (see Material and methods). F. Relative intravascular volume change (). G. Non autonomic vascular resistance (MSVR). H. Maximal relative contractility ().
Panels F-H in Figs 5–8 present the results of model estimation for three hidden parameters that are crucial to elucidate the causes for cardiovascular shock. A hypovolemic state is described by a negative (the intra-vascular volume change divided by the difference between the maximal and minimal unstressed venous volume ΔVvo). This magnitude does not depend on the patient weight or blood volume, has units of , and thus allows comparison between patients with different sizes. The vasodilatory state is represented by MSVR (non-autonomic vascular resistance) which is expected to be negative in a distributive shock state. The cardiac function is described by the Maximal relative contractility: , where Kmax is the maximal contractility that the ventricle can generate (see Eq 21), Cven is the compliance of the ventricle, and Ca is the arterial compliance. For constant arterial and ventricular compliances, is proportional to the maximal heart contractility. is expected to decrease in case of cardiogenic shock or cardiac dysfunction.
Each gray circle in panels F-H in Figs 5–8 corresponds to the optimal parameters which are obtained in a given time segment. The black line in each panel is a moving average of the results of different segments (see Material and methods).
Fig 5 illustrates model inference for a sample patient with post-operative bleeding. This is a neonate that presented with significant intravascular volume loss after cardiac surgery. Panels A-E show decrease in the arterial and venous blood pressures, while SVR increases and pulse pressure decreases, as expected. Of note is that for an unknown reason, the heart rate of this patient does not reflexively increase, as one would expect for a bleeding patient. Panel F shows that even without the expected heart rate increase, our model correctly identifies the negative intravascular volume change (Fig 5F). The maximal relative contractility seems stable or even increases (for comparison see the changes in for a patient with cardiac dysfunction in Fig 7).
Fig 6 illustrates the model inference results for an adolescent patient which was clinically identified as having distributive shock, specifically a vasoplegic state seen sometimes after solid organ transplantation. In the time window shown here, the mean venous and arterial pressures increase, as also the heart rate and the pulse pressure. Yet, the mean arterial pressure is low and the heart rate is high, supporting the clinical observation of a shock condition. We can further see that the peripheral resistance decreases over time. Our model indeed correctly identifies the cause for shock: MSVR which reflects the non-autonomic component of the vascular resistance is negative, in agreement with the clinically detected distributive shock. Moreover, the estimation process yields roughly zero , meaning that no intra-vascular volume change is detected, as well as demonstrating stable maximal relative contractility.
Fig 7 presents the estimation process for an unstable patient with a combined hemorrhagic, cardiogenic, and distributive shock. This was a neonate with massive bleeding post cardiac surgery, combined with left ventricular dysfunction and inadequate systemic vascular resistance. The patient was extremely unstable despite multiple interventions and underwent a chest opening procedure in the ICU (grey arrows in Fig 7) which resulted in marked improvement in cardiac function and bleeding cessation. We see that arterial blood pressure increases after the chest opening. In addition, the venous pressure decreases—compatible with an improvement in the heart function (Fig 7A and 7B). Note how the model correctly identifies the marked improvement in the maximal contractility and immediately after the procedure.
Finally, a neonate control patient that was admitted after thoracic surgery and had an uneventful post-operative course is presented in (Fig 8). As expected, the estimated intravascular volume change () and the non-autonomic SVR modulation (MSVR) are not negative, and there are no marked fluctuations in the maximal relative contractility ().
In Fig 9 we present results summarizing the model’s inference for the 10 patients. Panel A depicts the ability of our model to correctly identify intravascular volume loss: the patients marked in blue are those that had documented post-operative bleeding. Note that for these patients mean change in blood volume is more negative than their counterparts. Conversely, panel B depicts the mean non-autonomic component of systemic vascular resistance modulation. As noted above, a negative value indicates a maladaptive response with vasodilatation, as one can see in distributive shock states (for example sepsis or vasoplegia). Since the parameter Kmax is not comparable between different patients, this parameter is not presented in Fig 9. The patients marked in red were identified by expert clinicians as having abnormally low SVR. Indeed, overall, these patients were identified by our model as having lower values of MSVR. Additionally, the patients marked in pink were treated by vasoconstrictors, drugs that directly affect vascular smooth muscle to increase systemic vascular resistance, an indirect indication extracted from the electronic health record of at least a somewhat distributive state.
For each subject, the model is applied for a time window of 25 minutes (see Methods). The analysis is done over segments of 300 seconds and the average values of (panel A) and MSVR (panel B) are presented. A. The average value is presented for all sample patients. Blue—patients which are documented with post-operative bleeding, black—no post operative bleeding is documented (average value of the hypovolemic group is −0.016 min−1, average value of the non-hypovolemic group is 0.0003 min−1, Pvalue = 0.026). B. The average MSVR value of all sample patients is presented, ordered differently than in A. As there is no objective criteria for a distributive state, unlike blood loss, we resort to identifying such patients by either expert opinion or indirect evidence from the electronic medical record. In red—patients who were identified by clinical experts as suffering from distributive shock states, in pink—patients who were treated by vaso-constrictors (drugs that directly affect vascualr smooth muscle to increase systemic vascular resistance), in black—patients who were not documented with distributive shock states (average value of the distributive group (red) is −0.43, average value of the non-distributive group is −0.07, Pvalue = 0.077). Error bars in A and B represent standard error of the mean.
As an additional verification of model estimation performance on real data, we also compare between the measured observables: , and and reconstructed values for the observables which are obtained using model inference. The observables are reconstructed well by the iCVS model, with a Pearson correlation of 0.95–0.99 (see supplementary information). Furthermore, as is elaborated in the supplementary to this article, we demonstrate the robustness of our findings to the estimation algorithm.
Discussion
In this proof-of-concept work we detail a novel mechanistic model of the cardiovascular system: iCVS. We developed iCVS with a translational goal in mind—enabling estimation of the hidden cardiovascular state of critically-ill patients, in real-time, using only routinely available physiological signals at the bedside. iCVS does not preemptively assume knowledge of almost any parameter, except age and weight, maintaining flexibility for deployment in a variety of clinical scenarios. Unlike many similar models, iCVS considers the autonomic control and its specific components, imposing coupling constraints on inferred parameters and capturing common clinical phenotypes. An additional novel strength of this model is its ability to use slow and fast timescales of the arterial blood pressure waveform to gather more information about the cardiovascular state without increasing model complexity, by performing the main inference on the clinically relevant scale of minutes. In this manuscript, we detail iCVS and a suggested estimation process. This is a first modest step towards a clinically useful tool, and we provide initial indications to this extent by testing on a set of artificial patients and on data from a sample of critically ill children. We note that there are regulatory challenges facing the deployment of model-based systems at the patient bedside, including data privacy, safety, reliability and the need to continuously monitor and modify the system after deployment; however, progress has been made on this front, including the FDA’s 2021 ML-Based Software as a Medical Device (SaMD) action plan [18] and follow-up drafts for regulating ML systems at the bedside. We believe these pioneering efforts will assist the community in deploying complex systems, including systems borne out of the proof of concept model presented here.
We demonstrate the use of iCVS on a dataset of ten critically-ill children which was collected in a pediatric intensive care unit, illustrating the identification of bleeding, distributive states, cardiac dysfunction, and their combination. A key challenge of evaluating models such as iCVS on data acquired from real patients and not in a clean laboratory setting, is the lack of experimentally-defined cardiovascular hidden states. The ground-truth labels for the patients presented in this study were acquired either by information that was available from the electronic health record but which the model was not privy to (e.g. recorded blood loss), or by expert clinician assessment preformed pro- and retrospectively. In the current study, the model shows higher accuracy in the evaluation of a hypovolemic shock state than in the evaluation of a distributive shock state (see Fig 9). This may arise partly from a larger variability in the range of solutions of MSVR, (see S4–S6 Figs), and from the fact that these patients were already being treated with vasopressors, which is expected to (at least partly) offset the negative value of MSVR. In future studies, we aim to acquire data from the acute phase, before treatment is initiated.
In addition, we explore model identifiability and estimation process robustness. We did so by simulating several artificial patients with various clinical shock states and examining the model’s ability to correctly recover the preset hidden parameters, demonstrating the model’s ability to capture and track changes in these parameters (shown for three sample patients in the supplementary). We recognize that our model, like many underdetermined mechanistic physiological ones, suffers from an inherent degree of nonidentifaiblity. However, we do believe that both the demonstration with real sample patients and simulated ones indicates that our model is able to capture the salient and clinically important variations in the hidden parameters to be estimated. As noted, this is a proof-of-concept study and we plan a thorough,prospective validation and quantification of the model’s accuracy on a larger clinical set of critically-ill patients.
While we illustrated the ability of iCVS to identify hidden cardiovascular states that can not be easily identified by simple examination of the physiological recordings, it is obvious that further assessment, ideally prospectively, on a much larger patient cohort is warranted. Generation of a larger dataset, containing a diverse population of patients with various cardiovascular states and their evolution will not only allow model validation and refinement, but also to define cut-off values for shock-state identification. For example, while iCVS can identify ongoing blood loss or maladaptive vasodilatation, ideally we would wish to be able to add an additional layer for a decision support tool that would alert clinicians to the presence of vasoplegic shock or hemorrhagic shock, defined by certain thresholds relating blood pressure to MSVR or , respectively. An additional potential limitation of iCVS that can be properly quantified and addressed in future studies is its ability to deal with noisy and artefact-ridden recordings. Moreover, a prospectively collected dataset (that we are currently working towards acquiring) will also allow for precise recording of the various bedside interventions and thus test their effects on the cardiovascular hidden state and model estimation. Finally, despite its high-explanatory power, iCVS is still a relatively simple mechanistic model; it does not take into account various physiological contexts, and therefore cannot address for example the effects of positive pressure ventilation, pericardial effusion, increase in capillary wall permeability which may affect both SVR and intra-vascular volume due to fluid leakge to the tissues [19], or even interactions between the heart chambers and pulmonary and systemic circulations. An additional limitation to be acknowledged is the need for several assumptions in order to enable model simplification and estimation. We attempted to minimize their number and use only ones that are widely accepted, we stated clearly in the methods where such assumptions were made, and also cited their source, and where relevant, the literature supporting these assumptions. We also believe, that iCVS can be iteratively developed to include modules that address these limitations, potentially with addition of data-driven machine-learning modules as detailed below. An exciting avenue for future exploration is to combine iCVS with machine learning models to create robust hybrid clinical decision support systems that can exploit both the power of data-driven analytics and mechanistic models [20, 21]. In the future, we envision such hybrid models deployed at the bedside, providing the clinician with information regarding hidden cardiovascular states and optimal treatment strategies. For example, estimating volume and contractility states, and helping in guiding choices regarding intravenous fluid administration or inotropic support, a current unmet clinical need [22]. The model presented here is the first step towards this ambitious goal.
Materials and methods
Ethics statement
The study was approved with waiver of informed consent by the Research Ethics Board at The Hospital for Sick Children, approval number 1000048904.
Mechanistic model
The goal of this section is to detail iCVS: a mechanistic model of the cardio-vascular system that links measurable variables to clinically relevant physiological parameters that can not be measured directly. We use a lumped parameter model with several compartments: the circulatory system is represented as a closed loop with two compliant vessels (the arterial and venous compartments) and one pure resistance vessel (the organs) [5, 23]. The heart contains one chamber that acts as a pump. The venous-arterial flow IH is generated by the heart and is given by the product of the heart rate Hr and the stroke volume SV: (3)
The arteriolar and capillary IC flows are modeled as a linear resistor: (4) where Pa and Pv are the arterial and venous mean pressures, respectively, and R is the peripheral resistance. The evolution of arterial and venous volumes (Va and Vv respectively) is given by the following differential equations: (5) (6) where Iex is the intra-vascular volume change—either loss of volume due to bleeding or extra-vascular capillary leak or conversely intravenous fluid administration.
Eqs 5–6 relate between blood volumes dynamics and other hidden cardio-vascular variables. The challenge however is that blood volumes cannot be directly measured in patients. We therefore reformulate the mechanistic model in terms of blood pressures (which are measurable); this allows us to infer hidden cardio-vascular parameters based on observable data.
The pressure in the venous and arterial compartments is derived from the vessel’s volume and compliance as follows: (7) Where α ∈ {a, v} represents either the venous or arterial components, respectively, Cα the compliance, and the unstressed volume.
For the arterial compartment, we assume constant unstressed volume [5]. The venous unstressed volume is modulated by an autonomic activation Stot (see below): (8)
Assuming constant Ca [17], differentiating Eq 7 and rearranging, we obtain: (9)
Substitute Eqs 3, 4 and 9 in Eqs 5 and 6 we get: (10) (11) where we assume that the unstressed arterial volume is constant [5] and thus , and also use the approximation [17].
Autonomic Nervous System control.
The model contains elements of autonomic control which affect different physiological quantities: heart rate Hr, peripheral resistance R, unstressed venous volume and cardiac contractility K (see below). The autonomic control contains two components: the baro-reflex and an external control unit. The baro-reflex, marked by Sb, depends on the arterial blood pressure (see also [5]): (12) where kb describes the sensitivity of the baro-reflex and Pset is the set point of the baro-reflex. In this work kb = 0.1838 mmHg−1 [5].
The autonomic control also contains an additional component, marked by S which does not depend on the current arterial pressure. This component is not included in the work of [5]. The overall autonomic activation is a sigmoid function of the summation of the two components: (13)
Both Sb and S take values between 0 and 1. In addition, the shape of Stot is designed such that the total autonomic activation Stot is also ranged between 0 and 1, and the value of Stot is approximately the mean of Sb and S. The two different components of the autonomic activation represent two complementary functions: one is the baro-reflex which is part of a negative feedback loop aimed at maintaining blood pressure and cardiac output homeostasis while the other (S) represents the combined effect of sympathetic and parasympathetic activation on the cardiovascular system that is independent of current blood pressure. This division allows us to capture two different behaviours of the cardiovascular dynamics that are commonly observed in critically-ill patients: fluctuations which are negatively correlated between the arterial BP and heart rate, and synchronous fluctuations of the arterial pressure and other components such as heart rate.
Heart rate.
The heart rate is determined by the autonomic activation, and can vary between minimal and maximal values, Hrmin and Hrmax respectively: (14)
The minimal heart rate (Hrmin) and maximal heart rate (Hrmax) are parameters which are approximated by the optimization algorithm with a certain range that depends on the patient age (see Table 4).
Peripheral resistance.
The peripheral resistance varies between two values, Rmin and Rmax, and depends on the autonomic activation and on a local modulation, denoted MSVR, which does not depend on the autonomic nervous system: (15) where: (16) The value of MSVR can vary between −1 to 1.
Cardiac contractility.
The model contains one heart chamber. Of note, Zenker et al. [5], modeled the end systolic volume as a function of the end diastolic volume as follows: (17) Where Ves and Ved are the end systolic and end diastolic chamber volumes respectively, is the unstressed volume of the chamber, Ped is the end diastolic pressure of the chamber and K is the contractility. Eq 17 contains information about the heart volumes. Since heart volumes are not directly measurable, we modify the approach adopted by Zenker et al. and derive below an expression that relates heart contractility and blood pressures. We model the heart chamber as a compliance vessel, and estimate the difference between the end diastolic volume (Ved) and the unstressed ventricular volume as follows: (18) where Cven is the compliance of the ventricle. The intra-ventricular pressure (Ped) cannot be measured directly. Since the heart in the model contains only one chamber, we approximate the end diastolic pressure as the venous pressure: (19)
We acknowledge that this approximation is less accurate when the venous pressure is much larger than the pressure during diastole. This can happen for example in case of pulmonary hypertension, valve stenosis, or abnormal heart anatomy. In addition, the stroke volume (SV = Ved − Ves), can be approximated using the pulse pressure and the arterial compliance [17]: (20)
The contractility is modulated by the autonomic activation and is expressed as follows: (21) where Kmin and Kmax are the minimal and maximal contractility, respectively. Combining Eqs 17–21 we obtain: (22) The approximated expression for Ca ⋅ Pp(t) which is presented in Eq 22 can be substituted in IH which appears in Eq 10 to obtain a system of differential equations for Pa, Pv.
We have therefore introduced a model that relates between blood pressure dynamics and other cardiovascular parameters. We show below how this mechanistic model can be exploited in order to estimate hidden physiological parameters from routinely bedside collected data.
Simulations of the model
In Fig 3 and S1 Fig simulations of the iCVS model are presented. The hidden parameters are fixed and the observables are simulated according to Eqs 10–16 and 22. The following values for the hidden parameters are used: Pset = 50 mmHg, Kb = 0.1838mmHg−1, Hrmax = 180 Bpm, Hrmin = 100 Bpm, , , , , , , ΔVv0 = 65 ml. In Fig 3 S = 0.2. The initial conditions of the mean arterial and venous pressure are set to 48 mmHg and 4 mmHg respectively.
Extracting the observable variables from the raw data
The raw data which is used in this work is a continuous (sampled at 125 Hz) measurement of the venous and arterial pressure waveforms. The data analysis is performed as follows:
- Using the algorithm developed by [16, 17], we extract the following measurements from the raw data (see Table 2): (mean arterial pressure), (mean venous pressure) (heart rate), (pulse pressure) and , which is an approximation for the peripheral resistance multiplied by the arterial compliance and scaled by a constant (αRC): .
- Outliers are found and replaced—Outliers are defined as points which are smaller than the Lower threshold value—three local scaled median absolute deviation below the the local median within a sliding window of 20 sec, or higher than the Higher threshold value—three local scaled median absolute deviation above the the local median within a sliding window of 20 sec. The outliers are replaced by the following: Lower threshold value for elements smaller than the lower threshold value, and upper threshold value for elements larger than the upper threshold value.
- A linear interpolation is performed in order to obtain a data-set in 10 Hz frequency.
- A low pass filter is applied and the frequency components above Hz are removed.
- The data of each subject is divided into segments of 300 sec. Only Valid segments are used for further analysis. A valid segment should satisfy the following criteria: , , , , .
Sample subjects analysis
For each subject in Fig 9, the estimated parameters of 12 segments of 5 minutes are averaged. Each segment is 300 seconds long, and the segments are partly overlapped—each segment starts 100 seconds after the previous one. Segments which are not valid are not analyzed.
For the analysis in Fig 9, for each subject the first recorded 25 minutes time window that satisfies the follows is analyzed: the subject did not receive intra-venous fluids bolus during the time window, 15 minutes before and 15 minutes after.
In Figs 5, 6, 7 and 8, the above described 25 minutes time windows of patients 1,5,4,9 are presented, respectively.
The estimation procedure
Definition of the cost function.
Based on the iCVS model, we define a cost function which consists of a set of terms that link between the observables (Table 2) and the hidden parameters (Table 1). We first define the following expressions: , , (relative contractility), , , , and .
The cost function is defined as follows: (23) where: (24) (see Eq 14), (25) (see 15), (26) (see Eq 22), (27) (see Eq 11), (28) (see Eqs 10 and 15).
Rmodulation(t) is defined in Eq 16, and Sb(t), S(t), MSVR(t) are defined below: (29) (30) (31)
The parameters {NH, NR, NP, NI, Na} in Eq 23 normalize the different terms such that they have no units, and are adjusted in each time interval of the analysis: , ,, , and NI is the upper bound of . The parameters {wH, wR, wP, wI, wa} adjust the relative contribution of each term to the cost function. Throughout this work we use the following values: wH = 1, wR = 0.2, wP = 0.2, wI = 0.2, wa = 10. The values of {wH, wR, wP, wI, wa} are chosen such that the contributions of the terms , and to the cost function have the same order of magnitude, and the contributions of and (the terms that contain derivatives) are larger. The optimal parameters are searched within a certain physiological range (see Table 4).
Minimal and maximal values of Pset are the 10th and 90th percentile of the mean arterial pressure respectively, adapted from [11]. Bounds of minimal heart rate are the 5th percentile and 25th percentile (adapted from [10]). Bounds of maximal heart rate are the 75th percentile and 95th percentile (adapted from [10]). In addition, and are the values of minimal and maximal arterial compliance, respectively. See [24].
Estimation process.
Table 1 contains a list of the estimated parameters, where each parameter is bounded within a certain range (see Table 4). The values of the parameters Hrmin, Hrmax, Rmin, Rmax are updated according to the measurements. For example, if for a given patient a value of Hr = X is measured, for this specific patient we assume that the minimal heart rate cannot be larger than X and the maximal heart rate cannot be smaller than X in the next iterations. The estimation process is done separately in segments which are defined between two time points: [t, t + T], where T = 300 sec. During the estimation procedure, the hidden parameters are estimated by minimizing the time integral over the cost function f(t) within a given interval. For some of the parameters (see θ* below) we assume that the time dependent change is smooth, and consider also the previous estimated value of the parameters.
The optimization problem can be stated as: (32) where are the parameters for which previous estimators are considered, and Bi is the upper bound of the parameter θi (see Table 4). Throughout this work β = 0.1, τ = 500 sec.
The optimal solution (the set of optimal parameters ) in each segment is chosen using the function fmincon.m of Matlab with the option globalsearch, which searches a global solution for a constrained nonlinear optimization problem across different initial conditions [25]. We supply a random start point to the optimization function.
Supporting information
S1 Fig. Demonstration of the model for time-varying autonomic control.
S1 Fig demonstrates the effect of S (independent autonomic control) on the cardio-vascular dynamics in a simulation of the iCVS model. In response to an increase in S, arterial pressure, venous pressure, vascular resistance, and pulse pressure increase as well. The hidden parameters are fixed and the observables are simulated accordingly (see Methods). A. Time dependent intra vascular volume change—zero in this simulation. B. Time course of MSVR—constant in this simulation. C. Time dependent magnitude of the independent autonomic control (S) which does not depend on the cardio-vascular state. D-H The resulting observables: D. Arterial pressure (Pa), E. Venous pressure (Pv), F. Peripheral resistance multiplied by arterial compliance (RC). G. Heart rate (Hr), H. Pulse pressure (Pp). I. Time dependent magnitude of the baro-reflex (Sb) which depends on the arterial blood pressure.
https://doi.org/10.1371/journal.pcbi.1010835.s001
(TIFF)
S2 Fig. Reconstruction of the observables based on the optimal parameters which are found by the model.
A-C. Comparison between the measured value (horizontal axis) and estimated value (vertical axis). Each point represents the value of a given measurement (A—heart rate, B—peripheral resistance, C—pulse pressure) in a specific time point. The time points which are presented are a randomly chosen subset that constitute 1% of the whole data-set which is used in Fig 9. D—F. Comparison between observed data and reconstructed values during a segment of 500 seconds. Green—patient number 2, blue—patient number 3 and red—patient number 6 (see Table 3). Continuous line—observed data: (D), (E) and (F). Dashed line—reconstructed value: (D), (E) and (F). G-I. Histograms of the absolute values of the differences between the measured and estimated variables of heart rate (G), peripheral resistance (H) and pulse pressure (I). The dataset is same as presented in panels A-C. The arrows denote the mean differences for each of the patients presented in panels D-F using the same color scheme.
Using the estimated parameters which are found by the optimization process, we reconstruct the heart rate, peripheral resistance and pulse pressure, and compare the reconstructed values to the observed values. The reconstruction is based on Eqs 24–26, where we define: (33) (34) (35) where , and are the reconstructed values of the heart rate, vascular resistance and pulse pressure respectively.
S2A–S2C Fig show the reconstruction for a randomly chosen 1% of the whole data-set which is used in Fig 9. Each dot is a single time point of one patient (all time points of a given patient are represented by the same color). Horizontal axis represents the observed value (, , in panels A, B, C respectively), and vertical axis represents the reconstructed value (, , in panels A,B, C respectively). The Pearson correlation coefficient is 0.95 for Heart rate reconstruction, 0.98 for RC reconstruction, and 0.99 for Pulse pressure reconstruction. S2D–S2F Fig shows examples of reconstructed values vs. observed values. Panels D-F present the reconstruction of the observables for three patients, in a single interval of estimation. Each panel represents a given observable, and each patient is represented by a specific color across panels D-F. Continues line denotes for extracted observables, and dashed line for the reconstructed values. One can see that the iCVS model can reproduce the observables in different scales. The fast fluctuations are not reproduced, because MSVR and S are modeled as linear in each interval. S2G–S2I Fig shows histograms of the absolute value of the differences between measured and estimated values, for the random subset of the data which is shown in S2A–S2C Fig. The arrows show the mean differences in the interval which are presented above, each patient is represented by the same color as in S2D–S2F Fig.
https://doi.org/10.1371/journal.pcbi.1010835.s002
(TIFF)
S3 Fig. Robustness of the optimization process in a real dataset.
S3 Fig shows the results of the optimization process for different start points supplied to the optimization function. The illustration is done for patient number 1 which is identified in a hypovolemic shock state and is presented in Fig 5. Panels A-E present the extracted observables. Panels F-H present the optimal parameters obtained for realizations of the optimization process with different start points which are supplied to the optimization function. Grey lines denote for different start points, black lines are the results that appear in Fig 5. Purple lines are the average of all the realizations with different start points. It can be seen that the optimization algorithm which is used in this work to fit the iCVS model is affected by the supplied start points. However, when averaging across different initial conditions the results are concordant with the labeled shock state (hypovolemic in this case). In addition, in each single realization the average Iex is negative. We used Matlab version 2018b, running on 8 i7–8550U CPU @ 1.8GHz, 16.0GB RAM. The running time for running our optimization routine on a single interval with multiple initialization sets of parameters is 6 to 8 minutes. However, we believe that the running time can be shorter using a more optimal code, a stronger computer and maybe a different programming language. Results for patient number 1 (see Table 3) are presented. A-E.—The observables: Mean arterial pressure (A), mean venous pressure (B), heart rate (C.), RC—the peripheral resistance multiplied by the arterial compliance (D) and pulse pressure (E). F-H. Lines represent a smoothed version of the estimated parameters which are obtained in each interval (see Methods). Black line—the results which are presented in Fig 5. Purple line—average over all realizations. F. Relative intra vascular volume change (). G. Non autonomic vascular resistance (MSVR). H. Maximal relative contractility ().
Exploring identifiability: examining the ability of the estimation process and the optimization function to recover predefined simulated model patients and parameters. S4–S6 Figs show the results of the estimation process for the three sample artificial patients which appear in Fig 3 in the main text. An artificial patient with a hypovolemic shock state is presented in S4 Fig, an artificial patient with a distributive shock state is presented in S5 Fig, and an artificial patient with a combined hypovolemic and distributive shock state is presented in S6 Fig. In each figure, we depict the results of repeat runs of the estimation process with varying random initialization conditions and optimization function weights. While this is not a complete characterization of identifiability or estimation robustness, these results do provide some indication of our model’s and estimation process’s ability to correctly recover the hidden values and states. Examining these figures shows that while there is some variability in the estimations which stems from varying the initialization values of the estimator, the model still reproduces the hidden preset parameters used to generate the simulated patients fairly well. The results are more variable for the larger magnitude of Wa (Panels B,E,H in S4–S6 Figs), at least in these particular simulations. Also, it seems that the variability is reduced and accuracy is better when the hidden parameters to be estimated are time-varying. This probably stems from our choice to adaptively modulate the limits of “allowed” ranges for some parameters (see Methods) according to the patient’s specific history.
https://doi.org/10.1371/journal.pcbi.1010835.s003
(TIFF)
S4 Fig. Demonstration of the estimation process on simulation with predefined model parameters—An artificial patient with a hypovolemic shock state.
A-C.—The original parameters of the simulation chosen arbitrarily: relative intravascular volume change: (A), MSVR (B), and Kmax (C.). D-I. Each grey line represents a smoothed result of an estimation run obtained for each time interval (see Methods). The purple line is an average over all repeated runs (N = 5, different initialization conditions). We performed the estimation process using two different sets of weights for the optimization function to demonstrate the relative robustness of the estimation process to weight choices: Each column (D-F and G-I) represents the results that are obtained with a given set of weights as marked below. D,G Estimated relative intra vascular volume change (). E,H Estimated non autonomic vascular resistance (MSVR). F,I estimated maximal relative contractility ().
https://doi.org/10.1371/journal.pcbi.1010835.s004
(TIFF)
S5 Fig. Demonstration of the estimation process on simulation with predefined model parameters—An artificial patient with a distributive shock state.
A-C.—The original parameters of the simulation chosen arbitrarily: relative intravascular volume change: (A), MSVR (B), and Kmax (C.). D-I. Each grey line represents a smoothed result of an estimation run obtained for each time interval (see Methods). The purple line is an average over all repeated runs (N = 5, different initialization conditions). We performed the estimation process using two different sets of weights for the optimization function to demonstrate the relative robustness of the estimation process to weight choices: Each column (D-F and G-I) represents the results that are obtained with a given set of weights as marked below. D,G Estimated relative intra vascular volume change (). E,H Estimated non autonomic vascular resistance (MSVR). F,I estimated maximal relative contractility ().
https://doi.org/10.1371/journal.pcbi.1010835.s005
(TIFF)
S6 Fig. Demonstration of the estimation process on simulation with predefined model parameters—An artificial patient with a combined hypovolemic and distributive shock state.
A-C.—The original parameters of the simulation chosen arbitrarily: relative intravascular volume change: (A), MSVR (B), and Kmax (C.). D-I. Each grey line represents a smoothed result of an estimation run obtained for each time interval (see Methods). The purple line is an average over all repeated runs (N = 5, different initialization conditions). We performed the estimation process using two different sets of weights for the optimization function to demonstrate the relative robustness of the estimation process to weight choices: Each column (D-F and G-I) represents the results that are obtained with a given set of weights as marked below. D,G Estimated relative intra vascular volume change (). E,H Estimated non autonomic vascular resistance (MSVR). F,I estimated maximal relative contractility ().
https://doi.org/10.1371/journal.pcbi.1010835.s006
(TIFF)
Acknowledgments
We thank Ori Linial, Tomer Meir and Ron Teichner for useful discussions. We also thank Robert Greer, Tiffany Ramanathan and Will Dixon for helping with dataset consolidation.
References
- 1. Ellwein L, Pope S, Xie A, Batzel J, Kelley C, Olufsen M. Patient-specific modeling of cardiovascular and respiratory dynamics during hypercapnia. Mathematical biosciences. 2013;241(1):56–74. pmid:23046704
- 2. Olufsen MS, Ottesen JT. A practical approach to parameter estimation applied to model predicting heart rate regulation. Journal of mathematical biology. 2013;67(1):39–68. pmid:22588357
- 3. Fanelli A, Vonberg FW, LaRovere KL, Walsh BK, Smith ER, Robinson S, et al. Fully automated, real-time, calibration-free, continuous noninvasive estimation of intracranial pressure in children. Journal of Neurosurgery: Pediatrics. 2019;24(5):509–519. pmid:31443086
- 4. Chase JG, Preiser JC, Dickson JL, Pironet A, Chiew YS, Pretty CG, et al. Next-generation, personalised, model-based critical care medicine: a state-of-the art review of in silico virtual patient models, methods, and cohorts, and how to validation them. Biomedical engineering online. 2018;17(1):1–29.
- 5. Zenker S, Rubin J, Clermont G. From inverse problems in mathematical physiology to quantitative differential diagnoses. PLoS computational biology. 2007;3(11):e204. pmid:17997590
- 6. Denaï MA, Mahfouf M, Ross JJ. A hybrid hierarchical decision support system for cardiac surgical intensive care patients. Part I: Physiological modelling and decision support system design. Artificial intelligence in medicine. 2009;45(1):35–52. pmid:19112012
- 7. Li N, Cruz J, Chien CS, Sojoudi S, Recht B, Stone D, et al. Robust efficiency and actuator saturation explain healthy heart rate control and variability. Proceedings of the National Academy of Sciences. 2014;111(33):E3476–E3485. pmid:25092335
- 8. Stromberg D, Carvalho K, Marsden A, Mery CM, Immanuel C, Mizrahi M, et al. Standard CPR versus interposed abdominal compression CPR in shunted single ventricle patients: comparison using a lumped parameter mathematical model. Cardiology in the Young. 2022;32(7):1122–1128. pmid:34558399
- 9. Sproul A, Simpson E. Stroke volume and related hemodynamic data in normal children. Pediatrics. 1964;33(6):912–918. pmid:14169636
- 10. Eytan D, Goodwin AJ, Greer R, Guerguerian AM, Laussen PC. Heart rate and blood pressure centile curves and distributions by age of hospitalized critically ill children. Frontiers in pediatrics. 2017;5:52. pmid:28367430
- 11. Roberts JS, Yanay O, Barry D. Age-based percentiles of measured mean arterial pressure in pediatric patients in a hospital setting. Pediatric Critical Care Medicine. 2020;21(9):e759–e768. pmid:32740191
- 12. Sangkum L, Liu GL, Yu L, Yan H, Kaye AD, Liu H. Minimally invasive or noninvasive cardiac output measurement: an update. Journal of anesthesia. 2016;30:461–480. pmid:26961819
- 13. Thomsen KK, Kouz K, Saugel B. Pulse wave analysis: basic concepts and clinical application in intensive care medicine. Current Opinion in Critical Care. 2023;29(3):215–222. pmid:37078625
- 14. Greiwe G, Balfanz V, Hapfelmeier A, Zajonz TS, Müller M, Saugel B, et al. Pulse Wave Analysis Using the Pressure Recording Analytical Method to Measure Cardiac Output in Pediatric Cardiac Surgery Patients: A Method Comparison Study Using Transesophageal Doppler Echocardiography as Reference Method. Anesthesia & Analgesia. 2022;135(1):71–78. pmid:35452017
- 15. Goodwin AJ, Eytan D, Greer RW, Mazwi M, Thommandram A, Goodfellow SD, et al. A practical approach to storage and retrieval of high-frequency physiological signals. Physiological measurement. 2020;41(3):035008. pmid:32131060
- 16. Fazeli N, Hahn JO. Estimation of cardiac output and peripheral resistance using square-wave-approximated aortic flow signal. Frontiers in physiology. 2012;3:298. pmid:22934049
- 17. Erez E, Mazwi ML, Marquez AM, Moga MA, Eytan D. Hemodynamic Patterns Before Inhospital Cardiac Arrest in Critically Ill Children: An Exploratory Study. Critical care explorations. 2021;3(6). pmid:34151279
- 18. Artificial Intelligence and Machine Learning based Software as a Medical Device (SaMD) action plan. Food and Drug Administration. 2021; p. 2021–06.
- 19. Ursino M, Innocenti M. Modeling arterial hypotension during hemodialysis. Artificial organs. 1997;21(8):873–890. pmid:9247177
- 20.
Linial O, Ravid N, Eytan D, Shalit U. Generative ode modeling with known unknowns. In: Proceedings of the Conference on Health, Inference, and Learning; 2021. p. 79–94.
- 21.
Teichner R, Eytan D, Meir R. Enhancing Causal Estimation through Unlabeled Offline Data. In: 2022 7th International Conference on Frontiers of Signal Processing (ICFSP). IEEE; 2022. p. 147–158.
- 22. Messina A, Bakker J, Chew M, De Backer D, Hamzaoui O, Hernandez G, et al. Pathophysiology of fluid administration in critically ill patients. Intensive Care Medicine Experimental. 2022;10(1):1–15. pmid:36329266
- 23.
Keener J, Sneyd J. Mathematical physiology: II: Systems physiology. Springer; 2009.
- 24. Ursino M, Antonucci M, Belardinelli E. Role of active changes in venous capacity by the carotid baroreflex: analysis with a mathematical model. American Journal of Physiology-Heart and Circulatory Physiology. 1994;267(6):H2531–H2546. pmid:7810748
- 25. Ugray Z, Lasdon L, Plummer J, Glover F, Kelly J, Martí R. Scatter search and local NLP solvers: A multistart framework for global optimization. INFORMS Journal on computing. 2007;19(3):328–340.