Abstract
One of the clinical examinations performed to evaluate the autonomic nervous system (ANS) activity is the tilt test, which consists in studying the cardiovascular response to the change of a patient's position from a supine to a head-up position. The analysis of heart rate variability signals during tilt tests has been shown to be useful for risk stratification and diagnosis on different pathologies. However, the interpretation of such signals is a difficult task. The application of physiological models to assist the interpretation of these data has already been proposed in the literature, but this requires, as a previous step, the identification of patient-specific model parameters. In this paper, a model-based approach is proposed to reproduce individual heart rate signals acquired during tilt tests. A new physiological model adapted to this problem and coupling the ANS, the cardiovascular system (CVS), and global ventricular mechanics is presented. Evolutionary algorithms are used for the identification of patient-specific parameters in order to reproduce heart rate signals obtained during tilt tests performed on eight healthy subjects and eight diabetic patients. The proposed approach is able to reproduce the main components of the observed heart rate signals and represents a first step toward a model-based interpretation of these signals.
1. Introduction
Heart rate variability (HRV) represents one of the most efficient indicators to characterize the modulation of the cardiovascular system (SCV) by the autonomic nervous system (ANS) [1]. In fact, HRV can be easily extracted from the electrocardiogram (ECG), which is a common noninvasive clinical examination reflecting the electrical activity of the heart. However, the interpretation of HRV measurements can be difficult because of the complex mechanisms involved in the autonomic regulation. Time domain and frequency domain methods have been developed to assist the signal analysis and to estimate the levels of vagal (parasympathetic) and sympathetic activites [1]. Although these classical indicators provide useful information and have been widely used in clinical practice, a model-based approach can be particularly useful to complement this information and to ease its interpretation, as these mathematical models directly represent the interactions between the ANS and the CVS [2]. Such a model could also assist in the prediction of the patient’s response to different physiological conditions or therapeutic strategies. However, a necessary step for such model-based interpretation methods is the creation of a patient-specific instance of the model, characterized by individualized parameters.
In this paper, we propose a model-based approach to reproduce patient-specific HRV measurements acquired during tilt tests. After a brief description of the underlying physiology of ANS regulation, a new model of the short-term autonomic regulation of the cardiovascular system is presented in Section 2. The model is constituted of the following subsystems: (i) the cardiac mechanical activity, (ii) the circulatory system and (iii) the autonomic baroreflex loop, including afferent and efferent pathways. The conditions for the simulation of a tilt test using this model and the proposed identification procedure are also described in Section 2. Section 3 presents results of the identification procedure applied to the analysis of tilt tests of healthy and pathologic subjects. Finally, Section 4 presents the conclusions of this work.
2. Methods
2.1. Brief Description of the Underlying Physiology
The cardiovascular system is composed of the heart and two closed systems of vessels known as systemic and pulmonary systems. The primary objective of the CVS is to transport the blood to bring the oxygen from the lung to the organs that need it, and to carry important substances such as hormones and nutriments. The heart is a muscular pump system that pushes blood to all parts of the body. It is divided into four chambers: the two top chambers are called atria, and the lower chambers are called ventricles. The atria collect the blood that enters the heart and push it to the ventricles which eject blood out of the heart into the arteries. These heart chambers alternate periods of relaxation, called diastole, and periods of contraction, called systole.
At the scale of cardiac cells, the contraction is due to the shortening and lengthening of sarcomeres which are the elementary mechanical contractile elements. This mechanical activity is under the influence of an electrical activity, since the variation of the calcium concentration during the action potential (electrical activation of excitable cells) allows the development of force. This variation of force in thecardiac muscle fibres allows the delivering of a sufficient ventricular pressure.
The ANS is responsible for the short-term regulation of the SCV. The baroreflex is initiated by the stimulation of the baroreceptors which are sensory receptors that respond to variations of pressure that are mainly located in the wall of atria, vena cava, aortic arch, and carotid sinus. Other pressure receptors, called cardiopulmonary receptors, are found in veins and atria. Changes in blood pressure are translated into corresponding effects on the efferent sympathetic and parasympathetic pathways. The sympathetic system has a global excitatory effect, increasing heart rate, ventricular contractility, peripheral vascular resistance, and so forth, during situations like hunger, fear, and physical activity. The parasympathetic system presents generally an opposite effect. The main effectors are the heart rate, myocardium contractility, peripheral resistance, and venous blood volume.
The head-up tilt test allows the analysis of a patient’s variations on heart rate and blood pressure during a controlled postural change from a supine to a head-up position. During tilt, approximately 300 to 800 mL of blood may be shifted into the lower extremities, leading to a reduction of venous return and hence of stroke volume. In normal subjects, a decrease in the mean arterial blood pressure (MABP) causes the unloading of arterial baroreceptors, providing a sympathetic activation and a vagal inhibition that leads to an increase on heart rate, ventricular contractility, and peripheral vasoconstriction. A balance is established between heart rate and contractility to maintain the cardiac output and MABP in physiological levels. Some works have shown that slight differences in this response can be observed in patients suffering from diabetes mellitus [1].
2.2. Model Description
The proposed model of the SCV represents the main components of the baroreflex described in the previous section, namely, (i) the ventricles, (ii) the circulatory system, and (iii) the short-term regulation by the ANS. The ventricular model includes a simplified representation of the electromechanical processes involved in the ventricular activity. Moreover, the short-term autonomic regulation of the CVS is taken into account by the modulation of cardiovascular variables (heart rate, etc.). The bond graph formalism can be particularly useful for modelling physiological systems that often include various energy domains, for example, for the design of a bond-graph-based controller for muscle relaxant anesthesia [3] or for modelling the musculoskeletal structure [4]. Models of the vascular system [5–7] are especially interesting, since they take into account different energy phenomena (hydraulic, mechanic, chemical, etc.). the appendix presents the basics of bond graph modelling. A description of each one of the proposed model components is presented in the following sections.
2.2.1. Ventricles
A variety of mathematical models of the ventricular function has been proposed in the literature in order to represent explicitly, at different levels of detail, the cardiac electrical activity [8–10], the excitation-contraction coupling [11–13], the mechanical activity [14, 15], and the mechano-hydraulic coupling [16, 17]. Complete models of ventricular activity are developed from a combination of these different energy domain descriptions [18, 19]. The most detailed approaches represent a fine-grained description of the ventricular activity (at the cellular or subcellular levels), which is necessary for the analysis of regional myocardial dynamics [17, 20, 21]. However, these approaches require significant computational resources and are characterized by an important number of parameters. These aspects reduce the model identifiability, make more difficulty to couple these models with models of other physiological systems (e.g., the ANS) and thus limit their application to our problem. On the other hand, the simplest ventricular models are based on a time-varying elastance [22, 23], which can give realistic simulations of ventricular pressure and volume and require low computational resources. However, the influence of heart rate, calcium concentration variations, and autonomic modulation during the contraction process are not taken into account in these models.
The proposed model can be seen as an improvement of global elastance models and includes a simplified description of the excitation-contraction process. The well-known Beeler and Reuter [24] model of the cardiac action potential (BR model) has been chosen, as it presents a basic description of the intracellular calcium dynamics, while keeping a low level of complexity.
The ventricular mechanical activity is usually described as a function of its active and passive properties. Active properties are under the influence of an electrical activity, since the variation of the calcium concentration during the action potential allows the development of force. Passive properties are mainly due to myocardium organization (fibre orientation, collagen density, etc.). Myocardial tension is usually expressed as the sum of active and passive tensions.
The calcium concentration variable of the BR model is used as input to a model of active mechanical activity. In the present work, the active tension is inspired from the works of Hunter et al. [12] and is defined as where is the value of the tension at the calcium concentration at 50% of the isometric tension, n is the Hill coefficient determining the shape of the curve and is the myofilament “cooperativity.” Assuming that the ventricle is supposed to be made of circumferential muscular fibres, the fibre strain can be expressed as a function of ventricular radius; where r and R are the ventricular radii in the deformed and undeformed states, respectively. In order to obtain the total active force developed by the ventricle, the resulting tension is multiplied by the myocardial wall surface S: Passive tensions are derived from the pole-zero constitutive law presented in [25] to have a relation between the circumferential strain and the tension in the fibre axis.
(i)For axial tension : (ii)For axial compression : The total passive force developed by the ventricle is then calculated by multiplying the passive tension and wall surface S: The total force developed by the ventricle is the sum of active and passive ones: So the force in the fibre axis has been defined as a sum of two functions of ventricular radii. In the Bond Graph formalism, the force developed by the ventricle, which is defined by the previous relations, can be described by two capacitive elements representing the active and passive properties. A 1-junction is used to connect the two capacitive elements because the total force in the fibre axis is defined as the sum of two forces (Figure 1).
During the ventricular contraction, the rise of the force and the variation of the fibre length lead to variations of ventricular pressure. Assuming that the ventricle is supposed to be made of circumferential muscular fibres, the relation between fibre force and ventricular pressure proposed in [6] can be used. In this approach, the ejected volume V is defined as a function of the ventricular radius by where the values of A and n are empirically defined. This approach has shown to provide simulation results that are coherent with physiology [6]. As a result, the change of energy domain from mechanics to hydraulics can be described by a transformer implementing (4).
2.2.2. Circulation
The systemic and pulmonary circulations are composed of different kinds of vessels called arteries, capillaries, and veins. Windkessel models are often used to represent the vascular system [26, 27] by using an electrical analogy. In fact, each part of a vessel can be represented by a set of equations relating its volume V, flow Q, and pressure P.
(i)Capacity : as blood vessels are characterized by their elastic properties, a relation between the volume and the variation of pressure can be defined: , where is the unstretched volume.(ii)Resistance : the resistive properties of vessel are characterized by a relation between the flow and the variation of pressure: .(iii)Inertance : the mass of blood brings some inertial effects through the relation .
Model of vessels can be directly represented in a Bond Graph model by a parallel capacitance, a resistance, and an inertance in series. Models of the whole circulation are defined in a global way by considering groups of vessels (e.g., pulmonary and systemic circulation) in an equivalent lumped model (Figure 2). In order to simulate the effect of a postural change from a supine to a head-up position, the influence of gravity on different parts of the body should be taken into account. The model of the systemic circulation has thus been divided in three parts (the head, the abdomen, and the legs). The heart valves are modelled as nonideal diodes using modulated resistances. The atria are modelled as constant capacitances. The ventricular model described in the previous section is used for the left and right ventricles. Tables 1, 2, 3 and 4 summarizes the values used for the parameters of the circulation model.
2.2.3. Autonomic Nervous System
The existing models of ANS, which are based on a closed-loop representation, can be classified into three main categories.
(i)Behavioral models are based on signal processing and identification theories, such as autoregressive, moving-average models with exogenous input (ARMAX) [35]. Although they allow the reproduction of experimental signals, the physiological interpretation of the parameters of these models is difficult as there is not a direct structural relationship between the physiology and the model components and parameters.(ii)Global models present unphysiological descriptions of the system dynamics [36].(iii)Representation models consist in modelling the different subsystems that can be associated with an entity of the cardiovascular control [37, 38]. Many of them [33, 34] are based on the structure proposed originally by Wesseling and Settels [39], which is composed of delays and first-order filters. This kind of formalism allows the representation of the global neurotransmitter dynamics for a particular efferent pathway and the description of the different response times of the sympathetic and the parasympathetic branches.
A new ANS model is presented to describe the activity of the baroreflex and the cardiopulmonary reflex (Figure 3). The baroreceptor input is the arterial pressure (Pa) and its dynamical properties are represented by a first-orderfilter (Figure 4(a)). The cardiopulmonary receptors are represented by the difference between instantaneous and mean venous pressures (Pv and Pvmean, respectively) (Figure 4(b)).
(a)
(b)
Four variables are controlled in the model by means of different efferent pathways: heart rate, cardiac contractility, systemic resistance, and venous volume. Heart rate depends on the action of both the sympathetic and the parasympathetic systems. The cardiac contractility, the systemic resistance, and the venous volume are only on the influence of the sympathetic system. The same structure, based on a normalized function, a delay and a first-order filter, is used for each one of the modelled efferent pathways. The normalized function is the sigmoid input-output relationship defined in [34]: The generic parameter x represents heart rate, contractility, peripheral, and venous volume regulation. The ANS model is coupled to the CVS by injecting in the latter the four previous controlled variables in the following way.
(i)Heart rate: the model of heart rate regulation is composed of two parts for the vagal and the sympathetic pathways (Figure 5(a)). Each branch is composed of a delay (Rs and Rv are the sympathetic and par-asympathetic delays, respectively), and a first-order filter characterized by a gain (Ks and Kv for the sympathic and the vagal gains, respectively) and a time constant (Ts and Tv). The output signal of the heart rate regulation model (Fc) is continuous and is obtained by adding the contributions from the sympathetic and parasympathetic branches and a basal (unmodulated) heart rate (Io). To obtain pulsating blood pressure, an integral pulse frequency modulation (IPFM) model is used as it transforms a continuous input signal into an event series [40] (Figure 5(b)). The input of the IPFM model is the output signal of the heart rate regulation model. The output of the IPFM allows the excitation of the model of the electrical activity [24] and each emitted pulse (Fc) brings an augmentation of calcium concentration.(ii)Contractility: the Tref parameter of the active tension of cardiac fibres can be considered to be an indicator of the cardiac contractility. In this sense, the Tref definition is replaced by the sum of the basal value and the output signal of the contractility regulation model: (Figure 6).(iii)Peripheral resistance: the systemic resistances are equal to the sum of a constant value and a component, which depends on the ANS regulation: (Figure 7). Pb and Pcp are, respectively, the baroreceptors and the cardiopulmonary receptors outputs.(iv)Venous volume: the constitutive relation of the venous capacity depends on the unstretched volume . The regulated part of the unstretched volume corresponds to the differences with the basal value: (Figure 8). Pb and Pcp are, respectively, the baroreceptors and the cardiopulmonary receptors outputs.
(a)
(b)
2.3. Simulation of a Tilt Test with the Proposed Model
The tilt test produces significant variations of blood pressure in the different parts of the body. In order to take into account the effect of gravity, different pressure levels are imposed on the hydraulic capacities of the model [26]. In the Bond Graph formalism, these pressures are introduced by a 1-junction (Figure 9(b)): where is the fluid pressure in the supine position and has been defined in [26] by the relation where is the table angle represented as a ramp function from 0 to the maximal angle is the onset of tilt, is time to maximum angle and is the pressure due to gravity, which is equal to where is the fluid density and g is the gravitational constant. The parameter h corresponds to the mean distance between one part of body and the heart level. For an adult, the mean value for h is equal to −30 cm for the upper part, 20 cm for the abdomen, and 80 cm for the legs. So it is possible to determine the variation of the pressure due to the gravity in function of the part of the body that is considered for the upper part, for the abdomen, and for the legs.
(a)
(b)
These results have been confirmed experimentally by measuring the pressure in the finger during three testsrealized sequentially on the same normal subject (Figure 9(a)). The arterial pressure has been measured using the Task Force Monitor (CNSystems, Graz, Austria), which is a non-invasive system for continuous beat-to-beat evaluation of cardiovascular variables. Blood pressure is recorded byusing the finger plethysmography technique and is regularly calibrated by an oscillometric blood pressure measurement with a cuff-based device.
During the first, second, and third tests, the finger sensor is, respectively, placed at the level of the abdomen, in an up-per position and in a lower position. When the finger is located at the level of the abdomen or the upper leg (Figures 9(a)(1), and 9(a)(3)), the blood pressure rises when changing from a supine to a head-up position, because the effect of the gravity brings an augmentation of blood volume in the lower part of the body. Oppositely, the measurement done in an upper position (Figure 9(a)(2)) results in a decreased blood volume at the sensor level during the tilt.
These pressure values are applied by implementing (11) in the circulation model for the three distinct parts defined for the systemic circulation. The introduction of time-varying gravity pressures is at the origin of the nonstationary conditions required to simulate the postural change of the patient.
2.4. Identification Algorithm
Most of the model parameters have been obtained from the literature and are presented in Tables 1, 2, 3, and 4. However, in order to determine patient-specific responses of the HR signal, parameters and Tv have been identified, by minimizing an error function between simulated and experimental heart rate signals. As this error function is not differentiable with respect to the model parameters and can have multiple local optima, evolutionary algorithms (EA) have been applied for parameter identification. EA are stochastic search techniques, inspired by the theories of evolution and natural selection, which can be employed to find an optimal configuration for a given system within specific constraints [41]. In these algorithms, each “individual" of a “population" is characterized by a set of parameters (or chromosomes). An initial population is created, usually from a set of random chromosomes, and this population will “evolve," improving its global performance, by means of an iterative process. During this process, each individual is evaluated by means of an error function, and a new generation is produced by applying mutation and crossover operators on selected individuals that present low error values, with probabilities pm and pc, respectively. Convergence and robustness properties of EA have been largely studied in the literature [42–44]. These properties depend upon (i) adequate individual coding, (ii) proper definition of the error function, and (iii) selection of appropriate genetic operators for crossover and mutation.
2.4.1. Individual Representation and Initial Population
Each individual represents an instance of the whole model and is characterized by the 6 real-valued parameters to identify ( and Tv). In order to reduce the search space, parameters values were bounded to physiologically plausible intervals that have been defined from physiological knowledge. The identification has been realized on heart rate signals of 200 heartbeats during tilt tests.
2.4.2. Error Evaluation
The selection process consists in associating each individual with a selection probability computed from its error value, which has to be minimized. The error function has been defined as the absolute value of the difference between each sample of experimental and simulated heart rate signals. This criterion has been chosen to obtain simulated strain as close as possible to real signals.
2.4.3. Selection Method
Once the error function has been evaluated for eachindividual, selection is carried out by means of the “Roulette Wheel" method, adapted for function minimization, in which the probability of selecting a given individual depends on the value of its error function divided by the sum of all the error values of the population. Only standard genetic operators, defined for real-valued chromosomes, have been used in this work: “uniform crossover" which creates two new individuals (offspring) from two existing individuals (parents) by randomly copying each allele from one parent or the other, depending on a uniform random variable and “Gaussian mutation," which creates a new individual by randomly changing the value of one allele (selected randomly), based on a Gaussian distribution around the current value.
3. Results and Discussion
Head-up tilt test has been applied to eight healthy subjects and eight type 2 diabetic patients. This pathology is characterized by high levels of sugar in the blood (hyperglycemia) due to metabolic disorders on the body’s response to insulin or to insulin deficiency.
For each tilt test, classical indicators for HRV analysis were estimated and the model parameters were identified for each subject. Figures 10 and 11 show the acquired heart rate for the eight healthy subjects and the eight diabetic patients during the tilt test (in gray) and the output of the model after parameter identification (in black). In most cases, it is possible to observe that the heart rate increases abruptly after tilt, applied at the end of the first minute of recording, and slowly decreases as the blood pressure approaches its physiological values.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
Simulation results show that the model is able toreproduce the global individual heart rate response to the tilt test for the 16 subjects. The low frequency variations of heart rate are consistent with experimental signals and the first rebound after the heart rate increase is well reproduced in most of the cases. The lack of high frequency components, which can be seen in the experimental heart rate, can be partly explained by the absence of the influence of respira-tion and environmental conditions (such as temperature or external modulations through the central nervous system) on the proposed model. Moreover, in this work, the identification algorithm has been applied in batch mode, in order to obtain one set of parameter values characterizing each patient. This is justified as the identified parameters (related to neurotransmitter densities and efferent temporal delays) are not supposed to vary during a tilt test.
Baroreflex control during the tilt test can be analyzed using the cardiovascular model. Figure 12 illustrates how heart rate is affected by the vagal and the sympatheticpathways for one diabetic patient. Two rebounds (R1 and R2) can be observed both on the experimental heart rate and the simulated signals. The first rebound (R1) is explained by the onset of the vagal inhibition just after the beginning of the tilt test. The sympathetic system is slowly activated and can explain the presence of the second rebound (R2). These estimations of the parasympathetic and sympathetic activities can be useful for a better interpretation of HRV signals during tilt tests.
(a)
(b)
In order to test the robustness of the identification method, the algorithm has been repeated 17 times on one healthy subject. Figure 13 shows boxplots of the values obtained for parameters: and Ts for each of the 17 identification processes. These results show the stability of the solution obtained by the identification procedure.
(a)
(b)
One common time domain indicator for heart rate signal analysis is the standard deviation of the RR intervals between normal beats (SDNN). This indicator reflects the global variability of the heart rate and translates the total power of spectral energy of the HRV signal [1]. The SDNN has been calculated for each subject and compared with the identified model parameters. The correlation coefficient between SDNN and the parasympathetic gain has been calculated and has been found to be equal to 0.9411. Figure 14 shows the linear regression between and SDNN. The correlation between and SDNN can be explained by the presence of the vagal activity both in the low and high frequency components.
4. Conclusion
A new model of the cardiovascular system composed of several coupled subsystems (ventricles, circulation, and ANS) has been proposed and applied to the analysis of heart rate variability signals. The use of Bond Graphs to model ventricular activity and the vascular system has shown to facilitate the coupling of different simple components to form a complex system using the same formalisms for different energetic domains (hydraulic and mechanic). However, although this formalism seems to be particularly adapted to the description of the circulation and global mechanical activity of the ventricles, the marked nonlinearities involved in the genesis of the cardiac action potential has been modeled by a set of ordinary differential equations (the BR model) and coupled with the global Bond Graph model.
The proposed model has been used to reproduce patient-specific heart rate signals during tilt tests on eight healthy subjects and eight diabetic patients. Evolutionary algorithms have been applied as the identification method. After parameter adaptation, simulated signals reproduce the low-frequency components of the observed heart rate signals, which accounts for most of the energy on HRV signals acquired during a tilt test. Moreover, the analysis of identified model parameters has shown a high correlation between the parameter and the global variability calculated on experimental heart rate. In order to reproduce the high frequency components on the HRV signal, new additions to the model should be made, in particular the coupling of the respiratory system. Recursive identification of some model parameters left fixed in this work can be also applied in order to better represent the nonstationarities generated by the tilt test.
A statistical analysis differentiating the two studied populations by means of classical or model-based indicators was out of the scope of this paper, as the number of sub-jects included is still low. So, no conclusions can be done for the moment on the model parameters altered ondiabetic patients. However, this represents one of the future directions of our work. The results presented in this paper areencouraging for the use of this model-based approach in computer-aided diagnosis, and for testing different therapeutic scenarios with a patient-specific model.
Appendix
The Bond Graph Formalism
The Bond Graph (BG) formalism is a diagram-based method that is particularly powerful to represent multienergy systems, as it is based on the representation of power exchanges [45]. Actually, the terminology, the rules, and the construction of Bond Graph models are the same for all energy domains. For example, in the mechanical domain, the effort variable e is the force and the flow variable f is the rate; whereas, in the hydraulic domain, the effort variable e is the pressure and the flow variable f represents flow. The power is the product of the effort and the flow: . The elements of the Bond Graph language can be classified in the following.
(a) Passive elements and :
(i)resistive element : the resistive element is used to describe dissipative phenomenon and can represent electrical resistors, dashpots, or plugs in fluid lines;(ii)capacitive element : the capacitive element is used to describe energy storage and can represent springs or electrical capacitors;(iii)inertial element : the inertial element is used to model inductance effects in electrical systems and mass or inertia effects in mechanical or fluid systems.
(b) Active elements: Se and Sf:
(i)An effort source is an element, which produces an effort, independently of the flow, and a flow source is an element that produces flow independently of the effort.
(c) Junction elements: 0, 1, TF, GY:
(i)0 junction: the 0 junction is characterized by the equality of the efforts on all its links, while the corresponding flows sum up to zero, if power orientations are taken positive toward the junction;(ii)1 junction: 1 junction is characterized by the equality of the flows on all its links, and the corresponding efforts sum up to zero with the same power orientations;(iii)transformer (TF): the transformer TF conserves power and transmits the factors of power with scaling defined by the transformer modulus; it can represent an ideal electrical transformer or a mass-less lever;(iv)gyrator (GY): a gyrator establishes a relationship between flow to effort and effort to flow andconserves the power. It can represent a mechanicalgyroscope or an electrical dc motor.