# An Autonomic Nervous System Model Applied to  the Analysis of Orthostatic Tests

Authors: Virginie Le Rolle; Alfredo I. Hernández; Pierre-Yves Richard; Guy Carrault
Journal: Modelling and Simulation in Engineering
 (2008)

---

## 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.


---

## Body

## 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 Section2. 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 [Ca2+]i 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 (1)Ta=Tref[Ca2+]in[Ca2+]in+Ca50n[1+β(λ−1)], where Tref is the value of the tension at λ=1,Ca50 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; (2)λ=rR, 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: (3)Fa=S⋅Ta. Passive tensions Tp 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(λ>1): (4)Tp=k1ε11(a1−ε11)b1[2+b1ε11a1−ε11],whereε11=12(λ2−1).(ii)
For axial compression(λ<1): (5)Tp=−2⋅k2ε22(a2−ε22)b2[2+b2ε22a2−ε22],whereε22=12(1λ−1). The total passive force developed by the ventricle is then calculated by
multiplying the passive tension and wall surface S: (6)Fp=S⋅Tp. The total force developed by the ventricle is the sum of active and
passive ones: (7)F=Fa+Fp. 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).Figure 1
Bond graph model of
the left ventricle. The calcium concentration variable of the BR model is used
as input to a model of the active mechanical activity, which is described by active and passive capacitive elements implementing (3) and (6), respectively. A 1-junction (7) relates active and passive capacitances. This junction is
linked to a transformer that describes the change of energy domain from
mechanics to hydraulics, as presented in (8).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 (8)V=Arn, 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(C): as blood vessels are characterized
by their elastic properties, a relation between the volume and the variation of
pressure can be defined: V=V0+CP, where V0 is the unstretched volume.(ii)
Resistance(R): the resistive properties of vessel
are characterized by a relation between the flow and the variation of pressure: Q=ΔP/R.(iii)
Inertance(I): the mass of blood brings some
inertial effects through the relation P=I⋅dQ/dt.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 (Figure2). 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.Table 1
Ventricular model.
Parameters valueSourceVentriculeTranformatorA100[6]N3[6]Active capacityTref940 mmHg[12]B01.45[12]B11.95[12]n_ref4.25[12]pC50_ref5.33[12]B20.31[12]S50 cm2EstimatedPassive capacityk11.5 mmHg[12]a10.22[12]b11[12]k20.45 mmHg[12]a20.41[12]b22.5[12]S50 cm2EstimatedTable 2
Systemic circulation.
Parameter valueSourceUnstretched volumeSourceAortaC0.219 mL/mmHg[27]0 mL[28]R0.067 mmHg⋅s/mL[27]I0.00082 mmHg⋅s2/mL[27]ArteriesHead-armsC0.73 mL/mmHg[29]240 mL[29]AbdomenC0.21 mL/mmHg[30]274 mL[30]LegsC1.28 mL/mmHg[29]151 mL[29]ArteriolsHead-armsR3.6 mmHg⋅s/mL[26]AbdomenR3.307 mmHg⋅s/mL[26]LegsR3.9 mmHg⋅s/mL[26]VeinsHead-armsC8 mL/mmHg[26]650 mL[26]R0.23 mmHg⋅s/mL[31]AbdomenC43.11 mL/mmHg[31]1121 mL[30]R0.038 mmHg⋅s/mL[26]LegsC19 mL/mmHg[26]350 mL[26]R0.3 mmHg⋅s/mL[26]Table 3
Pulmonary circulation.
Parameter valueSourceUnstretched volumeSourcePulmonary arteryC0.605 mL/mmHg[32]0 mL[28]R0.1 mmHg⋅s/mL[33]I0.00075 mmHg⋅s2/mL[27]ArteriesC0.02 mL/mmHg[32]90 mL[26]ArteriolsR0.0338 mmHg⋅s/mL[27]I0.0018 mmHg⋅s2/mL[27]VeinsC39.99 mL/mmHg[27]490 mL[26]R0.01875 mmHg⋅s/mL[27]Left AtriaC6.665 mL/mmHg[27]25 mL[28]Table 4
ANS.
Valeur parametersSourceBaroreceptorsT2 s[33]Heart rateVagal activityKIdentifiedTIdentifiedRIdentifiedax0[34]bx1[34]NxO110[34]Tx−0.04[34]Sympathic activityKIdentifiedTIdentifiedRIdentifiedax0.3[34]bx0.7[34]NxO100[34]Tx0.09[34]ContractilityK1EstiméT10 s[33, 34]R3 s[33, 34]ax0.3[34]bx0.7[34]NxO110[34]Tx0.04[34]Peripheral ResistanceBaroreflexK1EstiméT6 s[33, 34]R3 s[33, 34]ax0.3[34]bx0.7[34]NxO110[34]Tx0.04[34]Cardiopulmonary ReflexK0.06[26]T6 s[33]R3 s[33]Venous volumeBaroreflexK1EstiméT6 s[33, 34]R3 s[33, 34]ax0.3[34]bx0.7[34]NxO110[34]Tx0.04[34]Cardiopulmonary ReflexK22[26]T6 s[33]R3 s[33]Figure 2
Proposed Bond Graph model
of the circulation. Groups of vessels are represented by several Bond Graph
elements (capacitance, a resistance, and an inertance) to represent the
pulmonary and the systemic circulation. The latter has been divided into three
different parts (the head, the abdomen, and the legs).


#### 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 (Figure3). 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)).Figure 3
Diagram showing the
coupling between the model of the ANS and the models of the ventricles and the
circulatory system.Model of the baroreceptors (a) and the cardiopulmonary receptors (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]: (9)Fx(t)=ax+bxeTx(PB(t)−Nx,o)+1. 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: Tref=Tref0+ΔT (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: R=R0+ΔR (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 V0:ΔP=(V−Vo)/C.
The regulated part of the unstretched volume corresponds to the differences
with the basal value: V0=V0basal−ΔV (Figure 8). Pb and Pcp are, respectively, the baroreceptors
and the cardiopulmonary receptors outputs.Model of heart rate
regulation (a) and IPFM model (b).
(a)(b)Figure 6
Model of
contractility regulation.Figure 7
Model of peripheral
resistance regulation.Figure 8
Model of venous volume
regulation.


### 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)): (10)P=P0+Ptilt, where P0 is the fluid pressure in the supine position
and Ptilt has been defined in [26] by the relation (11)Ptilt={Pgravity⋅sin(α(t)),t0≤t≤t0+ttilt,Pgravity⋅sin(αmax),t>t0+ttilt, where α(t) is the table angle represented as a ramp
function from 0 to the maximal angle αmax,t0 is the onset of tilt, ttilt is time to maximum angle and Pgravity is the pressure due to gravity, which is equal
to (12)Pgravity=ρ⋅g⋅h, 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 Pgravity=−20mmHg for the upper part, Pgravity=15mmHg for the abdomen, and Pgravity=60mmHg for the legs.(a) Blood pressure
acquired at the level of the abdomen (1), the legs (2), and the upper body(3). (b)
Effort sources introduced in our model to take into account the gravity effects
induced by the tilt test.
(a)(b)These results have been confirmed experimentally by measuring the
pressure in the finger during three testsrealized sequentially on the same
normal subject (Figure9(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 (Figures9(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 Tables1, 2, 3, and 4. However, in order to determine patient-specific
responses of the HR signal, parameters Rs,Rv,Ks,Kv,Ts, 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 (Rs,Rv,Ks,Kv,Ts, 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.


## 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 [Ca2+]i 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 (1)Ta=Tref[Ca2+]in[Ca2+]in+Ca50n[1+β(λ−1)], where Tref is the value of the tension at λ=1,Ca50 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; (2)λ=rR, 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: (3)Fa=S⋅Ta. Passive tensions Tp 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(λ>1): (4)Tp=k1ε11(a1−ε11)b1[2+b1ε11a1−ε11],whereε11=12(λ2−1).(ii)
For axial compression(λ<1): (5)Tp=−2⋅k2ε22(a2−ε22)b2[2+b2ε22a2−ε22],whereε22=12(1λ−1). The total passive force developed by the ventricle is then calculated by
multiplying the passive tension and wall surface S: (6)Fp=S⋅Tp. The total force developed by the ventricle is the sum of active and
passive ones: (7)F=Fa+Fp. 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).Figure 1
Bond graph model of
the left ventricle. The calcium concentration variable of the BR model is used
as input to a model of the active mechanical activity, which is described by active and passive capacitive elements implementing (3) and (6), respectively. A 1-junction (7) relates active and passive capacitances. This junction is
linked to a transformer that describes the change of energy domain from
mechanics to hydraulics, as presented in (8).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 (8)V=Arn, 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(C): as blood vessels are characterized
by their elastic properties, a relation between the volume and the variation of
pressure can be defined: V=V0+CP, where V0 is the unstretched volume.(ii)
Resistance(R): the resistive properties of vessel
are characterized by a relation between the flow and the variation of pressure: Q=ΔP/R.(iii)
Inertance(I): the mass of blood brings some
inertial effects through the relation P=I⋅dQ/dt.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 (Figure2). 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.Table 1
Ventricular model.
Parameters valueSourceVentriculeTranformatorA100[6]N3[6]Active capacityTref940 mmHg[12]B01.45[12]B11.95[12]n_ref4.25[12]pC50_ref5.33[12]B20.31[12]S50 cm2EstimatedPassive capacityk11.5 mmHg[12]a10.22[12]b11[12]k20.45 mmHg[12]a20.41[12]b22.5[12]S50 cm2EstimatedTable 2
Systemic circulation.
Parameter valueSourceUnstretched volumeSourceAortaC0.219 mL/mmHg[27]0 mL[28]R0.067 mmHg⋅s/mL[27]I0.00082 mmHg⋅s2/mL[27]ArteriesHead-armsC0.73 mL/mmHg[29]240 mL[29]AbdomenC0.21 mL/mmHg[30]274 mL[30]LegsC1.28 mL/mmHg[29]151 mL[29]ArteriolsHead-armsR3.6 mmHg⋅s/mL[26]AbdomenR3.307 mmHg⋅s/mL[26]LegsR3.9 mmHg⋅s/mL[26]VeinsHead-armsC8 mL/mmHg[26]650 mL[26]R0.23 mmHg⋅s/mL[31]AbdomenC43.11 mL/mmHg[31]1121 mL[30]R0.038 mmHg⋅s/mL[26]LegsC19 mL/mmHg[26]350 mL[26]R0.3 mmHg⋅s/mL[26]Table 3
Pulmonary circulation.
Parameter valueSourceUnstretched volumeSourcePulmonary arteryC0.605 mL/mmHg[32]0 mL[28]R0.1 mmHg⋅s/mL[33]I0.00075 mmHg⋅s2/mL[27]ArteriesC0.02 mL/mmHg[32]90 mL[26]ArteriolsR0.0338 mmHg⋅s/mL[27]I0.0018 mmHg⋅s2/mL[27]VeinsC39.99 mL/mmHg[27]490 mL[26]R0.01875 mmHg⋅s/mL[27]Left AtriaC6.665 mL/mmHg[27]25 mL[28]Table 4
ANS.
Valeur parametersSourceBaroreceptorsT2 s[33]Heart rateVagal activityKIdentifiedTIdentifiedRIdentifiedax0[34]bx1[34]NxO110[34]Tx−0.04[34]Sympathic activityKIdentifiedTIdentifiedRIdentifiedax0.3[34]bx0.7[34]NxO100[34]Tx0.09[34]ContractilityK1EstiméT10 s[33, 34]R3 s[33, 34]ax0.3[34]bx0.7[34]NxO110[34]Tx0.04[34]Peripheral ResistanceBaroreflexK1EstiméT6 s[33, 34]R3 s[33, 34]ax0.3[34]bx0.7[34]NxO110[34]Tx0.04[34]Cardiopulmonary ReflexK0.06[26]T6 s[33]R3 s[33]Venous volumeBaroreflexK1EstiméT6 s[33, 34]R3 s[33, 34]ax0.3[34]bx0.7[34]NxO110[34]Tx0.04[34]Cardiopulmonary ReflexK22[26]T6 s[33]R3 s[33]Figure 2
Proposed Bond Graph model
of the circulation. Groups of vessels are represented by several Bond Graph
elements (capacitance, a resistance, and an inertance) to represent the
pulmonary and the systemic circulation. The latter has been divided into three
different parts (the head, the abdomen, and the legs).


### 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 (Figure3). 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)).Figure 3
Diagram showing the
coupling between the model of the ANS and the models of the ventricles and the
circulatory system.Model of the baroreceptors (a) and the cardiopulmonary receptors (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]: (9)Fx(t)=ax+bxeTx(PB(t)−Nx,o)+1. 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: Tref=Tref0+ΔT (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: R=R0+ΔR (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 V0:ΔP=(V−Vo)/C.
The regulated part of the unstretched volume corresponds to the differences
with the basal value: V0=V0basal−ΔV (Figure 8). Pb and Pcp are, respectively, the baroreceptors
and the cardiopulmonary receptors outputs.Model of heart rate
regulation (a) and IPFM model (b).
(a)(b)Figure 6
Model of
contractility regulation.Figure 7
Model of peripheral
resistance regulation.Figure 8
Model of venous volume
regulation.


## 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 [Ca2+]i 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 (1)Ta=Tref[Ca2+]in[Ca2+]in+Ca50n[1+β(λ−1)], where Tref is the value of the tension at λ=1,Ca50 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; (2)λ=rR, 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: (3)Fa=S⋅Ta. Passive tensions Tp 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(λ>1): (4)Tp=k1ε11(a1−ε11)b1[2+b1ε11a1−ε11],whereε11=12(λ2−1).(ii)
For axial compression(λ<1): (5)Tp=−2⋅k2ε22(a2−ε22)b2[2+b2ε22a2−ε22],whereε22=12(1λ−1). The total passive force developed by the ventricle is then calculated by
multiplying the passive tension and wall surface S: (6)Fp=S⋅Tp. The total force developed by the ventricle is the sum of active and
passive ones: (7)F=Fa+Fp. 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).Figure 1
Bond graph model of
the left ventricle. The calcium concentration variable of the BR model is used
as input to a model of the active mechanical activity, which is described by active and passive capacitive elements implementing (3) and (6), respectively. A 1-junction (7) relates active and passive capacitances. This junction is
linked to a transformer that describes the change of energy domain from
mechanics to hydraulics, as presented in (8).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 (8)V=Arn, 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(C): as blood vessels are characterized
by their elastic properties, a relation between the volume and the variation of
pressure can be defined: V=V0+CP, where V0 is the unstretched volume.(ii)
Resistance(R): the resistive properties of vessel
are characterized by a relation between the flow and the variation of pressure: Q=ΔP/R.(iii)
Inertance(I): the mass of blood brings some
inertial effects through the relation P=I⋅dQ/dt.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 (Figure2). 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.Table 1
Ventricular model.
Parameters valueSourceVentriculeTranformatorA100[6]N3[6]Active capacityTref940 mmHg[12]B01.45[12]B11.95[12]n_ref4.25[12]pC50_ref5.33[12]B20.31[12]S50 cm2EstimatedPassive capacityk11.5 mmHg[12]a10.22[12]b11[12]k20.45 mmHg[12]a20.41[12]b22.5[12]S50 cm2EstimatedTable 2
Systemic circulation.
Parameter valueSourceUnstretched volumeSourceAortaC0.219 mL/mmHg[27]0 mL[28]R0.067 mmHg⋅s/mL[27]I0.00082 mmHg⋅s2/mL[27]ArteriesHead-armsC0.73 mL/mmHg[29]240 mL[29]AbdomenC0.21 mL/mmHg[30]274 mL[30]LegsC1.28 mL/mmHg[29]151 mL[29]ArteriolsHead-armsR3.6 mmHg⋅s/mL[26]AbdomenR3.307 mmHg⋅s/mL[26]LegsR3.9 mmHg⋅s/mL[26]VeinsHead-armsC8 mL/mmHg[26]650 mL[26]R0.23 mmHg⋅s/mL[31]AbdomenC43.11 mL/mmHg[31]1121 mL[30]R0.038 mmHg⋅s/mL[26]LegsC19 mL/mmHg[26]350 mL[26]R0.3 mmHg⋅s/mL[26]Table 3
Pulmonary circulation.
Parameter valueSourceUnstretched volumeSourcePulmonary arteryC0.605 mL/mmHg[32]0 mL[28]R0.1 mmHg⋅s/mL[33]I0.00075 mmHg⋅s2/mL[27]ArteriesC0.02 mL/mmHg[32]90 mL[26]ArteriolsR0.0338 mmHg⋅s/mL[27]I0.0018 mmHg⋅s2/mL[27]VeinsC39.99 mL/mmHg[27]490 mL[26]R0.01875 mmHg⋅s/mL[27]Left AtriaC6.665 mL/mmHg[27]25 mL[28]Table 4
ANS.
Valeur parametersSourceBaroreceptorsT2 s[33]Heart rateVagal activityKIdentifiedTIdentifiedRIdentifiedax0[34]bx1[34]NxO110[34]Tx−0.04[34]Sympathic activityKIdentifiedTIdentifiedRIdentifiedax0.3[34]bx0.7[34]NxO100[34]Tx0.09[34]ContractilityK1EstiméT10 s[33, 34]R3 s[33, 34]ax0.3[34]bx0.7[34]NxO110[34]Tx0.04[34]Peripheral ResistanceBaroreflexK1EstiméT6 s[33, 34]R3 s[33, 34]ax0.3[34]bx0.7[34]NxO110[34]Tx0.04[34]Cardiopulmonary ReflexK0.06[26]T6 s[33]R3 s[33]Venous volumeBaroreflexK1EstiméT6 s[33, 34]R3 s[33, 34]ax0.3[34]bx0.7[34]NxO110[34]Tx0.04[34]Cardiopulmonary ReflexK22[26]T6 s[33]R3 s[33]Figure 2
Proposed Bond Graph model
of the circulation. Groups of vessels are represented by several Bond Graph
elements (capacitance, a resistance, and an inertance) to represent the
pulmonary and the systemic circulation. The latter has been divided into three
different parts (the head, the abdomen, and the legs).


## 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 (Figure3). 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)).Figure 3
Diagram showing the
coupling between the model of the ANS and the models of the ventricles and the
circulatory system.Model of the baroreceptors (a) and the cardiopulmonary receptors (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]: (9)Fx(t)=ax+bxeTx(PB(t)−Nx,o)+1. 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: Tref=Tref0+ΔT (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: R=R0+ΔR (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 V0:ΔP=(V−Vo)/C.
The regulated part of the unstretched volume corresponds to the differences
with the basal value: V0=V0basal−ΔV (Figure 8). Pb and Pcp are, respectively, the baroreceptors
and the cardiopulmonary receptors outputs.Model of heart rate
regulation (a) and IPFM model (b).
(a)(b)Figure 6
Model of
contractility regulation.Figure 7
Model of peripheral
resistance regulation.Figure 8
Model of venous volume
regulation.


## 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)): (10)P=P0+Ptilt, where P0 is the fluid pressure in the supine position
and Ptilt has been defined in [26] by the relation (11)Ptilt={Pgravity⋅sin(α(t)),t0≤t≤t0+ttilt,Pgravity⋅sin(αmax),t>t0+ttilt, where α(t) is the table angle represented as a ramp
function from 0 to the maximal angle αmax,t0 is the onset of tilt, ttilt is time to maximum angle and Pgravity is the pressure due to gravity, which is equal
to (12)Pgravity=ρ⋅g⋅h, 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 Pgravity=−20mmHg for the upper part, Pgravity=15mmHg for the abdomen, and Pgravity=60mmHg for the legs.(a) Blood pressure
acquired at the level of the abdomen (1), the legs (2), and the upper body(3). (b)
Effort sources introduced in our model to take into account the gravity effects
induced by the tilt test.
(a)(b)These results have been confirmed experimentally by measuring the
pressure in the finger during three testsrealized sequentially on the same
normal subject (Figure9(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 (Figures9(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 Tables1, 2, 3, and 4. However, in order to determine patient-specific
responses of the HR signal, parameters Rs,Rv,Ks,Kv,Ts, 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 (Rs,Rv,Ks,Kv,Ts, 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.


## 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 (Rs,Rv,Ks,Kv,Ts, 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. Figures10 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.Experimental (gray curve) and simulated (black
curve) heart rate (beat/second) for healthy subjects.
(a)(b)(c)(d)(e)(f)(g)(h)Experimental (gray curve) and simulated (black curve) heart rate
(beat/second) for diabetic patients.
(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. Figure12 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.Comparison between
simulated and experimental heart rates for one diabetic patient (a). Simulation
of the activity of sympathetic and parasympathetic pathways after parameter
optimization (b).
(a)(b)In order to test the robustness of the identification method, the
algorithm has been repeated 17 times on one healthy subject. Figure13 shows boxplots of the values obtained for parameters: Kv,Ks,Rv,Rs,Tv, and Ts for each of the 17
identification processes. These
results show the stability of the solution obtained by the identification
procedure.Boxplotsof parametersKv,Ks,Rv, and Tv (a), and Rs and Ts (b)
obtained after 17 applications of the identification procedure. The vertical
axes represent the parameter value.
(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 Kv has been calculated and has
been found to be equal to 0.9411. Figure 14 shows the linear regression between Kv and SDNN. The correlation between Kv and SDNN can be explained by the
presence of the vagal activity both in the low and high frequency components.Figure 14
Linear regression betweenKv and SDNN for the eight healthy subjects and the eight patients.


## 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 parameterKv 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.



---

*Source: 427926-2008-06-16.xml*