A high-capacity coverless image steganography method based on double-level index and block matching

1.   Introduction

Hemorrhagic shock is a form of hypovolemic shock determined by rapid and large loss of intravascular blood volume. It may progressively lead to decreases in cardiac output and oxygen and nutrient delivery, decreased tissue perfusion, organ damage, and eventually death following brain or myocardial hypoperfusion ^[1].

In general, traumas and, more in particular, the hemorrhagic shock, are the leading cause of death in the population aged 1 to 45 years, and the third leading cause of death worldwide both for civilians and the military ^[2]. Given that hemorrhagic shock leads to death in a short time if not promptly treated, the interest in finding the most effective protocol for quickly restoring definitive homeostasis is clear ^[3].

Depending on the amount of blood lost, different clinical manifestations may be present and, for this reason, the Advanced Trauma Life Support (ATLS) doctrine defines 4 classes of hemorrhage ^[4], based on the amount of blood lost, with associated symptoms for each class ^[5].

Despite the fact that ATLS therapeutic guidelines are widely considered the "standard" for treatment of hemorrhagic shock, the optimal method of resuscitation has not yet been established with certainty ^[6], since there is no consensus within the medical community for identifying the best approach.

Moreover, individual victims, having different morphological and functional characteristics, react differently to external adverse factors and the severity, extent and type of trauma depend in fact on a multitude of factors ^[7]. This heterogeneity makes it difficult to reliably assess the overall state of the patient and define functional parameters following traumatic injuries accompanied by hemorrhage ^[7].

This difficulty stems in part from the fact that the clinician must simultaneously consider a host of factors and the impact each may have, individually or in combination, on the survival of the victim. Even if experimental data were available, the wide range of physiological variables to be considered and their complex dynamical relationship make it hard to have a precise and complete evaluation of the situation.

An effective way to study the complex cardiovascular system response to hemorrhage in humans is via mathematical modeling. With this approach, it is possible to integrate in a single (albeit complex) mathematical model the different functional sub-systems that are relevant to the response, focusing in particular to the faster-acting components of the circulatory system: it does becomes feasible to forecast the effect of combinations of external injuries and therapeutic maneuvers.

Probably, the most influential and complex mathematical model of the circulatory system that has been developed to date is the Guyton model, which has evolved over three stages ^[8]. Guyton first started with an extremely simplistic model of the entire cardiovascular system as a closed hydraulic loop ^[9]. It highlighted short-term regulations of cardiac output focusing mainly on the role of blood volume and vascular capacity ^[10].

In the following stage, Guyton sought to achieve long-term regulation of arterial pressure, for which the following slow-acting mechanisms were deemed crucial: (1) the capacity for wide fluctuations in urinary output with minimal changes in arterial pressure (the renal function curve in the Guyton model) and (2) long-term vascular autoregulation (via vasoconstriction or vasodilation) so as to adapt blood flow to tissue oxygen demands. At this second stage, the model afforded a better representation of the transient dynamics of renal hypertension as well as its steady state ^[11,12].

During the third stage of its evolution Guyton's model was further enriched with a variety of additional endocrine and neural mechanisms, which exert control on renal function ^[13]. Due to its complexity, the main application of Guyton's model is currently as a conceptual construct, in particular for educational purposes: recent and continuing advances in physiological understanding as well as in computing and simulation software have not yet been incorporated into a complete human circulation model ^[14].

Despite being the most widely known model of long-term blood pressure control (with several software implementations having been disseminated), no published report of any of these implementations has provided the details of all equations, parameter values, data used for model identification, or the model identification process and results. Indeed, these implementations, even those of the relatively simple 1972 version, had avoided formal identification until Kofránek et al. ^[15,16,17]. These latter authors not only formalized Guyton's model but also debugged it. These facts attest to the largely impenetrable and generally irreproducible nature of Guyton's seminal work. According to Beard, the Guyton-Coleman model focuses on the axiom that arterial blood pressure is, regardless of the circumstances, regulated primarily by the kidney ^[18]. As such, the inherent flaw in the Guyton-Coleman model is that it was not designed to consider alternative hypotheses ^[18].

It is clear that a model as complex as Guyton's was (and still is) difficult to use in a practical clinical setting, especially when the investigator might wish to estimate model parameters from experimentally available observations. In these situations, the use of a relatively more small and compact model, with few free parameters, yet still able to quantitatively represent the time course of the main variables of interest, is preferable. This need gave rise, in the eighties and nineties, to the development of a few alternative models. In the last decade, several models have been proposed to take into account the relevant phenomena of the heart and its interactions with the cardiovascular system. Major interdisciplinary research efforts have strived to better understand and model the pathophysiology of the cardiocirculatory system in multi-dimensional (0D, 1D, 3D) ^[19,20] and multi-scale (e.g., molecular, cellular, organ, systemic) ^[21,22] representations.

Among the more recent 0D and compact models, Zenker's model would seem particularly useful in the attempt to represent the chain of compensation mechanisms following hemorrhage: Zenker's mathematical model uses a simplified representation of the cardiovascular system with its inherent autoregulatory control, also incorporating principles of Bayesian inference, to simulate volume loss (exemplified by acute hemorrhage) ^[23].

Notwithstanding its simplicity, it offers the ability to portray the time evolution of most variables of interest to healthcare professionals caring for a wounded victim (e.g., mean arterial pressure, heart rate and cardiac output).

It would therefore seem that Zenker's model ^[23], among all currently available models, offers the best compromise between simplicity of the formulation and relevance of the forecasts for the clinician.

The goal of the present work is to compare the clinical predictions obtained via the extensive Guyton model ^[13] to the compact Zenker model ^[23] following simulated hemorrhage, so as to highlight any significant differences concerning key physiological variables and assess their overall practical application in hemorrhagic shock ^[14]. The choice fell on these two models, among the vast panorama of existing cardiovascular mathematical models, because they seem to be at opposite extremes as regards the complexity and the detail of information that they analyzed.

The performances of the two models are evaluated considering the qualitative response to bleeding at different durations of hemorrhage (ranging from a few minutes to a few hours) and at different volumes of blood loss (from minor to fatal). The strengths and weaknesses of the two models are discussed, with the objective of assessing their potential validity in reliably predicting the evolution of the physiologic conditions of trauma victims.

The paper is organized as follows. We begin in Section 2 with a detailed description of the two investigated models in which we report the limits and the potential of each models. To analyze the response of cardiocirculatory models to hemorrhage, a series of simulations were performed at different degrees of blood loss and for different durations of bleeding. Simulations details and results are presented in Section 3 whilst the discussion is in Section 4. Finally, conclusions and comments are summarized in Section 5.

2.   Materials and method

Given the level of complexity that often characterizes natural phenomena, the development of precise mathematical models usually proves quite challenging. The key concept is that of a system, which is a whole whose overall behavior and functioning depend on the interactions of the constituent parts rather than on the sum of the actions of each part taken separately ^[24]. Therefore, the process of the physics modeling obtained by applying idealizations, simplifying hypotheses and approximations is an inevitable but necessary practice for mathematical modeling of physical systems to ensure some degree of reproducibility and comparability ^[25]. The level of complexity of a mathematical model is directly proportional to its accuracy in representing natural phenomena ^[25]. Starting from these considerations, in the panorama of the mathematical models that describe cardiac circulation, we have decided to take into consideration the two that we believe are at the extremes as regards the detail of the information that they analyze: the Zenker model and the Guyton model. The following offers a brief description of the mathematical models considered, also including the relative historical contexts.

2.1.   Guyton model

The Guyton, Coleman and Granger model ^[13] is probably the most widely known and discussed model in the history of cardiovascular research. Arthur Guyton's research covered the topics of short- and long-term regulation of arterial pressure and body fluid circulation. His seminal research still represents a mainstay of cardiocirculatory physiology and related disorders.

Terms such as cardiac output and venous return entered the medical lexicon thanks to his work spurring research into the role of the renal system in local regulation of blood flow and whole-body circulatory autoregulation (i.e., pressure-natriuresis mechanism). His greatest intuition was to apply engineering control system analysis in order to model the cardiovascular-renal system. He provided one of the first comprehensive cardiovascular models implementable on a computer and his mathematical framework remains as a reference for generations of cardiovascular researchers.

Trying to synthesize Guyton's accomplishments in the field, spanning over three decades, it is possible to identify three different generations of mathematical models ^[8]:

● in his first works, Guyton developed a simple model of the cardiovascular system focusing on the role of blood volume and vascular capacity by obtaining a closed hydraulic loop which accounted for short-term regulations of stroke volume (i.e., cardiac output) ^[10,26];

● in a second phase, Guyton focused his attention on studying long-term regulation of blood pressure and, in particular, he devoted himself to two fundamental slow-acting mechanisms: (1) the curve of renal function and (2) long-term vascular autoregulation. The first mechanism is able to determine large fluctuations in urinary output with minimal variations in arterial pressure, while the second defines the extent of vascularization through constriction or dilation of vessels ^[11,27,28].

● the third version of Guyton's model ^[13] (hereinafter referred to as G-72) includes endocrine and neural mechanisms that complete the control action of renal function ^[8,12].

Although the model has a purely hierarchical organizational structure, it has always been historically seen as a simple set of differential equations. This latest version of the model consists of 18 modules (comprising over 350 elementary components), including approximately 160 variables, 36 of which represent variables of state and nearly 500 model variables, parameters and constants ^[17].

Because the G-72 model is an extensive mathematical model of human circulatory physiology, extended over time to include many complex control processes, it allows for a wide range of time scales in the individual physiological subsystem, ranging from $5\cdot 10^{{-}4}$ min (autonomic control) to $10^{4}$ min (heart hypertrophy or deterioration) ^[29]. Specifically, the renal control of blood pressure obtained by means of long-term regulatory effects on cardiovascular activity (such as the renin-angiotensin-aldosterone system) adopt time constants measured in hours or days, while short-term regulation (mainly baroreceptor reflexes) uses time constants in the order of seconds.

The G-72 model is able to perform several simulations mimicking effects resulting from different types of stresses on the cardiovascular system. In fact the possible physiological behaviors predict by these simulations are obtained in a shorter time than the experiments conducted in vivo on animals. For more details and for a working version in Simulink of the G-72 model readers are referred to http://patf-biokyb.lf1.cuni.cz/wiki/workshopy/guyton.

Furthermore, since the system cannot take into account all cardiovascular pathologies, it has served to highlight its potential and limitations, thus helping in planning further research.

This model was further refined, making it even more complex and improved, by different authors. Below we mention only some of the most significant works. In 2010 Kofránek proposes a model development based on Guyton's original diagram for research into the causes of cardiovascular pathogenesis and for educational purposes ^[17]. Ikeda ^[30] included acid-base regulation parameters and variables (i.e., inputs/outputs), contemplating a greatly expanded list of solutes and metabolites. Verily, the respiratory subsystem, which is functionally interrelated with the acid-base balance subsystem, was derived from the detailed model provided by Grodins et al. ^[31,32].

Nowadays many researchers are still studying the Guyton model, always finding in its ideation of the cardiovascular system new hints and reasons for open discussion for deeper research on the subject ^{[33,34,35,36]}.

2.2.   Zenker model

In 2007, Zenker et al. proposed a lumped parameter model to describe the cardiovascular system, in which {the} heart (left ventricle) is modeled like a one-way pump that pushes the blood along the network of blood vessels ^[23]. The authors, in order to predict the hemodynamics throughout the major arterial and venous systems, relied on the electrical analogy with the principles of fluid mechanics and represented the elastic blood vessels by electrical capacitors (according to the Windkessel model), while the arterial pressure was kept under control through baroreflex feedback.

The resulting model is composed by five ordinary differential equations (ODEs), which are reported below:

The first two equations represent the temporal trend of the end-systolic volume, $V_{ES}$ (2.1), and the end-diastolic volume, $V_{ED}$ (2.2); the third and fourth equations describe the evolution of arterial volume, $V_a$ (2.3), as well as the venous volume, $V_v$ (2.4) and the last equation, the temporal evolution of the output, $S$ (2.5), comes from baroreflex central processing.

Note that equations (2.3) and (2.4) have been written assuming conservation of volume at the nodes, which represents the electric analogue of Kirchhoff's first law, while the term that concerns any blood losses (e.g., due to hemorrhages) or blood increments (e.g., from transfusions) is $I_{external}(t)$ in (2.4), i.e., we always have a venous bleeding.

The baroreflex control system has the task of maintaining cardiovascular homeostasis. This crucial regulatory mechanism is implemented through a logistic function and a linear system in a cascade configuration ^[37,38]. A blood pressure value, $P_{a_{set}}$, fixed a priori, establishes the equilibrium point of the baroreflex feedback loop. A first order low-pass filter represents the linear function of the central baroreflex control system and, moreover, the reaction times depend on the unstressed venous volume control mechanisms.

$\tilde{V}_{ES}$ is the current estimate of end-systolic ventricular volume ^[23]:

where $k_{VED0}$ is the empirical value of volume at which the passive intraventricular pressure is 0 mmHg ^[39] and $P_{ED}$ is the intraventricular pressure at the end of diastole, experimentally characterized by the following exponential relationship ^[39]:

where $k_{ELV}$ and $P_{0LV}$ are constants describing the ventricular pressure/volume relationship estimated from experimental data for the left ventricle ^[39].

The expression of $\hat{V}_{ES}$ is ^[23]:

where $c_{PRSW}$ is the preload recruitable stroke work (a contractility index representing the relationship between stroke work and diastolic filling) ^[39].

$\tilde{V}_{ED}$ is the current estimate of end-diastolic ventricular volume ^[23]:

where $\hat{V}_{ED}$ is given by:

where $T_{Sys}$ is the duration of systole and the function $V(t)$ is given by Equation ^[39]:

where

$R_{valve}$ is a valve opposing ventricular filling, $P_{CVP}$ is the central venous pressure and $P_{ES}$ is the intraventricular pressure at the end of systole:

From the above ODE system it is possible to derive the trend of other hemodynamic variables:

The use of this mechanistic model allows us to observe, with an acceptable degree of approximation, the time evolutions of the pressure and flow within the arterial and venous compartments, minimizing the number of physiological variables and parameters. However, the study of the dynamics that regulate the complex interactions of fluids within the human body still represents an extremely challenging issue.

In the Zenker model the total volume of circulating blood, denoted $V_{total}(t)$, is given by the sum of arterial and venous volumes, which are in a ratio of approximately $1$ to $4$, respectively, ($V_{total}(t) = V_a(t) + V_v(t)$). Furthermore, both venous and arterial blood are in turn subdivided into a stressed component and an unstressed one: $V_{v}$, $V_{v_0}$ $V_{a}$ and $V_{a_0}$, respectively ^[23].

The system devised by Zenker offers the possibility of analyzing the temporal trends of four important variables, namely $P_a$ (arterial pressure), $P_v$ (central venous pressure), $f_ {HR}$ (heart rate) and $I_ {CO}$ (cardiac output), during severe hemorrhage and/or in response to a possibly life-saving transfusion.

The first prerequisite of any modeling approach consists in obtaining reliable measurements that more directly mirror key determinants of clinically significant outcomes than traditional physiologic variables. Zenker's model was subjected to such measurements, subsequently validated experimentally via laboratory (swine) animal models under conditions of controlled hemorrhage. The above process is referred to as the calibration of model parameters ^[40] and all the parameters of the Zenker model and their sources are available in ^[23]. The ultimate goal of cardiocirculatory mathematical models is to effectively predict the responses to variations of key physiological parameters, whether in experimental animal protocols or in simulated clinical trials, for the virtual assessment of therapeutic strategies.

3.   Results

In order to analyze the behavior of the studied models and their response to a hemorrhagic event, first of all, it was necessary to define a simulation setup which is reported in the first part of this section, while, at the end of this section, the results of the simulations are presented for both the models taken into consideration.

3.1.   Simulation setup

The defined simulation setup can be schematized in the following operations:

1. identification of the variables of interest;

2. definition of hemorrhage classes;

3. definition of durations of hemorrhage events.

3.1.1.   Variables of interest

During a hemorrhage, loss of blood may result in insufficient ventricular filling during diastole such that the heart is no longer able to provide optimal blood flow to cells and tissues. In response to acute blood loss, neurohumoral activation at effector organs provides immediate countermeasures to maintain blood pressure and regional flow to vital districts (i.e., cerebral, hepatic and cardiac), thus assuring gas exchanges to the most oxygen-dependent organs ^[41]. However, when these compensatory mechanisms are surpassed severe cellular hypoxia and organ damage may occur, without early and effective therapy, ultimately leading to death ^[1,41]. In the setting of hemorrhagic shock, the major challenge consists in designing an adequate model to accurately mimic the clinical situation, while maximizing reproducibility and standardization. In essence, three distinct approaches can be used in animal models of hemorrhagic shock:

● fixed-volume hemorrhage, whereby a certain percentage of the total calculated blood volume is drawn over a predetermined time interval ^[41];

● fixed-pressure hemorrhage, whereby bleeding is protracted until a desired drop in arterial pressure is achieved, then maintained constant via controlled hemorrhage or fluid infusion, as needed ^[41];

● uncontrolled hemorrhage, induced by inflicting an injury (e.g., a standardized crush/organ laceration, or amputation) ^[41].

Fixed-volume hemorrhage models offer the opportunity to study hemodynamic responses and compensatory mechanisms following fixed bleeding volumes, therefore, we opted for such a model to evaluate the performance of the Zenker and Guyton models in hemorrhagic shock.

According to the literature ^[41], the most common hemodynamic variables that are usually measured in a animal fixed-volume hemorrhage shock models are:

● Mean arterial pressure;

● Heart rate;

● Stroke volume;

● Cardiac output;

● Blood/plasma volume;

● Shock index;

● Cardiac index;

● Pulmonary capillary wedge pressure.

In particular in our simulations, we consider only the first five variables because they were defined in both models analyzed (i.e., Guyton and Zenker models), whereas the last three are not present in the Zenker model. Although the stroke volume is also not defined in the Guyton model, it is possible to calculate it indirectly from the cardiac output and the heart rate values. In Table 1 the acronyms of the variables in the models are reported.

As shown in the last column of the Table 1, a univocal notation for the variables of interest of both models is used, with the aim of facilitating the interpretation of the results. Instead, the last row of Table 1 shows that, for both models, there is no variable that represents the total blood volume in the cardiac system, although it can be calculated as the sum of other two variables. In particular, in the Guyton model, the total blood volume ($VTot$) is given as the sum of the plasma volume ($VP$) and the total red cell volume ($VRC$) and, for the Zenker model, the total blood volume is the arterial volume ($V_a$) plus the venous volume ($V_v$). In the Guyton model, cardiac output is defined as the volume of blood being pumped by the left ventricle per unit time, instead Zenker defined it as the heart rate times the difference between the end-diastolic and end-systolic volumes (2.18).

3.1.2.   Class and duration of hemorrhages

After that the variables of interest were identified for the considered mathematical models, we proceeded to define and identify the type of bleeding to be performed.

According to Advanced Trauma Life Support ^[4], it is possible to define four classes of hemorrhage based on the physiological effects of bleeding by observing clinical signs. This kind of classification is useful for estimating the percentage of volume of acute blood loss, as reported below:

● Class Ⅰ Hemorrhage: < 15% Blood Volume Loss;

● Class Ⅱ Hemorrhage: 15-30% Blood Volume Loss;

● Class Ⅲ Hemorrhage: 31-40% Blood Volume Loss;

● Class Ⅳ Hemorrhage: > 40% Blood Volume Loss.

Using this information, we performed the simulations considering the four classes of hemorrhage and we compared the two models under identical durations of bleeding times for the 4 standard bleeding classes. For class Ⅰ the bleeding times are: 5, 10, 15 and 30 min since the amount of blood drawn is only 10% of the total circulating volume and longer durations would be less realistic as a possible cause of hemorrhagic shock. For class Ⅱ we have bleeding times of 10, 15, 30 and 60 min inasmuch as the amount of blood loss is increased (to 20%). The third class (with 35% of blood loss) includes even longer bleeding times: 15, 30, 60 and 120 min. When there is a blood loss of > 40% it is considered a class Ⅳ hemorrhage. We chose a blood loss equal to half of the total volume, with times ranging between 30 min and 180 min, also including intermediate times of 1 and 2 h.

Note that, the percentage of blood volume loss is divided in a ratio of 6:4 of the $VP$ and $VRC$ in the Guyton model and in a ratio of 1:0 of $V_v$ and $V_a$ in the Zenker model, i.e., the volume loss is totally venous.

3.2.   Guyton simulations

Tables 2-5 comprise the five important variables taken into account (i.e., MAP, HR, SV, CO and VTot) determined for the 4 classes of hemorrhage. The values of the variables are calculated for each class of bleeding at set times: at the beginning of the bleeding ("At BB"), at the end of the bleeding ("At BE") and at 1, 7, 15 and 30 days after BE. It is important to note that the value of VTot variable, in the column indicating the end of the bleeding, does not equate to the difference between the initial VTot volume ("At BB") and the volume of hemorrhaged blood: this can be explained by the fact that Guyton's cardiovascular system presents a mechanism, with a slow but nevertheless effective action, which tends to restore blood volume (with more plasma) after a hemorrhage. This physiological compensatory mechanism, called transcapillary refill, commences at the onset of the bleeding and continues until homeostasis is reached.

The simulations were performed using MATLAB (MathWorks) programming environment. In particular, the ode23tb function was used to solve stiff differential equations using the trapezoidal rule plus backward differentiation formula. Max step size and min step size were set to auto, while the relative tolerance was set to 1e-3.

Figures 1, 3 and 5 show the temporal courses of the three important hemodynamic variables, i.e., MAP, HR and CO, for Class Ⅰ hemorrhages, while Figures 2, 4 and 6 for Class Ⅳ hemorrhages. All six figures refer to the case of a 30-min hemorrhage considering a time span of 30 days. The paired figures plot the temporal trends for the two extreme hemorrhagic classes (the least and most severe, Class Ⅰ and Class Ⅳ, respectively) for each hemodynamic variable considered, over a time interval of one month.

The diagrams clearly show that in Guyton's mathematical model all the variables considered return, over time, to their initial conditions prior to the traumatic bleeding event. In fact, in Guyton's model, as in the human body, the mechanisms of cardiovascular control maintain arterial and central venous pressures as constant as possible, seemingly by design. No doubt, any drop in arterial or central venous pressure due to, for example, a loss of blood volume would tend to modify the pressure relationships at the capillary level so as to facilitate a net fluid transfer into the intravascular compartment of the cardiovascular system. Figures 7, 8 and 9 show the temporal trends of MAP, HR and CO, for the 4 hemorrhagic classes in the case of a 30 min hemorrhage. The initial conditions for the model variables were gleaned from the literature, especially considering the works of Guyton. With these values, the model is not in a steady-state condition at the onset, but requires a few minutes for the system to reach a state of equilibrium, as can be seen from the "ripple" present in the first minutes of the Figures 7-9.

In Figure 7, MAP is simulated for subjects who undergo a bleeding of 10% (Class Ⅰ), 20% (Class Ⅱ), 35% (Class Ⅲ) and 50% (Class Ⅳ) of their total intravascular volume, respectively. From the onset of the hemorrhagic event, all of the MAP curves abruptly fall. Unsurprisingly, the greater the blood loss, the steeper the drop (Figure 7). For classes Ⅲ and Ⅳ the minimum point is reached considerably before the bleeding ends. From that point on, the curves invert their downward trend, (due to the sum of the complex compensatory mechanisms integrated into the Guyton model). Once again, the higher the bleeding class, the earlier and the steeper the rise. For the other two curves (classes Ⅰ and Ⅱ) both the fall and the rise are less pronounced, more delayed and much more gradual. In particular, with a loss of 10% of the circulating blood, the MAP curve drops steadily almost until the end of the bleeding and then rises, with an inverted slope, immediately afterwards.

The variations in heart rate (HR) are illustrated in Figure 8, plotted by curves for the 4 different classes of blood loss. Throughout the 30-min hemorrhage, the curves ascend with increasingly steeper slopes, the higher the class.

Figure 9 shows the variations in cardiac output (CO). In the case of relatively minor bleeding (10%), we can observe the adequacy of the compensatory mechanisms for most of the bleed, until they start to falter several minutes before it ends. Although less than in Class Ⅰ, the homeostatic mechanisms compensate for hemorrhage even in Class Ⅱ, yielding only around halfway through the hemorrhagic interval. For more severe bleeding (35% and 50%) we observe progressively decreasing values of CO, up to about the 15 min mark from the onset the bleeding, followed by an inversion in the slope of the curves. In all 4 classes, once bleeding stops there is a marked recovery in CO, due to the subsequently unresisted compensatory mechanisms.

3.3.   Zenker simulations

Tables 6 to 9 represent the five important variables: MAP, HR, SV, CO and VTot determined for the 4 classes of hemorrhage. The values of the variables are calculated for each class, at set times: at the beginning of the bleeding ("At BB"), at the end of the bleeding ("At BE") and at 20 minutes to BE, respectively.

Note that, in contrast to the corresponding values for the Guyton model, in the At BE column (i.e., volume when bleeding is interrupted) the VTot variable equals the initial VTot volume ("At BB") minus the percentage of blood lost: this is due to the fact that in the Zenker model of the cardiovascular system the physiological mechanism of the transcapillary refill is unaccounted for.

The numerical solution of the system Equation (Eqs (2.1)-(2.5)), as well as all other algorithms used for this work, was implemented in the MATLAB (MathWorks) programming environment, using the ode15s solver for numerical integration. Ode15s was used to solve stiff differential equations and Differential-Algebraic Equations (DAEs) with the variable order method. Max step size and min step size were set to auto, while the relative tolerance was set to 1e-3.

Figures 10 to 12 show the temporal trends of the 3 important hemodynamic variables, i.e., MAP, HR and CO, for the 4 hemorrhagic classes, given a 30-min hemorrhage.

In Figure 10, MAP is simulated for subjects undergoing a blood loss of 10% (Class Ⅰ), 20% (Class Ⅱ), 35% (Class Ⅲ) and 50% (Class Ⅳ) of their total intravascular volumes, respectively. Throughout the hemorrhage and after a brief initial overshoot, the top three curves for MAP show a gentle yet increasingly downward slope, ranging from barely perceptible for Class Ⅰ to progressively more marked going from Class Ⅱ to Class Ⅲ.

Similarly, the temporal evolution of the fourth MAP curve, which represents a loss of 50% of the total blood volume, descends with a slightly downward slope (albeit lower than the first 3) until about the 20-min mark. Then it drops precipitously until the end of bleeding. Once the bleeding stops, from that moment onwards all 4 curves level off at their respective MAP values. This behavior is due to the absence of any cardiovascular compensatory mechanism.

The trend in heart rate (HR) is described in Figure 11, which illustrates the curves for the 4 different classes of blood loss. During the bleeding and throughout all its duration, the curves all climb. The greater the blood loss, the progressively steeper the rise, particularly marked for Class Ⅳ (50% volume loss): by the end of the latter, HR more than triples, which is far greater than the other curves. As soon as the bleeding ceases, all the curves plateau at the level reached. Likewise to the previous variable, there is also an overshoot at the beginning of the simulations for each curve.

As to the last variable (i.e., CO), it remains almost constant, indicating that the compensation mechanisms are adequate, up to roughly a third of the way through the bleeding, but then incurs a decrease. Although slightly less efficiently than in the first curve, the homeostatic mechanisms hold fairly well even for Class Ⅱ. In both classes Ⅱ and Ⅳ hemorrhages the CO decreases appreciably about midway through the hemorrhagic interval. For the more severe blood losses (i.e., 35% and 50%) we see a decrease in CO already after the 10-min mark from beginning of the bleeding, followed by a progressively downward slope, especially for Class Ⅳ. For all 4 different cases, there is an immediate "plateau" as soon as bleeding stops, due to the lack of the trancapillary refill mechanism mentioned above (Figure 12).

With reference to certain simulation results in the literature ^{[42,43,44,45,46,47,48,49,50]}, it seems interesting to compare the trends obtained in our simulations of the Guyton and Zenker models, with those obtained by applying these models to particular experiments conducted on lab animals. In Guyton and Crowell's experimental results on dogs ^[42,43], for example, it is interesting to analyze the diagram plotting the controlled bleeding experiments, of various degrees of acute hemorrhage, and the corresponding effects on the time course of mean arterial pressure. By contrasting Guyton's plots and their respective simulation values (Figure 7) with Zenker's (Figure 10), we can clearly see how the former simulations outperform the latter.

Considering other experimental data from the literature ^{[44,45,46,47,48,49,50]}, focusing on the other two hemodynamic variables (i.e., HR and CO), we observe that again the Guyton model simulations are more accurate in predicting the time course of these two variables throughout hemorrhagic events.

4.   Discussion

By comparing the simulation diagrams obtained for the two models considered in our study, we can see that their temporal courses present different dynamics. The diagrams of the Guyton model are more complex and, although not designed for a specific study on hemorrhagic shock, its responses are able to predict the temporal course of the main hemodynamic variables both in the short- (Figures 7-9) and long-term (Figures 1-6). In particular, since Guyton's model encompasses various compensatory mechanisms, two time scales were required. Time-wise, these mechanisms can be subdivided into three categories: rapid-onset, within seconds or minutes; intermediate-acting, that is over minutes to hours; and long-term arterial pressure regulation within days, months, or longer ^[43]. The first category includes the following mechanisms: baroreflex; ischemia-induced, via central nervous system; chemoreceptors. The second comprises the renin-angiotensin vasoconstrictor mechanism, vascular tone (i.e., degree of smooth-muscle contraction), and transcapillary refill. The last consists in the renal mechanisms for controlling blood volume and pressure over the long term ^[43].

On the other hand, the Zenker diagrams have a more linear trend and their predictions are limited to the short-term description of the hemodynamic variables (Figures 10-12), given its relative simplicity, essentially limited to the baroreflex-dependent compensatory mechanism.

In general, however, the two models represent very similar behaviors (trends) during the bleeding phase. This behavior is due to the fact that the Guyton model is much more complex and integrates, in addition to the baroreflex system, another important mechanism designed to maintain a constant circulating fluid volume within the cardiovascular system : transcapillary refill.

This mechanism involves a cascade of inter-related processes that altogether underlie this complex function. A key feature is the mobilization and redistribution of fluid reserves within the extra-vascular compartment. This entails that fluid readily available within the body can be partly tapped to restore volume of the intravascular compartment. At the kidney, volume depletion evokes renal conservation of salt and water to limit further losses, eventually bolstered by the thirst response, leading to additional fluid intake ^[51]. As long as the overall systemic effects of the variations in volume, pressure and blood flow can be simulated as the system attempts to maintain homeostasis its real-time output offers valuable data and information ^[51].

The principal factors regulating the exchange of fluid across the capillary walls were proposed by Starling ^[52]. The transcapillary refill flow has direction and rate that depend on the relationship between the gradient of capillary hydrostatic and colloid osmotic pressure ^[51,53,54]. In our discussion we suppose that the hemorrhage is "simple" or that it occurs with a negligible amount of tissue damage. Under these conditions, as the bleeding proceeds, we witness a biphasic response of the human body to simple hemorrhage ^[55]: a first phase is characterized by hemorrhage of up to 10-15% of total blood volume (Vtotal) ^[56] in which the arterial baroreceptor reflex maintains a constant mean arterial pressure (MAP) (via tachycardia and an increase in peripheral vascular resistance) and a second phase that usually occurs when blood loss exceeds 20% of the Vtotal (whose effect may be protective) in which there is a decrease in the MAP by a second reflex ^[56].

Contrary to the Guyton model, incorporating a sophisticated reflex control of the circulation ^[13] and a long-term control of blood pressure through the baroreflex arc and the renin-angiotensin system ^[13], in the Zenker model the baroreflex mechanism is modified through the fifth ODE of the system:

as is well-known, by increasing total peripheral resistance (via a positive inotropic effect) and by means of tachycardia, the baroreceptor reflex limits the drop in MAP ^[57]. Therefore, in the graph (Figure 10) and Tables 6-9 that represent the MAP as a function of time we observe a slight increase in its value at the end of the bleeding, but, during hemorrhage, venous return decreases and the tachycardia may not be sufficient in supporting cardiac output ^[57] and hence MAP.

Since the Zenker model does not include any transcapillary refill mechanism (which would involve an increase in MAP), nor the more complex 'depressive' reflex that describes the second phase of hemorrhage, we can conclude that, post-bleeding, the MAP rises minimally, remaining constant thereafter, regardless of the eventual, impending the death of the individual, even when MAP approaches zero.

In our simulations the Guyton model reproduces the phenomenon of hemorrhagic shock more precisely than the Zenker model, thus more closely mirroring reality, in accord with the experimental data from the literature ^{[42,43,44,45,46,47,48,49,50]}. However, Guyton's model is undeniably more complex than its "rival" due to the far greater number of requisite parameters. While added complexity makes for a more realist model, it can render the model more cumbersome to understand and analyze, often with overburdening computational demands.

In addition, as with all more complex models, the Guyton model includes a number of somewhat poorly known parameters (i.e., not all parameters of the models can be identified) ^[58]. A second pitfall consists in the dearth of available data to estimate model parameters, associated with the inherent uncertainties of empirical measurements.

Last but not least, nonlinear systems often admit a host of potential solutions. Accordingly, the number of parametrizations theoretically compatible with observed data can be huge ^[58]. This inherent ambiguity has been dubbed "ill-posedness", defined as the absence of a unique solution to the inverse problem ^[23].

In fact, a major challenge faced by modelers is to design sound models that minimize the impact of such ambiguity. In practice, the ultimate goal is striking an acceptable trade-off, cognizant of the inverse relationship between model complexity-accuracy and simplicity-user friendliness. In this context, the Zenker model has significantly fewer parameters, but an expectedly lower predictive accuracy for real-world scenarios.

5.   Conclusion

We have presented two very different mathematical models of the cardiovascular systems in order to evaluate their strengths and weaknesses. Guyton's model, the most well-known and studied, represents one of the most complete and complex mathematical representations of the circulatory system and is still a valid reference point, despite many improvements and corrections ^{[15,16,17,59]}. However, during the hemorrhagic event, the Zenker model also adequately captures the changes in the circulatory system, albeit not with the same degree of precision and closeness to physiological reality.

Depending on the requirements of the specific scenario considered and / or the availability of a powerful calculator, one may choose either of the two models. We speculate that the Zenker model, opportunely enhanced and integrated with the transcapillary refill mechanism, could describe the hemorrhagic shock related to traumatic events just as well, if not better. In any case, Zenker's model would nonetheless maintain the advantage of being leaner than the more cumbersome, but more realistic, Guyton model.

The following phase of our inquiry will be to design a relatively simple model, such as Zenker's, and enhance it by incorporating the mechanism of the transcapillary refill. Once operational, it must be subsequently validated by comparing the simulations obtained with experimental data obtained in vivo.

Acknowledgments

The present work was carried out within the MoSES project (Modeling Shock in the Experimental Setting), an Italy-Israel bilateral technological cooperation between The Italian Ministry of Defence, Secretariat General of Defence and National Armaments Directorate, Contract n. 1888/2016 of Land Armaments Directorate.

Conflict of interest

The authors declare no conflicts of interest in this paper.