University of Southern California Dissertations and Theses

Modeling autonomic peripheral vascular control

(USC Thesis Other)

Modeling autonomic peripheral vascular control

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content

MODELING AUTONOMIC PERIPHERAL VASCULAR CONTROL

by

Patjanaporn Chalacheva

A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(BIOMEDICAL ENGINEERING)

March 2014

Copyright 2014 Patjanaporn Chalacheva
ii

Dedication
To Mom, Dad, Soh and Saw

iii

Acknowledgements
I would like to thank, first and foremost, my advisor, Dr. Michael C.K. Khoo.
Thank you for your support and encouragement through all these years. I am truly
fortunate to have an opportunity to learn from and work with you. Thank you for always
being so kind, understanding, patient, very patient. While letting me explore the world of
science, you are always there to guide me back whenever I drift a little too far. For
everything you have done for me, Dr. Khoo, I thank you.
I would like to also express my gratitude to my doctoral committee: Drs. Thomas
Coates, Vasilis Marmarelis, John Wood, and Roberta Kato for your valuable suggestions.
Thank you for pointing out various aspects of my work that I have not thought of and for
being readily available, in person and emails, whenever I need your advice despite your
busy schedules.
I would like to acknowledge the USC Biomedical Engineering Department and
its staff for assistances. I would like to give special thanks to Mischal Diasanta, our sweet
graduate student advisor, for always making sure that I meet all my degree requirements.
I thank my fellow Khoo crew: Suvimol ‘Ming’ Sangkatumvong, Wenli Wang,
Flavia Oliveira, Winston Tran, Limei Cheng, Jasmine Thum, Alison Hu, Sadaf
Soleymani, Leonardo Nava-Guerra, John Sunwoo, and Wanwara ‘Toey’ Thuptimdang,
for making my life at USC eventful and enjoyable. Thank you for all your help, valuable
inputs and tireless discussion about research and all.
iv

I am especially grateful to Ming for putting in a good word for me with Dr. Khoo,
giving me an opportunity to get my foot in the door so that I eventually found my place in
the ‘Khool lab’. Thank you also for helping me settle down when I first arrived at USC.
Last, but above all, I am eternally thankful to my dearest family. My heartfelt
thanks goes to my sisters, Pitchayanuch ‘Soh’ and Patchara-on ‘Saw’ Chalacheva, for
always being there for me during the tough times, happy times, fun times. I would like to
give my deepest thanks to Mom and Dad, Supang and Wichai Chalacheva, for your
unwavering love and care. Thank you for giving me an opportunity to broaden my
horizons by sending me to study abroad. You are the source of my strength and
inspiration. Without your guidance and support, I would not be where I am today.
Thank you.

v

Table of Contents
Dedication ......................................................................................................................... ii
Acknowledgements .......................................................................................................... iii
Table of Contents .............................................................................................................. v
List of Tables .................................................................................................................... ix
List of Figures .................................................................................................................... x
Abbreviations ................................................................................................................ xvii
Abstract ............................................................................................................................ xx
Chapter 1. Introduction ................................................................................................... 1
1.1 Overview ................................................................................................................. 1
1.2 Purpose of Investigation ......................................................................................... 4
1.2.1 Motivations ................................................................................................. 4
1.2.2 Specific Aims .............................................................................................. 6
Chapter 2. Literature Review .......................................................................................... 8
2.1 Regulation of Blood Pressure – An Overview ........................................................ 8
2.2 Cardiovascular Variability ...................................................................................... 9
2.2.1 Heart Rate Variability (HRV) ..................................................................... 9
2.2.2 Blood Pressure Variability (BPV) ............................................................ 11
2.2.3 Assessment of Cardiovascular Variability – Spectral Analysis ................ 12
2.3 Assessment of Baroreflex Function ...................................................................... 13
2.4 Sigh-Vasoconstriction Response .......................................................................... 16
2.5 Quantification of Peripheral Vascular Resistance ................................................ 18
2.5.1 Peripheral Arterial Tonometry .................................................................. 18
vi

2.5.2 Laser Doppler Flowmetry ......................................................................... 21
2.5.3 Nexfin ....................................................................................................... 22
2.6 Quantitative Modeling of Baroreflex Control....................................................... 24
2.6.1 Structured (Parametric) Models ................................................................ 25
2.6.2 Minimal (Non-Parametric) Models........................................................... 30
Chapter 3. Methods ........................................................................................................ 42
3.1 Experimental Protocols and Data Processing ....................................................... 43
3.1.1 Experimental Protocols ............................................................................. 43
3.1.2 Data Processing ......................................................................................... 45
3.2 Simulation Model.................................................................................................. 48
3.2.1 Description of the Model .......................................................................... 48
3.2.2 Simulation Procedure ................................................................................ 59
3.2.3 Simulated Interventions ............................................................................ 60
3.3 Minimal Modeling Methodology .......................................................................... 61
3.3.1 Linear Model of Peripheral Vascular Conductance .................................. 61
3.3.2 Nonlinear Model of Peripheral Vascular Conductance ............................ 64
3.3.3 Model Estimation and Optimization ......................................................... 66
3.3.4 Physiological Interpretation of Nonlinear Dynamics ............................... 70
3.3.5 Frequency Domain Representation of the Kernels ................................... 73
3.3.6 Decorrelation of Model Inputs .................................................................. 77
3.3.7 Residual Analysis...................................................................................... 79
3.3.8 Estimation of Simulated Data ................................................................... 80

vii

Chapter 4. Results ........................................................................................................... 83
4.1 Subject Characteristics .......................................................................................... 84
4.2 Simulation Model: Model A ................................................................................. 87
4.2.1 Simulated Interventions ............................................................................ 88
4.3 Respiratory Effect on Peripheral Vascular Conductance ...................................... 92
4.3.1 Indirect vs. Direct Respiratory Effect on GPV ......................................... 92
4.3.2 Respiratory Input Signal ........................................................................... 93
4.3.3 Respiratory Coupling in Non-sigh vs. Sigh Data...................................... 95
4.4 Extended Simulation Model: Model B ................................................................. 97
4.4.1 Incorporating RPC into the Simulation Model ......................................... 97
4.4.2 Simulated Interventions ............................................................................ 99
4.4.3 Continuous vs. Intermittent Sigh-Vasoconstriction Mechanism ............ 102
4.5 Minimal Model Evaluation and Validation ........................................................ 109
4.5.1 Model Prediction Accuracy .................................................................... 109
4.5.2 Decorrelation of the Model Inputs .......................................................... 112
4.5.3 Nonlinear Model Kernel Contributions to GPV Variability ................... 113
4.5.4 Estimation of Simulated Data ................................................................. 115
4.6 Effect of Metabolic Syndrome and OSAS on GPV Control .............................. 127
4.6.1 Estimated Kernels ................................................................................... 127
4.6.2 Physiological Interpretation of Nonlinear Kernels ................................. 131
4.6.3 Effect of Orthostatic Stress on System Gains ......................................... 134
4.7 Effect of Blood Transfusion on GPV Control .................................................... 138
4.7.1 Estimated Kernels ................................................................................... 138
viii

4.7.2 System Gains .......................................................................................... 142
Chapter 5. Discussion ................................................................................................... 145
5.1 Respiratory Effect on Peripheral Vascular Conductance .................................... 146
5.2 Simulation Models .............................................................................................. 150
5.3 Minimal Models .................................................................................................. 152
5.4 Kernel Estimation Errors .................................................................................... 159
5.5 Physiological Interpretation of Nonlinear Dynamics ......................................... 161
5.6 System Gains and Autonomic Functions ............................................................ 164
5.7 Limitations of the Study...................................................................................... 168
5.7.1 Simulation Model.................................................................................... 168
5.7.2 Minimal Models ...................................................................................... 168
5.7.3 Surrogate Measures of TPR .................................................................... 170
Chapter 6. Future Work and Conclusions ................................................................. 174
6.1 Future Work ........................................................................................................ 174
6.1.1 Closed-loop Model of Cardiovascular Variability .................................. 174
6.1.2 Linear Time-varying Model .................................................................... 174
6.1.3 Principal Dynamic Modes....................................................................... 176
6.1.4 Multi-input Minimal Model of RPV ....................................................... 177
6.1.5 Reduction in Parameter Estimate Variability ......................................... 177
6.2 Conclusions ......................................................................................................... 180
References ..................................................................................................................... 183

ix

List of Tables
Table 3.1 Simulation Model A parameters ......................................................................49
Table 4.1 Characteristics of the obese pediatric subjects in control, and metabolic
syndrome and obstructive sleep apnea syndrome groups ................................84
Table 4.2 Characteristics of the sickle cell subjects before and after transfusion ...........86
Table 4.3 Baseline values of the key variables generated by Simulation Model A .........87
Table 4.4 Baseline values of the key variables generated by Simulation Model B .........99

x

List of Figures
Figure 1.1 Mechanisms involved in oscillations in sympathetic nervous activity
and blood pressure (Malpas et al. 2001) .........................................................1
Figure 1.2 Extension of closed-loop minimal model to include effects of stroke
volume variability and baroreflex control of the peripheral vasculature. .......6
Figure 2.1 Schematic diagram of the finger probe (Goor et al. 2004) ...........................18
Figure 2.2 An example of PAT traces in relation to ECG and ABP ..............................20
Figure 2.3 A scatter plot of PAT amplitude and Nexfin SVC during supine, cold
face stimulation, and standing in one subject ...............................................24
Figure 2.4 Schematic block diagram of the closed-loop model of circulatory
control (Belozeroff et al. 2002) .....................................................................35
Figure 2.5 Schematic block diagram of the closed-loop model of heart rate
variability (Belozeroff et al. 2003) ................................................................37
Figure 2.6 Schematic block diagram of the closed-loop diagram of circulatory
control. Surrogate cardiac output (SCO) was employed as an input to
circulatory dynamics (CID) instead of RRI (Chaicharn et al. 2009). ...........40
Figure 3.1 Example of the experimental data from an obese pediatric subject and
the extracted beat-to-beat parameters ...........................................................47
Figure 3.2 Example of the raw perfusion signal (PU) collected from a sickle cell
subject and the extracted beat-to-beat PU and G
PV
.......................................48
Figure 3.3 Simulation Model A: Simulation model of cardiovascular system.
Denoted in blue are continuous variables, green are beat-to-beat
variables, and red are key variables. .............................................................52
Figure 3.4 ABP waveform generated by the Windkessel model and the
corresponding notations of beat-to-beat parameters .....................................53
Figure 3.5 Schematic diagram of the one-input minimal model describing G
PV

fluctuations ....................................................................................................62
Figure 3.6 Schematic diagram of two-input minimal model describing G
PV

fluctuations ....................................................................................................63
xi

Figure 3.7 Schematic diagram of the nonlinear model describing GPV
fluctuations ....................................................................................................64
Figure 3.8 Illustration of how to obtain the simulated results from h
BPC
and h
2BPC
.......71
Figure 3.9 Illustration of how to obtain the simulated results from h
RPC
and h
2RPC
.......72
Figure 3.10 Illustration of how to obtain the simulated results from h
BPC,RPC
.................73
Figure 3.11 Step 1 of the model estimation process using the decorrelated inputs .........78
Figure 3.12 Step 2 of the model estimation process using the decorrelated inputs .........78
Figure 3.13 Step 3 of the model estimation process using the decorrelated inputs .........79
Figure 3.14 Configurations of how TPC, the output of the minimal model, was
obtained from the simulation model .............................................................81
Figure 4.1 Scatter plot of log(OAHI) and log(SI) of the control obese subjects
(blue circles) and metabolic syndrome and obstructive sleep apnea
subjects (red crosses). The dash lines show the median of log(OAHI)
and the median of log(SI) from the total subject pool. .................................85
Figure 4.2 Simulation Model A: simulated data under normal physiological
conditions ......................................................................................................90
Figure 4.3 Simulation Model A: simulated data with interventions: an abrupt
increase in ABP (A), an abrupt vasoconstriction (B), and a sigh (C) ..........91
Figure 4.4 Grouped bar chart of 1-input-model and 2-input-model NMSE (mean
± SE) across transfusion treatment groups. * Denotes p<0.001 for
Holm-Sidak post hoc test. .............................................................................92
Figure 4.5 Representative power spectra of the ILV, MAP, PATamp, the one-
input linear model residuals (dash-dotted magenta line), and the two-
input linear model residuals (dashed green lines). ........................................93
Figure 4.6 Grouped bar chart of NMSE of the two-input linear model using ILV,
V
T
and T
E
as the respiratory inputs (mean ± SE) across transfusion
treatment groups. * Denotes p<0.05 for Holm-Sidak post hoc test. .............94
Figure 4.7 Estimated h
RPC
from non-sigh and sigh segments of 7 sickle cell
subjects before receiving blood transfusion ..................................................95
xii

Figure 4.8 Grouped bar chart of the first ten seconds of the h
RPC
estimated from
non-sigh and sigh segments of the sickle cell subjects (mean ± SE). *
Denotes significant difference compared to non-sigh h
RPC
at 1 second
(p<0.05). ** Denotes significant difference compared to non-sigh h
RPC

at 1 and 2 second (p<0.05).
†
Denotes significant difference compared
to non-sigh h
RPC
at 1 second (p<0.05). .........................................................96
Figure 4.9 Simulation Model B: Simulation model after incorporating the sigh-
vasoconstriction mechanism (lower right corner). Denoted in blue are
continuous variables, green are beat-to-beat variables, and red are key
variables ......................................................................................................103
Figure 4.10 Simulation Model B: simulated data under normal physiological
conditions ....................................................................................................104
Figure 4.11 Simulation Model B: simulated data with interventions: an abrupt
increase in ABP (A), an abrupt vasoconstriction (B), and a sigh (C) ........105
Figure 4.12 Data during a sigh from the experimental data (A), the Simulation
Model A (B), and the Simulation Model B (C). Note that the y-scale of
the plot in the last row in panel A is different from panel B and C. ...........106
Figure 4.13 Normalized power spectra of data under normal physiological
conditions from the experimental data (A), the Simulation Model A
(B), and the Simulation Model B (C)..........................................................107
Figure 4.14 Data during a sigh from the experimental data (A), the Simulation
Model B with RPC operating continuously (B), and the Simulation
Model B with RPC operating intermittently (C). Note that the y-scale
of the plot in the last row in panel A is different from panel B and C. .......108
Figure 4.15 Grouped bar chart of linear-model and 2
nd
-order-nonlinear-model
NMSE (mean ± SE) across postures. * Denotes p<0.001 for Holm-
Sidak post hoc test. .....................................................................................109
Figure 4.16 Representative power spectra of the PATamp signal (solid blue lines),
the linear model residuals (dashed green lines) and the second-order
nonlinear model residuals (dotted red lines), for supine and standing
postures. From top to bottom: the spectra from subjects with the
largest, moderate and smallest NMSE reduction. .......................................110
Figure 4.17 Grouped bar chart of linear-model and 2
nd
-order-nonlinear-model
NMSE (mean ± SE) across treatment groups. * Denotes p<0.001 for
Holm-Sidak post hoc test. ...........................................................................111
xiii

Figure 4.18 Representative power spectra of the G
PV
signal (solid blue lines), the
linear model residuals (dashed green lines) and the second-order
nonlinear model residuals (dotted red lines), for pre- and post-
transfusion. The top and the bottom rows: the spectra from subjects
with the largest and smallest NMSE reduction. ..........................................112
Figure 4.19 Grouped bar chart of NMSE of the linear model and linear model with
input decorrelation (mean ± SE) across treatment groups. .........................113
Figure 4.20 Grouped bar chart of the ratio of the variance of the model prediction
by each kernel in the 2
nd
-order nonlinear model and the variance of
PATamp (mean ± SE) across postures (supine and standing).
* Denotes significant difference compared to full model (p<0.001).
** Denotes significant difference compared to full model and h
BPCRPC

(p<0.001). ....................................................................................................114
Figure 4.21 Grouped bar chart of the ratio of the variance of the model prediction
by each kernel in the 2
nd
-order nonlinear model and the variance of
G
PV
(mean ± SE) across treatment groups. ** Denotes significant
difference compared to full model and h
RPC
(p<0.05). ...............................115
Figure 4.22 Estimated h
BPC
from the simulated data with different levels of system
noise and measurement noise when TPC was taken before adding
system noise. Thin black tracing shows the simulated (“true”) h
BPC
. .........120
Figure 4.23 Estimated h
BPC
from the simulated data with different levels of system
noise and measurement noise when TPC was taken after adding
system noise. Thin black tracing shows the simulated (“true”) h
BPC
. .........121
Figure 4.24 Estimated h
RPC
from the simulated data with different levels of system
noise and measurement noise when TPC was taken before adding
system noise. Thin black tracing shows the simulated (“true”) h
RPC
. .........122
Figure 4.25 Estimated h
RPC
from the simulated data with different levels of system
noise and measurement noise when TPC was taken after adding
system noise. Thin black tracing shows the simulated (“true”) h
RPC
. .........123
Figure 4.26 Estimated h
BPC
and h
RPC
from simulated data without operating RPC
(left), with continuously operating RPC (middle), and RPC operating
only with the presence of a sigh (right), in comparison to the simulated
impulse responses (thin black lines) ...........................................................124

xiv

Figure 4.27 Estimated kernels (2
nd
-order nonlinear model) from the original
Simulation Model B “data” (left) and from “data” generated by
Simulation Model B with increased TPR baroreflex nonlinearity
(right). .........................................................................................................125
Figure 4.28 Relationship between simulated gain and system gain derived from the
estimated linear kernels at low- and high-frequency ranges. Top row:
TPR baroreflex gain (K
b,TPR
) vs. system gains of h
BPC
. Bottom row:
sigh-vasoconstriction reflex gain (K
RPC
) vs. system gains. of h
RPC
.
Note that the ordinate scales for the low-frequency RPC gains are
larger than the corresponding scales for the high-frequency RPC gains. ...126
Figure 4.29 Average linear BPC kernels estimated from the linear (blue solid
tracings) and nonlinear model (red dashed tracings), corresponding to
supine and standing postures in control and MetS + OSAS groups. ..........129
Figure 4.30 Average linear RPC kernels estimated from the linear (blue solid
tracings) and nonlinear model (red dashed tracings), corresponding to
supine and standing postures in control and MetS + OSAS groups. ..........129
Figure 4.31 Average 2
nd
-order BPC kernels, corresponding to supine and standing
postures in control and MetS + OSAS groups. ...........................................130
Figure 4.32 Average 2
nd
-order RPC kernels, corresponding to supine and standing
postures in control and MetS + OSAS groups. ...........................................130
Figure 4.33 Average 2
nd
-order cross-kernels, corresponding to supine and standing
postures in control and MetS + OSAS groups. ...........................................131
Figure 4.34 The linear (blue circles) and linear+nonlinear (red squares) BPC
steady-state response as a function of MAP step magnitude in control
(left) and MetS + OSAS (right) representative subjects. ............................132
Figure 4.35 The normalized frequency responses of the linear (top) and
linear+nonlinear (bottom) RPC component for different values of tidal
volume and breathing frequency in control (left) and MetS + OSAS
(right) representative subjects. ....................................................................133
Figure 4.36 Peak responses of the BPC-RPC cross-kernel to a 30 mmHg MAP
pulse train triggered at different times during the inspiratory/expiratory
cycles in control (left) and MetS + OSAS (right) representative
subjects. .......................................................................................................133
Figure 4.37 Nonlinear-model gains from supine to standing in both low- and high-
frequency ranges in obese pediatric subjects (N=49) .................................135
xv

Figure 4.38 Linear-model gains from supine to standing in both low- and high-
frequency ranges (median and IQR) across subject groups. .......................136
Figure 4.39 Nonlinear-model gains from supine to standing in both low- and high-
frequency ranges (median and IQR) across subject groups. .......................137
Figure 4.40 Average linear BPC kernels estimated from the linear (blue solid
tracings) and nonlinear model (red dashed tracings), corresponding to
pre- and post-transfusion.............................................................................138
Figure 4.41 Average linear RPC kernels estimated from the linear (blue solid
tracings) and nonlinear model (red dashed tracings), corresponding to
pre- and post-transfusion.............................................................................139
Figure 4.42 Average 2
nd
-order BPC kernels, corresponding to pre- and post-
transfusion. ..................................................................................................141
Figure 4.43 Average 2
nd
-order RPC kernels, corresponding to pre- and post-
transfusion. ..................................................................................................141
Figure 4.44 Average 2
nd
-order cross-kernels, corresponding to pre- and post-
transfusion. ..................................................................................................141
Figure 4.45 Linear-model gains from pre-transfusion to post-transfusion in both
low- and high-frequency ranges (median and IQR). ..................................142
Figure 4.46 Nonlinear-model gains from pre-transfusion to post-transfusion in both
low- and high-frequency ranges (median and IQR). ..................................144
Figure 5.1 Estimated h
BPC
from data generated by the Simulation Model B. Solid
blue trace shows h
BPC
estimated using MAP as the only input to the 1-
input minimal model. Dashed red trace shows h
BPC
estimated using
MAP and ILV as the inputs to the 2-input minimal model. Thin black
trace shows the simulated (“true”) h
BPC
. .....................................................153
Figure 5.2 Changes in R-R interval at different breathing frequencies and tidal
volumes (Hirsch and Bishop 1981; Eckberg 1995) ....................................162
Figure 5.3 Changes in muscle sympathetic nerve activity (left) after delivering
neck suction at different times in the respiratory cycle (Eckberg 1995). ...163
Figure 5.4 Changes in P-P interval after delivering neck suction at different times
in the respiratory cycle (Eckberg 1995) ......................................................163
Figure 5.5 Shift in TPR baroreflex during postural change .........................................165
xvi

Figure 5.6 Muscle sympathetic traffic measured at different diastolic pressure
levels during one breath (Eckberg 1995) ....................................................170
Figure 6.1 Procedure to generate surrogate data using AAFT method ........................178
Figure 6.2 Estimated h
BPC
and h
RPC
of a sickle cell disease subject. Thick red and
blue tracing are the estimated impulse responses using the actual data.
Thin grey lines are the estimations using the “new” output data, which
is generated using the AAFT method. ........................................................179

xvii

Abbreviations
AAFT amplitude-adjusted Fourier transform
ABP arterial blood pressure
ABR arterial baroreflex control of heart rate
ANOVA analysis of variance
ARX autoregressive model with exogenous input
AV atrioventricular
BAUS brachial artery ultrasound scanning
BMI body mass index
BPC baroreflex control of peripheral vascular conductance
BPV blood pressure variability
BRS baroreflex sensitivity
CFS cold face stimulus
CID circulatory dynamics
CO cardiac output
CO
2
carbon dioxide
CPAP continuous positive airway pressure
CV coefficient of variation
DBP diastolic blood pressure
DER direct effect of respiration
ECG electrocardiogram
FFT fast Fourier transform
FSIVGTT frequently sampled intravenous glucose tolerance test
xviii

G
PV
peripheral vascular conductance
HbS sickle hemoglobin
HCT hematocrit
HOMA homeostatic model assessment
HR heart rate
HRV heart rate variability
ILV instantaneous lung volume
IQR interquartile range
LDF laser Doppler flowmetry
MAP mean arterial pressure
MBF Meixner basis function
MDL minimum description length
MER mechanical effect of respiration on blood pressure
MetS metabolic syndrome
MSNA muscle sympathetic nerve activity
NMSE normalized mean squared error
O
2
oxygen
OAHI obstructive apnea hypopnea index
OSAS obstructive sleep apnea syndrome
PAT peripheral arterial tonometry
PATamp amplitude of peripheral arterial tonometry pulse
PATampN normalized amplitude of peripheral arterial tonometry pulse
PDM principal dynamic modes
PP pulse pressure
xix

PPG photoplethysmograph
PU microvascular perfusion
RATP right atrial transmural pressure
RCC respiratory-cardiac coupling
RPC respiratory-peripheral vascular conductance
R
PV
peripheral vascular resistance
RRI R-R interval
RSA respiratory sinus arrhythmia
SA sinoatrial
SaO
2
oxygen saturation
SBP systolic blood pressure
SCD sickle cell disease
SCO surrogate cardiac output
SE standard error
SI insulin sensitivity
SNA sympathetic nervous activity
SV stroke volume
SVC system vascular conductance
SVR system vascular resistance
T
E
expiration time
TPC total peripheral conductance
TPR total peripheral resistance
V
T
tidal volume

xx

Abstract
The regulation of peripheral vascular resistance (R
PV
) is believed to be largely
sympathetically-mediated. Thus assessment of R
PV
control would allow us to infer
valuable information regarding sympathetic nervous activity. Variability in R
PV
is
generally attributed to the baroreflex control of total peripheral resistance (TPR).
Although it is known that respiration affects sympathetic outflow and deep breaths, akin
to sighs, can lead to peripheral vasoconstriction, the respiratory modulation of TPR has
been little studied. In the present study, we utilized noninvasive surrogate measures of
R
PV
to examine the two mechanisms that influence its variability: the baroreflex control
of peripheral vascular resistance and the respiratory-peripheral vascular resistance
coupling. The first surrogate measure was obtained from peripheral arterial tonometry
(PAT). PAT measured the changes in volume at the finger tip, reflecting the
vasoconstriction response as the reduction in its signal amplitude. The other surrogate
measure was obtained from laser Doppler flowmetry, which monitors microvascular
perfusion. The results of this study suggest that R
PV
fluctuations were directly modulated
by respiration rather than through indirect effect of respiratory modulation of arterial
blood pressure (ABP). The simulation model developed based on previous literatures
pertinent to short-term blood pressure regulation could not reproduce the sigh-
vasoconstriction response as observed in the experimental data. The minimal modeling
approach was employed to estimate this respiratory coupling effect, which would be
incorporated into the simulation model. By means of both modeling approaches, we
xxi

demonstrated that only after the direct respiratory modulation mechanism was added to
the simulation model that a similar vasoconstriction response following a sigh could be
reproduced.
The linear and nonlinear dynamics as well as the interaction effect involved in the
modulation of R
PV
through changes in ABP and respiration were investigated in obese
pediatric subjects exposed to orthostatic stress and subjects with sickle cell disease before
and after blood transfusion treatment. In the obese pediatric subject group, we found that
the linear gains of both the TPR baroreflex as well as the respiratory coupling
mechanisms diminished as a result of orthostatic stress. The reduction in these gains
suggests that sympathetic modulation of TPR decreased in spite of a rise in sympathetic
tone. Orthostatic stress was found to lead also to a reduction in the strength of the
nonlinear behavior in obese pediatric subjects. Subjects with more severe degrees of
metabolic syndrome and obstructive sleep apnea syndrome showed larger reduction in
nonlinear TPR baroreflex gain. Transfusion therapy in the sickle cell disease subjects led
to an increase in nonlinear TPR baroreflex gain as well as the interaction between ABP
and respiration.
In conclusion, through a combination of the structured and the minimal modeling
approaches, we have developed an extended model of blood pressure variability that
incorporates the respiratory modulation effect on R
PV.
Taking this respiratory modulation
effect into account is important for achieving accurate TPR baroreflex estimation.
Finally, the system gains derived from the estimated kernels may constitute potentially
useful biomarkers of sympathetic nervous system function.
1

Chapter 1. Introduction
1.1 Overview
Sympathetic nervous activity (SNA) plays a major role in the regulation of the
cardiovascular system. To better understand how SNA regulates bodily functions, many
techniques have been developed to measure/quantify SNA. However, these techniques
often involve invasive measurements, require technical expertise and can be costly. As a
consequence, these techniques are not readily applicable for clinical usage. Since the
changes in total peripheral resistance (TPR), reflected in vasoconstriction response, is
known to be sympathetically mediated (Malpas et al. 2001), we could potentially employ
the detection of the vasoconstriction response as an indicator of changes in SNA.

Figure 1.1 Mechanisms involved in oscillations in sympathetic nervous activity and blood
pressure (Malpas et al. 2001)
2

Variability in TPR is generally attributed to the baroreflex control. While many
models of the baroreflex control of heart rate (HR) exist in the literature, there have been
very few quantitative studies of the baroreflex control of TPR. This is likely because
measuring TPR is not such a trivial task compared to measuring HR. In the present study,
we utilized noninvasive surrogate measures of peripheral vascular resistance (R
PV
)
obtained from the peripheral arterial tonometry (PAT) and laser Doppler flowmetry. PAT
measures, at the finger tip, the volume changes due to fluctuations in blood pressure.
Vasoconstriction would be reflected on the PAT signal as a decrease in its pulse
amplitude, and vice versa. Thus the amplitude of the PAT signal could be employed as a
measure of vasoconstrictive response. The laser Doppler flowmetry measures skin
microvascular perfusion. The vasoconstriction response would be reflected as the
reduction in microvascular perfusion. From our experimental data, we also observed
evident reduction in PAT amplitude and microvascular perfusion following a deep
inspiration, or the so-called sigh-vasoconstriction response. This suggests that R
PV
is not
only influenced by blood pressure but also by respiration. Thus, in addition to
characterizing the baroreflex control of R
PV
, we also seek to identify how respiration
modulates R
PV
– direct modulation on R
PV
or through the respiratory modulation of blood
pressure.
By employing a combination of structured and minimal modeling approaches, we
found that respiration is likely modulating R
PV
directly. Although the modulatory effect
only becomes apparent when there is a sigh, our findings indicate that this respiratory
modulation occurs with or without the presence of a sigh. Applying these models to
3

responses elicited from obese pediatric subjects with varying degrees of metabolic
syndrome (MetS) and obstructive sleep apnea syndrome (OSAS) during orthostatic stress
showed significant reduction in the gains of both the baroreflex control of peripheral
vascular resistance (BPC) as well as the respiratory-peripheral vascular coupling (RPC)
from supine to standing. This suggests that, in these subjects, their sympathetic
modulations on R
PV
decreased during standing in spite of an increase in sympathetic tone.
Additionally, we also investigated the nonlinear dynamics as well as the
interaction effect involved in the modulation of R
PV
through changes in ABP and
respiration in both obese pediatric subjects exposed to orthostatic stress and sickle cell
subjects before and after blood transfusion treatment. In both subject groups, the results
indicate that both BPC and RPC mechanisms exhibited nonlinear behavior. Furthermore,
there existed an interaction between BPC and RPC dynamics. We demonstrated that
including the nonlinear dynamics as well as the interaction effect, in addition to just the
linear dynamics, significantly decreased the model prediction error. Mechanisms with
large contributions to variability in R
PV
were the second-order nonlinear RPC dynamics
and the interaction effect between ABP and respiration term. In the obese pediatric group,
the nonlinear behavior diminished when exposed to orthostatic stress. Subjects with more
severe degrees of MetS + OSAS showed larger attenuation in nonlinear behavior of the
TPR baroreflex control. On the other hand, the blood transfusion treatment in sickle cell
group led to more pronounced nonlinear behavior in BPC as well as the interaction
between BPC and RPC dynamics compared to before receiving transfusion.
4

To conclude, this study demonstrated a potential model of the R
PV
regulation
using PAT and laser Doppler flowmetry measurement. It also illustrates the possibility of
direct influence of respiration on R
PV
but the mechanism behind this modulation remains
to be investigated. Further, BPC and RPC mechanisms exhibited nonlinear behavior as
well as interacted with each other. By including of these nonlinear and interaction
dynamics, the minimal model of R
PV
was able to provide sensitive biomarkers for the
detection of the changes sympathetic nervous system function in obese pediatric subjects
exposed to orthostatic stress as well as sickle cell subjects with transfusion treatment.
1.2 Purpose of Investigation
1.2.1 Motivations
Assessment of SNA has been an important focus of medical research, as
knowledge of sympathetic tone provides information not only about the underlying
autonomic physiology, but also about the clinical state of the subject being tested.
Various techniques have been developed for the assessment of SNA (Grassi and Esler
1999). Examples include measurement of the hemodynamic responses to various stimuli
in regional cardiovascular districts (Parati et al. 1985), measurements of noradrenaline in
urine or plasma (Esler et al. 1988), microneurography (Vallbo et al. 1979), radiotracer
technology for measurement of norepinephrine spillover (Lambert et al. 1998), power
spectral analysis of HR and ABP (Akselrod et al. 1985), and imaging techniques such as
positron emission tomography and single photon emission computed tomography
scanning (Goldstein et al. 1990). Each technique has its own advantages and
5

disadvantages. Microneurography allows direct and continuous measurement of
adrenergic activity in skeletal muscle circulation but this technique involves insertion of
electrode in a nerve, and the signals obtained are highly local or regional. Power spectral
analysis of HR and ABP technique is noninvasive but the interpretation of the low-
frequency variability is controversial as it contains both sympathetic and parasympathetic
influences. The imaging technique can provide an in vivo measure of noradrenaline
uptake in the heart but it, at the same time, is costly and implementation requires a
considerable degree of technical expertise.
In this study, we seek to develop a computational model of R
PV
regulation, and to
use this model to provide information about sympathetic modulation of the peripheral
vasculature. A key feature is that this model would rely only on measurements that can be
obtained noninvasively. We proposed adopting PAT and laser Doppler flowmetry as the
surrogate measures of R
PV
. Previous studies have demonstrated that peripheral
vasoconstriction response could be measured noninvasively by PAT (Schnall et al. 1999;
O'Donnell et al. 2002; Kuvin et al. 2003; Rubinshtein et al. 2010) and laser Doppler
flowmetry (Inwald et al. 1996; Aso et al. 1997; Allen et al. 2002; Mayrovitz and
Groseclose 2002). Interestingly, no study has yet, to our knowledge, employed PAT or
laser Doppler flowmetry measurements as a means of quantitatively assessing the
mechanisms behind the regulation of R
PV
. As well, only a handful of R
PV
regulation
models have, to date, appeared in the literature. Integrating the two, a new method of
assessing SNA and R
PV
regulation that is noninvasive could be developed. The developed
model of R
PV
regulation could later be also incorporated into the larger closed-loop
6

model of cardiovascular variability (Belozeroff et al. 2002; Khoo 2008; Chaicharn et al.
2009) to give a more comprehensive picture of the cardiovascular control system and
subsequently let us gain more insight into the dynamics of the regulation of this system.
The proposed extended model structure is shown in Figure 1.2.

Figure 1.2 Extension of closed-loop minimal model to include effects of stroke volume
variability and baroreflex control of the peripheral vasculature.

1.2.2 Specific Aims
1.2.2.1 To Model the regulation of Peripheral Vascular Resistance
Our first aim is to characterize the control of R
PV
using noninvasive
measurements collected from 1) obese pediatric subjects with varying degrees of
ABR
Dynamics
MAP(t)
Baroreflex
Control of
TPR
TPR(t)
MAP
Noise
Conversion to
“Surrogate Stroke
Volume”
Feature
Extraction
Arterial Blood
Pressure Contour PP(t)
Windkessel
Model
∆SSV(t)
∆RRI(t)
SQ(t)
SSV(t)
RSA
Dynamics
∆V(t)
Respiration
MER
Dynamics
RRI
Noise
RRI(t)
7

metabolic syndrome and obstructive sleep apnea during supine and standing postures, and
2) sickle cell anemia subjects before and after receiving blood transfusion treatment. PAT
and laser Doppler flowmetry measurements were employed as the surrogate
measurements of R
PV
. Through iterations of structured (parametric) modeling and
minimal (non-parametric) modeling, we would characterize the dynamics of R
PV

regulation and investigate possible nonlinear dynamics as well as interaction effect on the
variability of R
PV
. Further, we would develop and refine a closed-loop simulation model
of cardiovascular system with R
PV
regulation that would be able to reproduce
physiological effects observed in the experimental data.
1.2.2.2 To Explore the Effect of Respiration on Peripheral Vascular Resistance
Vasoconstriction as a response to deep inspiration, or a sigh, was observed in
PAT and laser Doppler flowmetry measurements, suggesting that respiration plays a role
in modulating R
PV
. Our aim is to determine whether respiration has a direct effect on R
PV

or whether changes in R
PV
are mediated by the baroreflexes as a consequence of the
respiratory modulation of ABP. We will also investigate whether the mechanism
responsible for vasoconstriction observed during sigh operates continuously but becomes
apparent only during a sigh; or it acts as a reflex that is triggered only by a sigh of
sufficient large amplitude.
8

Chapter 2. Literature Review
2.1 Regulation of Blood Pressure – An Overview
Blood pressure regulation is dependent upon a number of factors. These include
mechanisms such as the baroreflexes, chemoreflexes, hormonal regulation, and renal
regulation. The arterial baroreflex is a negative feedback mechanism that regulates short-
term changes in arterial blood pressure (ABP) by controlling heart rate (HR), cardiac
contractility, vascular resistance, and venous return. A change in ABP is detected by
stretch receptors called (high-pressure) baroreceptors located in the carotid sinuses and
the aortic arch. Impulses then travel along the glossopharyngeal nerve (if arisen in the
carotid sinus) and the vagus nerves (if arisen in the aortic arch) to the medulla where they
modulate the parasympathetic outflow to the heart and the sympathetic outflow to the
heart and the blood vessels (La Rovere et al. 1995; Levy et al. 2007), leading to change in
ABP in response to the initial change in ABP. In addition to the arterial baroreceptors, the
cardiopulmonary baroreceptors, sometimes called low-pressure baroreceptors, regulate
blood pressure by controlling blood volume. The chemoreflex system consists of
chemoreceptors that are sensitive to medullary pH, arterial CO
2
and O
2
concentrations.
For example, an increase in carbon dioxide concentration or exposure to hypoxia leads to
peripheral vasoconstriction. Other mechanisms are more involved in long-term blood
pressure regulation. Hormonal regulation maintains blood pressure releasing hormones
that will affect fluid retention and thus affect plasma volume. Renal regulation maintains
blood pressure through retention and excretion of extracellular fluid. An increase in renal
9

output results in decreased venous return. Thus cardiac output (CO) is decreased and so
does ABP subsequently. In this study, we would focus our attention to short-term ABP
regulation.
2.2 Cardiovascular Variability
While cardiovascular control mechanisms act to maintain homeostasis, blood
pressure and heart rate still continuously fluctuate over time. The sources of these
fluctuations are from both external perturbations as well as neural control mechanisms
opposing and interacting with each other in order to maintain blood pressure at the
desired level (Parati et al. 2006) and thus these fluctuations reflect the performance of our
cardiovascular health. In this section, we would discuss about mechanisms associated
with oscillations in ABP and HR. It has been established that there are two frequency
bands in HR and ABP with autonomic involvement – low-frequency (0.1 Hz oscillations)
and high-frequency (associated with respiration) oscillations.
2.2.1 Heart Rate Variability (HRV)
The autonomic nervous system regulates HR through cardiac sympathetic and
parasympathetic branches. An increase in parasympathetic activity causes a decrease in
HR whereas an increase in sympathetic activity causes an increase in HR as well as heart
contractility. These two branches generally work reciprocally to control HR. The
parasympathetic system operates by releasing acetylcholine (ACh) at the vagal nerve
terminals (Kandel et al. 2000). ACh acts on muscarinic receptors in the cardiocytes of SA
and atrioventricular (AV) nodes of cardiac muscles (Kandel et al. 2000). The effect of the
10

parasympathetic activity has short latency because the ACh directly activates special
potassium channels in the cardiac cell without the need of the second messenger
operation, allowing the rapid opening of the ion channels (Levy et al. 2007). The
sympathetic nerves innervate both nodal regions and the myocardium. The release of
norepinephrine increases the force of contraction by acting on beta-adrenergic receptors,
which activate the cyclic adenosine monophosphate (cAMP) second-messenger system
(Kandel et al. 2000). This causes an increase in the long-lasting calcium channel current
in the muscle. At the same time, the activation of beta-adrenergic receptors also decreases
the threshold for firing the cardiac pacemaker cells in the SA node, which in turn
increases HR (Kandel et al. 2000). These effects of norepinephrine can be potently
reinforced by circulating epinephrine released from the adrenal medulla (Kandel et al.
2000). Because of these cascade operations, the HR response to sympathetic activity has
longer latency as well as more gradual onset compared to the parasympathetic stimulation
(Levy et al. 2007). Baroreflex mechanism is also thought to influence HRV. Changes in
ABP detected by baroreceptor would lead to changes in autonomic activity to the heart,
which then produce changes in HR.
HR is also influenced by respiration through the mechanism called respiratory
sinus arrhythmia (RSA). HR increases during inspiration and decreases during expiration
(Neff et al. 2003; Levy et al. 2007). RSA is largely considered to be parasympathetically-
mediated since the dynamics of sympathetic modulation of the SA node is known to be
too slow to mediate the variations at respiratory frequency (Katona and Jih 1975;
Eckberg 1983). Moreover, vagal blockade with atropine abolishes the high-frequency
11

oscillations in HR (Eckberg 1995; Lanfranchi and Somers 2002). Origin of RSA is
central coupling or respiratory drive to cardiac vagal motor neurons as changes in HR
was observed at the onset of breathhold or right before the termination of a prolonged
breathhold, where there was an absence of respiratory movement (Malpas 2002).
2.2.2 Blood Pressure Variability (BPV)
There has been considerable debate for many decades as to whether low-
frequency oscillations in ABP are central in origin or they reflect a resonance due to
baroreflex feedback. In support of the central oscillator theory, experiments in cats
showed slow oscillations of preganglionic sympathetic activity in the absence of
concomitant changes in ABP and in animals with denervated baroreceptors (Preiss and
Polosa 1974). The competing hypothesis postulates that the combination of the time
delays involved in the baroreflex feedback loop gave rise to the 0.1 Hz oscillations.
deBoer et al. (1987) developed a simulation model using the information regarding time
delays in the system obtained from experimental data and found that the 0.1 Hz
oscillations could be reproduced. In addition, a unifying theory to account for low-
frequency oscillations in ABP was also suggested. BPV could also be influenced by
vasculature. Structural changes within the vasculature or other mediators of vasculature
tone may modulate the strength of the low-frequency oscillations (Malpas 2002).
The respiratory oscillation in blood pressure could be attributed to variation in
intrathoracic pressure. The effect of respiration on blood pressure was reported to be a
combination of the direct pressure transmission and variation in stroke volume (Scharf et
al. 1980; Toska and Eriksen 1993; Yang and Kuo 1999). During inspiration, the
12

intrathoracic pressure decreases. This lowering of the pressure affects the arteries, and
thus through this direct pressure transmission, ABP decreases. Moreover, respiration has
a direct effect on stroke volume. Intrathoracic pressure decreases during inspiration,
allowing thoracic veins to expand and increase venous return towards the right atrium;
thus increases preload. Through Frank-Starling law of the heart, stroke volume (SV)
increases in response to an increase in preload. At the same time, the lungs expand during
inspiration causing a pulmonary blood volume to increase. This subsequently leads to a
decrease in blood flow from the lungs to the left atrium. As a result, there is a reduction
in stroke volume during inspiration. Decrease in SV leads to decrease in CO.
Consequently, reduction in CO leads to reduction in ABP. However, one study (Yang &
Kuo, 1999) using pharmacological interventions to alter cardiac sympathetic activity,
suggests that there may be a significant autonomic component in the respiratory-related
variations in ABP.
2.2.3 Assessment of Cardiovascular Variability – Spectral Analysis
Spectral analysis is one of the methods used to assess the cardiovascular
variability. The power spectrum of HRV or BPV contains information about sympathetic
and parasympathetic modulation. For HRV, the low-frequency power (0.04-0.15 Hz)
reflects a combination of sympathetic and parasympathetic activities while the high-
frequency power (0.15-0.4 Hz) reflects largely parasympathetic activity (Task Force
1996; Berntson et al. 1997). In addition, the ratio of the low-frequency and high-
frequency spectral power of HR is used as an index of sympathovagal balance (Task
Force 1996). The fact that we can infer information about sympathetic and
13

parasympathetic activity from the frequency components is due to the difference in the
dynamics between the sympathetic and parasympathetic modulation of the SA node and
their synaptic transmission. ACh receptors are fast so parasympathetic traffic to the heart
gets transmitted quickly without much “filtering”. Therefore, the SA node acts as an
allpass filter for parasympathetic activity. On the other hand, sympathetic activity is
slower and so the SA node acts as a lowpass filter for sympathetic traffic. As a result, the
sympathetic traffic is filtered such that it is only present in the low-frequency component
while the parasympathetic traffic is not filtered and so it is present in both low- and high-
frequency components. For BPV, the low-frequency power reflects largely sympathetic
activities (Malpas 2002). The power spectra can be computed using autoregressive
approach or fast Fourier transform (FFT). The former method can be applied to non-
stationary data while the later method is restricted to an assumption of stationarity as a
tradeoff for more convenient computation (Eckberg 1995). Concerning the information
derived from the spectral method, one should note that these power indices represent
degrees of autonomic modulations rather than the level of autonomic tone and averages
of modulations do not represent an averaged level of tone (Task Force 1996).
2.3 Assessment of Baroreflex Function
Arterial baroreflex function is generally assessed by imposing externally-induced
changes to the ABP either pharmacologically or mechanically, or by investigating the
effects of spontaneous changes in ABP at baseline. Baroreflex sensitivity (BRS)
characterizes the power of baroreceptors and central nervous system to maintain and
minimize ABP fluctuations by adjusting heart rate and peripheral vascular resistance.
14

BRS is generally evaluated by four methods of measurements: 1) pharmacological
assessment using the vasoactive drugs, 2) the neck chamber technique, 3) the Valsalva
maneuver, and 4) analysis of spontaneous fluctuations in ABP (La Rovere et al. 1995).
Pharmacological assessment using the vasoactive drugs
The concept of this technique is to study how HR changes in response to
activation or deactivation of baroreceptors by inducing vasoactive drugs into our system
to cause fluctuations in ABP. Phenylephrine, a vasoconstrictor drug, is used to increase
the systolic arterial pressure (SBP). The prolongations of the R-R intervals (RRI) after
injection and the preceding changes in SBP are plotted together. The slope of the
regression line is taken as the BRS, quantifying an increase in RRI per a unit increase in
SBP. Conversely, a vasodilator drug such as nitroglycerin can be used to decrease the
SBP and the shortening of the following RRI is used to derive the BRS in a similar
fashion as described above. While this method is invasive, its advantage is that the
baroreceptors are activated or deactivated in the same direction.
The neck chamber technique
This technique activates or deactivates the carotid baroreceptors by applying
positive or negative pressure to the neck region. An increase in neck chamber pressure is
sensed by the carotid baroreceptors which decreases the ABP. This leads to inhibition of
parasympathetic drive to the heart and activation of sympathetic outflow to the blood
vessels. Neck suction can also be applied to generate the opposite effect on ABP and HR.
The advantages of this method are 1) the localized stimulus to the baroreceptor area
allows us to minimize changes in the response due to other mechanisms, and 2) it allows
15

observation of the response on both the heart and the circulation. However, data
collection from this technique is time-consuming and the neck chamber is not
comfortable to wear.
The Valsalva maneuver
The responses to increases and decreases in ABP mediated by both sympathetic
and parasympathetic drives are produced by abrupt transient voluntary elevation of intra-
thoracic and intra-abdominal pressures provoked by straining. The ABP increases and
immediately follows by reduction in HR at the onset of straining. The continuation of
straining causes a fall in ABP and a rise in HR. Upon releasing, ABP further decreases
and thus increases HR. Then an overshoot in ABP follows, which then reduces HR.
While the Valsalva maneuver is a non-invasive technique, it involves stimulations of
several autonomic receptors, causing high variability in the results.
Analysis of spontaneous fluctuations in ABP
This technique relies on spontaneous changes in ABP and RRI at the baseline
(non-perturbed) level to quantify the BRS. The time-domain approach, also known as the
sequence method, identifies sequences of concomitant changes in SBP and RRI. The
slopes of the regression lines fitted to these sequences are taken as the measure of BRS.
The frequency-domain approach originates from the idea that each spontaneous
oscillations in ABP cause oscillations in RRI of similar frequency. Thus, BRS can be
derived from the ratio of the amplitude of the oscillations in RRI and the ABP. Both
time-domain and frequency-domain approaches allow us to quantify BRS non-invasively.
The advantages of the time-domain approach are 1) the causality in the changes of ABP
16

and RRI is taken into account and 2) it allows the non-baroreflex responses to be
separated from the baroreflex responses. The frequency-domain approach does not
discreetly impose causality but the BRS can be obtained with simple data processing and
computation.
2.4 Sigh-Vasoconstriction Response
Vasoconstriction response following deep inspiration was observed during the
experiments. In 1936, Bolton et al. used digit plethysmography to measure changes in
volume of the digits in all four limbs. They found that the decrease in volume occurred 2-
3 seconds after the commencement of the sigh. Through a series of experiments, their
results suggested that the decrease in digital volume was due to sympathetic activation
rather than an alteration in cardiac output. This was because the vasoconstriction
response was absent in the denervated and sympathectomized limb while remain present
in a normal limb occluded by pressure cuff. From this experiment, it could be inferred
that the decrease in digital volume was produced by vasoconstriction. They also found
that the case of the vasoconstriction was likely dependent on inspiratory phase rather than
the gaseous content of the blood and that the expansion of the chest was the important
factor. Since diminution of digital volume was observed during thoracic breathing but
became absent during abdominal breathing, they hypothesized that there may be an
afferent stimulus from the chest wall that causes this vasoconstriction response. Another
study by Ruttkay-Nedecký (1962) also investigated the finger volume following a deep
breath. The measure of the vasoconstrictor response included the area of the downward
deflection of the plethysmographic curve as well as the duration of the response. In this
17

study, the correlations between the size of the deep breath intervals and the area as well
as the duration of the downward deflection were found to be significant.
In 1969, Browse and Harwick (1969) investigated the effect of deep breath on the
venous tone. They found that there was an increase in vein pressure in the occluded hand
which began 2-5 seconds after the deep breath and that it reached the maximum at 25
seconds. The also investigated the effect of restricting the chest movements. When the
chest and diaphragm were restricted and consequently reduced inspiration volume from
2000 mL to 700 mL, the response was only slightly smaller than the unrestricted case.
Further, they found that during the Valsalva and Muller maneuvers, there was an increase
of occluded hand vein pressure, suggesting that stretch receptors in the lungs may not be
the source of the reflex. However, they noted that there were still changes in intrathoracic
pressure so it was possible that instead of stretch receptors, pressure receptors in the
lungs or blood vessels were the source of the reflex. Later, other studies utilized laser
Doppler flowmetry (LDF) (Baron et al. 1996; Inwald et al. 1996; Aso et al. 1997; Galland
et al. 2000; Allen et al. 2002; Rauh et al. 2003) and photoplethysmograph (PPG) (Allen et
al. 2002; Rauh et al. 2003) in measuring the blood flow response following deep
inspiration. They also found that under normal breathing, there was minimal variation of
cutaneous blood flow. However, deep inspiration would result in detectable variations in
skin blood flow.
Deep breath-vasoconstriction response draws growing interest as a method of
assessing the sympathetic activities. This test had been applied to infants and was found
that it would be a promising as well as convenient tool in assessment of autonomic
18

assessment in infancy (Inwald et al. 1996; Galland et al. 2000). Inwald et al. (1996) found
that this deep breath-vasoconstriction response was not present at birth but appeared to
develop in the first few months of life in both term and preterm infants. As well, the
development of this mechanism appeared to be dependent on post natal rather than post
conceptional age. However, a larger longitudinal study in term and preterm infants should
be investigated. Other studies applied deep breath-vasoconstriction test to type 2 diabetic
patients to evaluate their peripheral sympathetic function in their feet (Aso et al. 1997;
Takahashi et al. 1998). Using LDF, they found that the vasoconstriction response to deep
inspiration was significant decreased in diabetic patients compared to healthy subjects
and that the vasoconstriction response was positively correlated with the duration of
diabetes (Aso et al. 1997).
2.5 Quantification of Peripheral Vascular Resistance
2.5.1 Peripheral Arterial Tonometry

Figure 2.1 Schematic diagram of the finger probe (Goor et al. 2004)
19

Peripheral arterial tonometry (PAT) is developed based on a finger
plethysmograph. It is designed to capture a beat-to-beat plethysmographic recording of
the finger arterial pulse wave amplitude (Schnall et al. 1999; Grote et al. 2003; Kuvin et
al. 2003). The PAT finger probe has a sensor cap that transmits pressure filed and
exhibits a clamp-like effect on the surface of the distal phalanx, measuring pulsatile
changes in volume (Kuvin et al. 2003). A schematic diagram of the figure probe is
illustrated in Figure 2.1. The finger probe consists of a system of inflatable latex air cuffs
connected by pneumatic tubes to an inflating device. A constant counter-pressure was
applied through the air cushion to avoid venous pooling and thus avoid venoarteriolar
reflex vasonconstriction responses (Henriksen and Sejrsen 1976; Rubinshtein et al. 2010).
The device was also designed to avoid occlusion of the arterial blood flow. The pulsatile
volume changes of the digit induce pressure changes in the finger cuff, which were
detected by pressure transducers. A decrease in the ABP volume in the distal fingertip
results in a decrease in pulsatile arterial column changes. This change is reflected as a
decrease in the measured PAT signal (Henriksen and Sejrsen 1976; Rubinshtein et al.
2010). However, one should note that PAT measurement is in arbitrary unit. An example
of the PAT measurement in relation to electrocardiogram (ECG) and continuous blood
pressure waveform is shown in Figure 2.2.
PAT has been used in different applications. Schnall et al. (1999) used PAT in
determining sleep apnea by means of measuring the peripheral circulatory response. By
measuring the peripheral vascular responsiveness at the finger tip, which has high density
of α-sympathetic innervation and high degree of blood flow rate lability, they found that
20

Figure 2.2 An example of PAT traces in relation to ECG and ABP
the blood flow patterns were associated with apneas and this may serve as a marker of the
occurrence and severity of the apneas. Another study investigated the effect of upper
airway obstruction and arousal on PAT in obstructive sleep apnea subjects (O'Donnell et
al. 2002). The amplitude of PAT was employed an indicator of vasoconstriction, which
was reported to be associated with an increase in sympathetic nerve activity and sleep-
disordered breathing events that terminates with arousal. They found that PAT amplitude
was significantly reduced as a consequence of airflow obstruction. In addition, periods of
induced airflow obstruction that resulted in arousal would even further decrease the
amplitude of the PAT signal.
Kuvin et al. (2003) assessed the peripheral endothelial functions by looking at the
pulse wave amplitude measured by PAT as abnormalities in pulse wave amplitude have
150 151 152 153 154 155
-5
0
5
ECG (Volt)
150 151 152 153 154 155
60
80
100
120
ABP (mmHg)
150 151 152 153 154 155
-1
0
1
Time (sec)
PAT (au)
21

been previously reported to be associated with atherosclerosis. Their study compared the
measurement of peripheral endothelial function using brachial artery ultrasound scanning
(BAUS) versus PAT. They found that both BAUS and PAT hyperemia ratio was more
impaired in subjects with cardiovascular risk factors and in subjects who were indicative
of having coronary artery disease. Another study used PAT to detect peripheral arterial
vasoconstriction in response to mental stress test as the associations between myocardial
ischemia induced during mental stress and peripheral arterial vasoconstriction have been
reported (Goor et al. 2004). In this study, they defined a reduction in PAT amplitude
more than 20% of the baseline as an indicator for myocardial ischemia. The results from
PAT were correlated to the measurements using equilibrium radionuclide angiography,
being 88% in concordance. Rubinshtein et al. (2010) examined whether endothelial
dysfunction, measured by PAT, could predict late cardiovascular events. To test the
endothelial function, they computed the reactive hyperemia index. This index was a ratio
of the pulse volume during reactive hyperemia, which was induced by pressure cuff, and
the baseline. They found that low reactive hyperemia index was associated with higher
incidence of adverse events and it was particularly effective in predicting significant
symptoms such as chest pain.
2.5.2 Laser Doppler Flowmetry
In 1975, Stern proposed the development of a noninvasive method of differentiate
flow in microvascular compartments by their differences in velocities (Stern 1975; Stern
et al. 1977). The operating principles for the measurement of blood flow are based on the
idea that the monochromatic laser beam scattered and then partially absorbed by the
22

tissue (Nilsson 1984). As the light hits the moving blood cells, the scattering event
generates shift in wavelength according to the Doppler Effect; however, the wavelength
of the light hitting static object would not be altered (Nilsson 1984). The light scatters
coming back from the microvascular bed have gained the Doppler shift. The frequency
and magnitude distribution of shifts in wavelength are related to the velocity and volume
of moving red blood cells in the tissue volume (Nilsson 1984). Use of laser Doppler
flowmetry has been validated by previous studies and is shown to be a useful tool for
microvascular perfusion monitoring (Arvidsson et al. 1988; Lindsberg et al. 1989;
Almond and Wheatley 1992; Eun 1995). Laser Doppler flowmetry has been employed to
measure the skin blood flow in various applications such as measuring vascular response
due to postural change (Jepsen and Gaehtgens 1993), measuring vascular function in
diabetic patients(Hannemann et al. 2002), as well as measuring the vasoconstriction
response to deep inspiration ) (Baron et al. 1996; Inwald et al. 1996; Aso et al. 1997;
Galland et al. 2000; Allen et al. 2002; Rauh et al. 2003). However, similar to PAT, its
measure of microvascular perfusion has no absolute measurement (Rajan et al. 2009).
2.5.3 Nexfin
Nexfin (BMEYE, Amsterdam, the Netherlands) is a noninvasive and continuous
blood pressure monitoring device. It was developed using the same principles as the
Finapres (Eeftinck Schattenkerk et al. 2009), the predecessor of Nexfin. Finapres
employed the method based on the development of the pulsatile unloading of the finger
arterial walls using an inflatable cuff with built-in photophlethysmograph (Imholz et al.
1998). However, unlike Finapres, Nexfin takes into account of wave reflection while
23

travelling from the aorta to the periphery (the SBP increases while DBP decreases due to
the resistance to the blood flow) and thus provides reconstructed brachial ABP (Bogert et
al. 2004; Eeftinck Schattenkerk et al. 2009). Using the pulse contour method, cardiac
output could be estimated from the finger arterial pressure (Wesseling et al. 1993; Truijen
et al. 2012). Previous studies validated the cardiac output estimated by Nexfin against
other existing methods such as pulmonary artery catheter thermodilution method
(Sokolski et al. 2011) and trans-cardiopulmonary thermodilution (Broch et al. 2012) and
showed promising results. Once the cardiac output is estimated, along with the available
ABP, the continuous TPR (denoted as systemic vascular resistance, SVR, by Nexfin) can
be calculated. Note that with the adjustment of the wave reflection measured at the finger,
the calculated SVR provided by Nexfin should ideally represent the systemic resistance
as opposed to the regional/local resistance measured by PAT or laser Doppler flowmetry.
However, one should keep in mind that the Nexfin SVR was calculated from the finger
arterial pressure. Validation of the accuracy in estimating TPR deserves further
investigation. Figure 2.3 shows a scatter plot of amplitude of PAT and Nexfin SVC, the
inverse of Nexfin SVR, from one subject during supine, cold face stimulation in supine
posture, and standing. Each marker on the scatter plot represents the value of Nexfin SVC
and amplitude of PAT at each beat.
In this study, we would employ PAT amplitude and microvascular perfusion
measured by laser Doppler flowmetry as the indicators of peripheral vascular resistance
(R
PV
). A reduction in PAT amplitude and microvascular perfusion is associated with an
increase in R
PV
and vice versa. Thus PAT amplitude and microvascular perfusion are, in
24

effect, the surrogate measures of peripheral vascular compliance (G
PV
), the inverse of R
PV

(G
PV
= 1/R
PV
).

Figure 2.3 A scatter plot of PAT amplitude and Nexfin SVC during supine, cold face
stimulation, and standing in one subject
2.6 Quantitative Modeling of Baroreflex Control
To better understand the physiological systems, mathematical models are often
employed as investigation tools as they allow us to computationally isolate the
mechanisms of interest from the rest of the other factors. Responses to various kinds of
stimuli can also be conveniently simulated when experiments on animals or human
subjects prove to be difficult. This section discusses structured (parametric) and minimal
(non-parametric) models of short-term cardiovascular regulations. Non-parametric
0.8 1 1.2 1.4 1.6 1.8 2
x 10
-3
0
0.5
1
1.5









⋅s dyn
cm
SVR Nexfin
1
5
PAT Amplitude (au)
+
o
Δ
Supine
Cold Face
Standing
25

models, in this case, differ from parametric models in a sense that the transfer functions,
the functions that describe the dynamics of the mechanisms of interest, are determined
from the data. Parametric models, however, are constructed from predetermined systems
of equations based on a priori knowledge about the system.
2.6.1 Structured (Parametric) Models
2.6.1.1 deBoer (1987)
This is a beat-to-beat, closed-loop model of the cardiovascular system designed to
study the spontaneous short-term variability of ABP and HR in humans at rest. The
model consists of four mechanisms: 1) control of HR and TPR by the baroreflex, 2)
properties of the systemic arterial tree represented by Windkessel model, 3) contractile
properties of the myocardium, and 4) mechanical effects of respiration on ABP. Each
mechanism is represented by difference equations. HR is determined by one weighting
factor for SBP of the current beat and another weighting factor for SBP of the previous
beats. This, in effect, simulates fast parasympathetic and slower sympathetic influence.
The weighting factors, also known as baroreflex sensitivity, are determined from
experimental studies. The TPR is also calculated in the similar fashion as HR. The
arterial compliance is assumed to be constant in the Windkessel model. Once HR is
determined, it then influences CO – longer pulse interval (lower HR) leads to increased
filling of the ventricle, which then results in a more forceful contraction. CO, together
with TPR, subsequently determines the ABP value of the following beat and thus closes
the loop. Respiration is assumed to have a direct effect on CO. This effect of respiration
is subsequently carried on to both the ABP and the HR. Spectral analysis of the simulated
26

data suggested that the respiratory sinus arrhythmia (RSA) is parasympathetically
mediated. The 0.1-Hz oscillation in ABP and HR was described as a resonance
phenomenon due to the delay in the sympathetic control loop of the baroreflex. deBoer et
al. pointed out the drawback of this beat-to-beat model that there is no direct relationship
with real time. This could complicate the use of the model if HR is greater than 75
beats/min because the current SBP does not affect the duration of the current pulse
interval but the following one. This problem arises due to the latency of the vagal
baroreflex.
2.6.1.2 Madwed (1989)
Madwed et al. proposed another closed-loop model of cardiovascular regulation.
It was developed to elucidate the low-frequency oscillations in ABP and HR. The model
incorporates four main features: 1) arterial baroreceptor feedback loops, each controls HR
and TPR, 2) vagal, beta-adrenergic, and alpha-adrenergic effector mechanisms which
control HR and TPR, 3) fixed beat-to-beat SV, and 4) Windkessel model, which
represents the peripheral circulation. The effector mechanisms are modeled as low-pass
filters connected in series with time delays. Vagal effector has the shortest time delay of
100 msec and reaches the maximum response instantaneously. β-adrenergic effector has a
time delay of 2.5 seconds and reaches the maximum HR after 7.5 seconds. α-adrenergic
effector is the slowest, having 5-second delay with a response time of 15 seconds. The
baroreflex feedback loops are modeled as constant gains given that the ABP is within the
operating range. Beyond this range, the resulting HR or TPR saturates at the lower or
upper limit. Respiration enters the model through the vagal effector mechanism. HR
27

increases linearly at a rate of 2-3 beats/min over 50 msec during the rising phase of
inspiration and decreases linearly at a rate of 2-3 beats/min over 50 msec during the
falling phase of expiration. The pulse pressure (PP), which is proportional to SV, is fixed
at 40 mmHg during baseline. The low-frequency oscillation is suggested to be due to 1) a
negative-feedback reflex mechanism and 2) a time-delay after the simulation of the reflex
mechanism. The α-adrenergic effector mechanism controlling TPR is the critical element
in determining the observed oscillations in ABP. It is worth noting that the precise
identification and definitions of each element in this model are not very important for the
generation of the low-frequency oscillations, which is the main purpose of this model.
2.6.1.3 Saul (1991)
Saul et al. developed a closed-loop model of cardiovascular regulation that
focuses on three mechanisms: RSA, HR baroreflex and the mechanical effect of
respiration on ABP. Respiration enters the model through two pathways: 1) centrally,
with a direct effect on sympathetic and parasympathetic efferent activities and 2)
mechanically, through the coupling of respiratory flow and the vasculature. The
feedforward effect of HR on ABP is modeled as a constant gain with a fixed delay of
0.42 seconds. The gain is allowed to change with posture. The baroreflex feedback from
ABP to sympathetic and parasympathetic activities is modeled as a constant gain with a
pure delay of 0.3 seconds. This baroreflex feedback is then combined with the central
effect of respiration to obtain the sympathetic and parasympathetic efferent activities to
the SA node. HR is controlled by distinct sympathetic and parasympathetic pathways to
the SA node where their dynamics are modeled as a single-pole lowpass filter. The
28

parasympathetic response is instantaneous while the sympathetic response has a delay of
1.7 seconds. The parameters used in this model were obtained from previous
experimental studies. In this model, the baroreflex does not control TPR, venous return or
cardiac contractility and does not show 0.1-Hz oscillations.
2.6.1.4 Ursino (1998)
This closed-loop model focuses on modeling the short-term ABP control by the
carotid baroreflex in the presence of pulsatile conditions. It consists of three main
components: 1) the vascular system, 2) the heart, and 3) the carotid baroreflex. The
vascular system is divided into eight compartments. Each compartment includes a
hydraulic resistance, a compliance, and an unstressed volume. The inertia effect is only
modeled in large artery compartments where blood acceleration is significant. The heart
model makes no distinction between the right and the left heart. The atrium is presented
by constant compliance and unstressed volume, which in effect neglects the contractile
effect. The ventricle is modeled as a series of time-varying compliance and resistance.
The carotid baroreflex, which is a focus of this model, is divided into different pathways
and effector responses. The afferent pathway, involving the carotid baroreceptors and the
sinus nerve, is represented as a first-order linear differential equation with a static gain
and a rate-dependent gain, connected in series with a sigmoidal static function. The
efferent sympathetic and parasympathetic pathways are represented as a monotonically
increasing and decreasing exponential curve, respectively. The central nervous system
processing time delay is considered as a part of the effector responses. Effector responses
to sympathetic activity are represented by a pure delay, a monotonic logarithmic static
29

function and a linear first-order dynamics. The response to the vagal stimulation includes
a pure delay, a monotonic static function and a first-ordered linear dynamics. Finally, the
heart period is determined by summing the changes in heart period induced by
sympathetic and parasympathetic stimulations and the intrinsic heart period. The
limitation of the model includes the lack of local autoregulation mechanisms which
controls peripheral systemic resistance. The model also determines the contractility of the
heart based solely on the sympathetic activation while some studies suggested that vagal
stimulation could also play an important role. Lastly, this model does not include vagal
afferents, especially cardiopulmonary baroreflex mechanism.
2.6.1.5 van de Vooren (2007)
This closed-loop model is designed to study short-term ABP control during
supine. This model consists of three sections: 1) a hemodynamic section, 2) a
baroreceptor section, and 3) an autonomic control section. The variables in the
hemodynamic section are beat-to-beat variables. The relationships among the variables
are modeled by difference equations. The heart is modeled as a one-chamber Starling
heart. SV is dependent on pulse interval, venous return volume and contractility volume.
PP is directly proportional to SV. DBP is simulated by Windkessel model, whose time
constant is directly associated with the TPR. Lastly, SBP is a sum of PP and DBP. The
baroreceptors are modeled linearly within a threshold and a saturation level. The SBP is
compared with a reference value, which serves as a set point, mimicking the process of
baroreceptor resetting. The variables in the autonomic control section are time-
continuous. SBP is converted into an afferent neural signal, which goes to three effectors:
30

sympathetic HR control, parasympathetic HR control, and sympathetic TPR control.
These effectors are modeled as frequency-dependent blocks with delays and constant
gains. The neural (time-continuous) and the hemodynamic (beat-to-beat) parts are
integrated together by an integral pulse frequency modulation, which simulates cardiac
pacemaker function. This model does not incorporate the effect of respiration. Moreover,
the dynamic control capabilities of cardiac contractility and venous return on CO were
not included.
2.6.2 Minimal (Non-Parametric) Models
2.6.2.1 Baselli (1988)
This model is one of the earliest closed-loop minimal models of cardiovascular
system based on an autoregressive formulation. It consists of feedforward and feedback
pathways representing the interaction between HR and ABP. The feedforward pathway
represents the non-neural effects of HR on ABP while the feedback pathway represents
the baroreflex control. Respiration enters the model as an exogenous input through both
ABP and HR. Lastly, the pressure loop represents all mechanisms involved in ABP
control that are induced by other variables besides HR. Delays are present in the HR
baroreflex pathway as well as in the pressure loop. This model assumes linearity in the
system and all transfer functions are modeled as polynomials of finite length. The
parameter identification is accomplished by means of the generalized least-squares
method. The optimal order of the model is chosen based on the Akaike’s information
criterion. It is worth noting that there was no distinction between sympathetic and
parasympathetic characteristics to the SA node in this model. Respiration has direct
31

effect on both ABP and HR but the results suggest that respiration affects ABP more
directly compared to HR.
2.6.2.2 MIT Group
The MIT group developed a closed-loop model based on autoregressive moving
average technique (Chon et al. 1997; Mullen et al. 1997; Mukkamala et al. 1999). Its goal
is to analyze fluctuations in HR, ABP, and instantaneous lung volume (ILV) around their
mean values by characterizing the transfer functions from experimental data. The
fluctuations in these variables are interrelated through five coupling mechanisms:
circulatory mechanics, HR baroreflex, SA node, ILV to HR (ILV HR), and ILV to ABP
(ILV ABP). The circulatory mechanics transfer function represents the cardiac
contraction and the generation of the ABP waveform. Pulsatile HR is an input to this
block. The HR baroreflex represents the autonomic coupling between the fluctuations in
HR and ABP. The SA node, not identified from the experimental data, is modeled as a
predefined nonlinear integrate-and-fire component. The ILV HR mechanism represents
the RSA, the autonomic coupling between respiration and HR. Finally, the ILV ABP
represents the mechanical effect of respiration on ABP due to alterations in venous return
and the filling of intrathoracic vessels. Moreover, sources of noise are added to the model
through HR and ABP. Noise added to HR represents other influences on HR that are not
caused by respiration or ABP fluctuations. Likewise, noise added to ABP represents other
influences on ABP that are not caused by respiration or HR fluctuations, for example,
fluctuations in peripheral resistance. The aforementioned mechanisms, except for SA
node, are modeled as linear time-invariant autoregressive moving average difference
32

equations and are solved by least squares fitting. The optimal model order is chosen
based on Rissanen’s minimum description length criterion (Rissanen 1982). The impulse
responses can then be computed from the identified model parameters.
This algorithm was applied on data collected from healthy, nonsmoking male
volunteers during supine and standing postures (Mullen et al. 1997). Control data were
collected prior to administration of atropine and propranolol in each posture. The model
was able to capture changes in cardiovascular control state. The impulse responses
representing autonomic mechanisms (HR baroreflex and ILV HR) were abolished by
double blockade whereas the impulse response representing mechanical effect of
respiration on ABP remained relatively unchanged. This model was also applied to data
collected from diabetic patients (Mukkamala et al. 1999). Subjects were divided into
three groups: minimal, moderate and severe autonomic neuropathy groups. The analysis
showed that the impulse responses of the HR baroreflex and the ILV HR mechanism as
well as the power spectrum of the noise source to HR had greater reduction with higher
severity of autonomic neuropathy. However, there was no significant change in
parameters characterizing circulatory mechanics and mechanical effect of respiration on
ABP as well as the power spectrum of the noise source to ABP.
In addition to the linear model, this research group also explored nonlinear system
identification methods in the cardiovascular system (Chon et al. 1997). Nonlinear
identification was of interest because previous studies (Chon et al. 1996; Kanters et al.
1996) have reported that there may be nonlinear components in the mechanisms
responsible for the fluctuations in HR. The second-order, time-invariant Volterra model
33

was used to represent the nonlinear two-input-one-output system. The two inputs are ILV
and ABP and the output is HR. The relationship between the HR and the two inputs are
represented by five kernels. Two kernels represent the linear contributions of each input
to the output. These two kernels can also be referred to as impulse responses. Another
two kernels represent the quadratic contributions of each input to the output. These four
aforementioned kernels are called self-kernels. The last kernel, called cross-kernel,
represents the effect of the second-order nonlinear interaction between the two inputs on
the output. Each kernel can be estimated as a weighted sum of Laguerre basis functions.
This second-order model showed that it was able to capture the contribution of HR
variability at frequencies below 0.08 Hz, which could not be achieved using a linear
model. Further investigation using neural network analysis with polynomial activation
functions to examine possibility of higher-order nonlinear component showed that a
third-order nonlinear model could explain a greater percentage of HR fluctuations than a
linear, a second-order nonlinear, or higher-order nonlinear models.
Mukkamala et al. (2003) proposed two identification algorithms for quantitating
the TPR baroreflex using non-invasive measurements of beat-to-beat ABP and CO. In
both algorithms, the TPR is regulated by arterial baroreflex and cardiopulmonary
baroreflex. The arterial baroreflex relates changes in ABP to changes in TPR while
cardiopulmonary baroreflex relates changes in right atrial transmural pressure (RATP) to
changes in TPR. However, the RATP cannot be measured noninvasively. The measured
CO thereby is used to compute the SV, which is then employed as a surrogate signal for
the RATP. This is under the assumption that small variations in SV can be characterized
34

by fluctuations in RATP through a linear time-invariant relationship, called the heart-lung
unit. To quantitate the TPR baroreflex, Mukkamala et al. (2003) formulated two system
identification algorithms: direct and indirect identification. The direct identification
involves estimating fluctuations of TPR directly from fluctuations in ABP and SV. The
pathway relating SV to TPR includes the cardiopulmonary baroreflex as well as the
inverse heart-lung unit connected in series. The indirect identification approach involves
two mechanisms relating CO and SV to ABP. The CO to ABP coupling encompasses the
arterial baroreflex and the systemic arterial tree dynamic properties. This means that an
increase in CO would cause an initial increase in ABP by systemic arterial tree. This
initial increase then stimulates arterial baroreflex/systemic arterial tree arc to decrease
TPR in order to maintain ABP level. Another pathway relating SV to ABP encompasses
the inverse heart-lung unit, the cardiopulmonary baroreflex, and the arterial baroreflex.
An increase in SV stimulates the cardiopulmonary baroreflex to decrease TPR. This in
turn stimulates the arterial baroreflex/systemic arterial tree arc to increase TPR in order to
maintain ABP level. Simulation results show that the indirect identification is more
reliable than the direct identification. However, it underestimates the static gain, the
integral of the impulse response, of the cardiopulmonary TPR baroreflex.
2.6.2.3 USC Group
The USC group developed a minimal model of cardiovascular control designed to
characterize the dynamic interrelationships among three key variables: respiration, HR,
and ABP (Belozeroff et al. 2002; Belozeroff et al. 2003; Jo et al. 2007; Khoo 2008;
Chaicharn et al. 2009; Sangkatumvong 2011).
35

Figure 2.4 Schematic block diagram of the closed-loop model of circulatory
control (Belozeroff et al. 2002)
A closed-loop linear model of circulatory control, with respiration as an external
input, was first proposed by Belozeroff et al. (2002). The schematic block diagram is
shown in Figure 2.4. The model can be divided into two parts: mechanisms controlling
fluctuations in HR and mechanisms controlling fluctuations in SBP. RRI, the inverse of
HR, are employed to represent fluctuations in HR. Fluctuations in RRI are assumed to be
produced by two mechanisms. The first mechanism is the arterial baroreflex (ABR) –
changes in SBP leads to changes in RRI. The second mechanism is the RSA – changes in
lung volume leads to changes in RRI. The second part of the model deals with
fluctuations in SBP, which is assumed to be produced by another two mechanisms. The
first one is the mechanical effect of respiration on blood pressure (MER) – changes in
intrathoracic pressure as a result of respiration. The other mechanism is called the
circulatory dynamics (CID), which is responsible for the feedforward effect of
fluctuations in RRI on SBP. Fluctuations in RRI lead to variation in CO as a consequence
of the Frank-Starling mechanism and Windkessel runoff and thereby cause fluctuations in
Σ
w
SBP
(t)
ΔSBP(t)
Σ
w
RR
(t)
h
cid
(t)
Circulatory Dynamics
Respiration

ΔV(t)
h
rsa
(t)
RSA dynamics
h
mer
(t)
Mechanical effects of respiration
h
abr
(t)
Baroreflex Dynamics
ΔRR(t)
36

SBP. For each mechanism, appropriate delays are assigned to impose causality in the
model. ABR has a minimum delay of 0.5 seconds or roughly 1 sampling interval. CID is
assumed to have a minimum delay of 1 second to ensure that a change in the current RRI
affects PP and thus SBP in the next beat. MER is assumed to be instantaneous. Lastly, the
non-causal relationship is assumed for the RSA mechanism and so the delay is allowed to
be negative. The explanation for this is that even though the neural modulation of RRI
and the drive to breathe are simultaneous, the mechanical inspiration takes effect later.
The impulse responses which characterize these four mechanisms are estimated using
kernel expansion technique based on the Volterra-Wiener theory of nonlinearity
(Marmarelis 2004). Each impulse response can be represented as a weighted sum of
Laguerre basis functions. The weighting factors, or also known as the expansion
coefficients, of the Laguerre functions can be estimated using least-squares minimization.
The optimal set of model parameters is chosen based on the global search for the
minimum description length (MDL). Furthermore, the optimal solution must satisfy the
condition that cross-correlations between the residual errors and past values of the inputs
are no significantly different from zero. This closed-loop model was applied to data
collected from obstructive sleep apnea syndrome (OSAS) subjects who had and had not
received long-term Continuous Positive Airway Pressure (CPAP) therapy. The data
indicate that long-term CPAP therapy greatly increases the magnitude of the impulse
responses of the autonomically controlled mechanisms: RSA and ABR. It also leads to
reduction in the magnitude of MER as well as CID impulse responses.
37

Figure 2.5 Schematic block diagram of the closed-loop model of heart rate
variability (Belozeroff et al. 2003)
The closed-loop model was modified to only focus on the fluctuations of the RRI
(Belozeroff et al. 2003) as shown in Figure 2.5. The fluctuations in RRI can be explained
by the ABR and the RSA mechanisms as well as the other influences not related to the
two aforementioned mechanisms. This reduced model was used to analyze OSAS
subjects. Data showed that the RSA gain was three times higher in healthy subjects
compared to subjects with OSAS. Likewise, the ABR gain was also higher in healthy
subjects.
Jo et al. (2007) extended the reduced model (Figure 2.5) to incorporate nonlinear
components of the cardiovascular system. The second-order linear time-invariant model
based on the Volterra-Wiener approach was developed to capture the linear and nonlinear
effects of respiration and SBP on RRI variability. It assumes that fluctuations in RRI are
generated from six mechanisms; five of which are autonomic-mediated mechanisms.
Two mechanisms are linear dynamics representing respiratory-cardiac coupling (RCC),
Respiratory 
Heart-Rate
Coupling
(“RSA”)
Arterial
Baroreflex
(“ABR”)
Other
Influences
Changes
In RRI
(ΔRRI)
Changes
Systolic Blood
Pressure (ΔSBP)
Respiration
(ΔV)
38

the coupling of respiration on RRI, and ABR. Another two mechanisms are the nonlinear
second-order effects of respiration on RRI (ILV
2
-RRI) and SBP on RRI (SBP
2
-RRI). The
fifth mechanism characterizes the interaction of respiration and SBP on RRI (ILV*SBP-
RRI). The last mechanism represents all other influences than cannot be explained by the
aforementioned five components. The model was applied to data collected from controls
and OSAS subjects. The ILV
2
-RRI showed slow dynamics involved in the nonlinear
effects of respiration and RRI. Its magnitude increased significantly from wakefulness to
sleep in both subject groups. This suggests that the nonlinear effects of respiration on
RRI play a more important role during sleep. The SBP
2
-RRI shows much fast dynamics
compared to the nonlinear respiration self-kernel. Its magnitude in OSAS group is
remarkably lower than control group. This supports the idea that both sympathetic and
parasympathetic functions involved in the nonlinear arterial baroreflex modulation of HR
are impaired in OSAS subjects. The ILV*SBP-RRI shows that the interaction between
the two inputs is reduced in OSAS subjects. This finding further suggests that these
subjects have impairment in both sympathetic and parasympathetic functions. The
residuals of this second-order model show that they were mostly broad-band, suggesting
that the model is able to capture most of the dynamics. However, it is worth noting that
although the nonlinear model may be able to capture more dynamics in the system, the
model complexity grows rapidly with increasing model order. The complexity includes
difficulty in representing the kernels graphically as well as difficulty in interpreting the
high-order (multi-dimensional) kernels.
39

Later, Chaicharn et al. (2009) employed a closed-loop model structure similar to
the model proposed by Belozeroff et al. (2002) but the assumption of stationarity was
relaxed to allow the model parameters to be time-varying. Fluctuations in RRI were
assumed to be mediated by the ABR and RCC mechanisms. Fluctuations in SBP were
assumed to be mediated by the direct effect of respiration (DER) and by circulatory
dynamics (CID). The DER represents the transmission between fluctuations in ILV and
SBP, which include the mechanical effect due to changes in intrathoracic pressure as well
as the direct effect of respiration on SV via respiratory-driven sympathetic modulation of
HRV. Unlike Belozeroff et al. (2002), surrogate cardiac output (SCO) was used as an
input to CID instead of RRI (Figure 2.6). SCO was defined as PP divided by RRI where
PP was equal to the difference between SBP and DBP. Thus, CID represents the
combined impedance properties of the heart and systemic vasculature.
Based on the previous study (Belozeroff et al. 2002), the minimum delay for
ABR, reflecting the latency in the baroreflex process, was 1.5 seconds. The delay for CID
was 0.5 seconds to ensure that the change in the current cardiac output can affect the
blood pressure of the next beat. The delay for RCC was allowed to be negative while the
effect of DER was assumed to be instantaneous. Moreover, ε
RRI
and ε
SBP
were included in
the model as well. They represent all other influences on RRI and SBP that are not
explained by the model. In estimating the impulse responses, generalized Laguerre
functions called Meixner functions were employed instead of Laguerre functions. This
was because the rising time of Meixner functions can be controlled by the parameter

40

Figure 2.6 Schematic block diagram of the closed-loop diagram of circulatory control.
Surrogate cardiac output (SCO) was employed as an input to circulatory dynamics (CID)
instead of RRI (Chaicharn et al. 2009).
called order of generalization (Asyali and Juusola 2005). The estimation algorithm of the
stationary impulse response was mostly similar to Belozeroff et al. (2002) but with an
additional step: decorrelate the two inputs to increase the orthogonality as this could
improve the estimation accuracy. This was achieved by using an autoregressive model
with exogenous input to filter out the effect of one input from another. Once the optimal
set of parameters (Meixner function order, delays, and order of generalization) for
stationary impulse responses were obtained, they were utilized in the recursive least
squares algorithm to estimate the autoregressive model coefficient at each time step. To
ensure that the estimation at the current time step was less sensitive to remote past input,
the forgetting factor technique was employed. Finally, the time-varying impulse
responses could be computed from the estimated parameters. This time-varying model
41

was used to analyze the cardiovascular consequences in controls and pediatric patients
with moderate to severe OSAS during cold face stimulation (CFS). The CFS activates
diving reflex, producing an increase in sympathetic drive as well as vagal drive. The data
showed that ABR gain increased steadily during CFS in both groups with similar time
course. RCC gain increased more rapidly after 30 seconds following the CFS onset in
control group. CID gain was higher in control during CFS but remained unchanged in
OSAS group. Lastly, DER gain increased progressively in control group during CFS but
only slightly increased above the baseline in OSAS group.
Sangkatumvong (2011) reduced the model structure proposed by Chaicharn et al.
(2009) to include only ABR, RCC, and ε
RRI
. The model was applied to data collected
from controls, non-transfused and transfused sickle cell disease (SCD) patients.
Stationary analysis on baseline data showed that even though the shapes of the impulse
responses were as expected, the RCC and ABR mechanisms did not differ among the
groups. Time-varying analysis on CFS data showed that CFS increased baroreflex
sensitivity and the ABR gain in control was higher than SCD subjects. Similarly, the
RCC gain increased in controls during CFS but remained unchanged in SCD subjects.
However, these observed changes still did not achieve statistical significance and this was
likely due to small sample size and high variability in subject responses.

42

Chapter 3. Methods
To assess the autonomic peripheral vascular functions, we collected noninvasive
physiological measurements from two subject groups – obese pediatric subjects and
sickle cell subjects. Different techniques employed in identifying/extracting autonomic
functions from the collected data are discussed in this chapter. Starting with the
description of the experimental protocols and data processing techniques that are used to
prepare data, the subsequent sections describe the two modeling approaches being
utilized to assess the autonomic functions. The first of these is referred to as the
“structured modeling” approach. Structured models consist of algebraic and differential
equations that represent the physico-chemical processes underlying the phenomena being
described; these models are generally built on “first principles” that reflect our prior
knowledge of the system under study. For instance, the simulation model of
cardiovascular variability that we have constructed is based on knowledge from the
existing literature. In contrast, the minimal modeling approach assumes little about the
underlying system, and characteristics of the system under study are estimated from only
input and output data. Lastly, we will also describe the methods for extracting
biophysical markers, which in turn are used to assess autonomic function.

43

3.1 Experimental Protocols and Data Processing
3.1.1 Experimental Protocols
Obese Pediatric Subjects
Obese pediatric subjects underwent autonomic function tests, metabolic tests, as
well as sleep studies. The study was conducted in a total of 49 subjects (35 males, age:
12.8 ± 0.30 years (mean ± SE), BMI > 95% for age) with mild to moderate OSAS.
Exclusion criteria were diabetes, systemic hypertension, and treatment for OSAS. The
experimental procedures included 1) noninvasive measurements of physiological signals:
respiratory airflow using pneumotachometer, continuous blood pressure using Nexfin
(BMEYE, Amsterdam, Netherlands), peripheral arterial tonometry or PAT (Itamar
Medical, Caesarea, Israel), and electrocardiogram or ECG during supine and standing
postures with a duration of 10 minutes per posture; 2) morning fasting blood samples,
followed by a frequently sampled intravenous glucose tolerance test (FSIVGTT); 3) dual
energy X-ray absorptiometry for adiposity assessment; and 4) polysomnography. The
procedures were conducted at Children’s Hospital Los Angeles. All physiological signals
collected in this study were sampled at a rate of 512 Hz. The signals were ensured to be
stabilized before recording could begin.
Twenty subjects were selected from the aforementioned subject pool for further
investigation of the effects of metabolic syndrome and obstructive sleep apnea (MetS +
OSAS) on the autonomic control of G
PV
(the inverse of R
PV
) – half of which were treated
as control obese subjects while the other half were treated as MetS + OSAS subjects. The
44

clustering the subjects into two groups were identified by their obstructive apnea
hypopnea index (OAHI) and the insulin sensitivity (SI). Insulin sensitivity describes how
sensitive the body is to the effect of insulin. Low insulin sensitivity means higher amount
of insulin is required to lower blood glucose levels. Low insulin sensitivity is found to be
associated with risk of type 2 diabetes (Fukushima et al. 1983). The control obese
subjects had OAHI of 1.56 ± 0.2 events/hr and SI of 6.67 ± 2.20 ×10
-4
min
-1
/ µU/ml. The
MetS + OSAS subjects had OAHI of 5.79 ± 1.15 events/hr and the SI of 1.83 ± 0.17 ×10
-
4
min
-1
/ µU/ml.
Sickle Cell Subjects
To assess the relationship between autonomic control of the G
PV
and the effects of
transfusion, the study was conducted in 14 sickle cell subjects (7 males, age: 18.4 ± 1.56
years) who were in chronic outpatient transfusion program. Each subject made two visits
to the hospital: 1) on the day of but prior to their regular transfusion, and 2) on the day
after transfusion or within 120 hours from the transfusion. It was necessary to have a wait
period between the two visits to allow the blood volume to stabilize after the transfusion,
avoiding abrupt changes in hemodynamic parameters such as cardiac output and blood
pressure that could bias the resulting measurements. In each visit, respiration, ECG,
continuous blood pressure (BMEYE, Amsterdam, Netherlands), and microvascular
perfusion (PU) using the PeriFlux laser Doppler perfusion monitoring system (Perimed,
Jarfalla, Sweden) were recorded in supine position. After all vital signs had stabilized, a
baseline recording of 5 minutes was made before the subject was asked to take 3
45

voluntary deep inspirations or “sighs”, each of which was spaced out by 1 minute. The
recording continued for another 5 minutes after the last sigh.
A subset of these sickle cell subjects were selected for the analysis of the
respiratory modulation of G
PV
during stable (non-sigh) breathing and during voluntary
sighs. A total of 7 subjects were selected as their baseline recordings showed stable
breathing sections without any spontaneous sighs.
3.1.2 Data Processing
The signals that would be utilized in this study were instantaneous lung volume
(ILV), mean arterial pressure (MAP), and G
PV
. Since G
PV
cannot be measured directly, it
would be represented by signal derived from the PAT signal in the obese pediatric dataset
or the PU signal in the sickle cell dataset. ILV was computed by integrating the airflow
signal. For the cardiovascular variables (MAP and G
PV
), the continuous signals were
converted to beat-to-beat signals. To obtain a beat-to-beat signal, R-waves from the ECG
must be identified as they would be used as the landmarks of the beginning and the end
of each beat. The duration between two consecutive R-waves is called the R-R interval
(RRI), the inverse of heart rate (HR). For each RRI i.e. the interval between the n
th
R-
wave and n+1
th
R-wave, the minimum and the maximum of the n
th
arterial blood pressure
(ABP) pulse waveform were identified and these values would correspond to DBP and
SBP, respectively. MAP of the n
th
beat was the average of the ABP waveform from DBP
of the n
th
beat to DBP of the n+1
th
beat.
Similar to how ABP was processed, the minimum and the maximum of the PAT
pulse waveform in each beat interval were identified. PATamp, the amplitude of PAT
46

signal, was defined as the difference between the maximum and the minimum value of
PAT within a beat interval. Further, because PATamp has an arbitrary unit, it was
normalized by its own baseline. The baseline mean and the standard deviation of
PATamp were computed from the 2
nd
minute (60
th
to 120
th
second) of the PAT recording
collected during supine posture. The normalization was achieved by subtracting the
baseline mean from PATamp then divided by the baseline standard deviation. These
baseline values (mean and standard deviation) were computed per individual and were
also used to normalize the respective PATamp collected during standing posture. The
normalized PATamp would be denoted as PATampN. Figure 3.1 (left panel) shows how
the ILV, beat-to-beat MAP and PATamp were extracted from the experimental data.
Next, the beat-to-beat signals were interpolated such that they became continuous
staircase-like signals as shown in Figure 3.1, right panel. Each of these processed signals
was demeaned (subtracting out the mean value) and subjected to 5
th
-order polynomial
trend removal. The ILV signal was then downsampled to 2 Hz while the beat-to-beat
signals were downsampled to 2 Hz by averaging that is for every 0.5 seconds, the
downsampled value representing that duration was the average of all samples within that
0.5-second interval of the original signal.
The n
th
-beat PU was defined as the median of the PU signal from the n
th
R-wave
to the n+1
th
R-wave of the ECG. Similar to PAT, PU has arbitrary unit and thus it was
necessary for the signal to be normalized. Due to rather noisy nature of the PU signal, the
modal value of the beat-to-beat PU was used for normalization. The normalized PU was
achieved by subtracting the modal value from the beat-to-beat PU then divided by the
47

modal value again. To obtain the surrogate measure of G
PV
, the normalized PU was then
divided by MAP. This transformation is analogous to Ohm’s law. If MAP is analogous to
voltage and normalized PU is analogous to current, then the resulting resistance is MAP
divided by normalized PU. Thus, the resulting conductance is the inverse of the resulting
resistance i.e. normalized PU divided by MAP. The beat-to-beat PU was interpolated into
continuous staircase-like signal then downsampled to 2 Hz by averaging. Figure 3.2
illustrates the extraction of the beat-to-beat PU from the raw PU signal (left) and the
corresponding normalized PU and G
PV
(right).

Figure 3.1 Example of the experimental data from an obese pediatric subject and the
extracted beat-to-beat parameters
123 124 125 126 127
-0.5
0
0.5
123 124 125 126 127
820
840
860
880
123 124 125 126 127
96
98
100
123 124 125 126 127
0.8
0.9
1
1.1
123 124 125 126 127
-0.5
0
0.5
123 124 125 126 127
-1
0
1
123 124 125 126 127
60
110
160
123 124 125 126 127
-0.5
0
0.5
1
Airflow ECG ABP PAT
ΔILV (L) RRI (ms) MAP (mmHg) PATamp (au)
Time (sec) Time (sec)
n-1 n+1 n
48

Figure 3.2 Example of the raw perfusion signal (PU) collected from a sickle cell
subject and the extracted beat-to-beat PU and G
PV

3.2 Simulation Model
3.2.1 Description of the Model
The simulation model of cardiovascular system developed in this study operates
in a closed-loop manner with pulsatile condition (Chalacheva and Khoo 2013). This
model would be denoted as Simulation Model A. The following properties of
cardiovascular system are included in the model: 1) Windkessel model, representing the
peripheral circulation; 2) mechanical effect of respiration on ABP (MER); 3) total
peripheral resistance (TPR) baroreflex control; 4) HR baroreflex control; 5) respiratory
sinus arrhythmia (RSA); 6) Sinoatrial (SA) node response characteristics; and 7) Frank-
24.5 25 25.5 26 26.5 27
5
6
7
8
24.5 25 25.5 26 26.5 27
60
100
140
180
24.5 25 25.5 26 26.5 27
110
115
120
125
24.5 25 25.5 26 26.5 27
-0.06
-0.04
-0.02
0
0.02
24.5 25 25.5 26 26.5 27
98
100
102
24.5 25 25.5 26 26.5 27
-6
-4
-2
0
2
x 10
-4
Time (sec) Time (sec)
ECG PU Beat-to-beat PU
Normalized PU MAP σ
PV
49

Starling effect. Respiration, represented by ILV, is treated as an external input. In
addition, sources of random noise are added to ILV, HR, stroke volume (SV) as well as
total peripheral conductance (TPC = 1/TPR), to represent noise in the system. The model
parameters are listed in Table 3.1. The schematic diagram of the model is illustrated in
Figure 3.3.
Table 3.1 Simulation Model A parameters
Arterial blood pressure control
C Arterial compliance 2 mL/mmHg
K
m
Gain of mechanical effect of respiration on
arterial blood pressure
10 mmHg∙sec/L
Baroreflex control of total peripheral resistance
K
bTPR
Gain of baroreflex control of total peripheral
resistance
0.8 sec/mL
TPR
max
Saturation level of total peripheral resistance 2.7 mmHg∙sec/mL
TPR
min
Threshold level of total peripheral resistance 1.3 mmHg∙sec/mL
ABP
bTPR
ABP value at the center of baroreflex control
of total peripheral resistance function
100 mmHg
K
pr
Gain of the smooth muscle response 0.5 mmHg∙sec/mL/Hz
T
pr
Time constant of the smooth muscle response 2 sec
D
pr
Time delay of smooth muscle response 5 sec
CV
TPR
Coefficient of variation in total peripheral
resistance
0.1 -
Baroreflex control of heart rate
K
bHR
Gain of baroreflex control of heart rate 1.4 bpm/mmHg
50

Baroreflex control of heart rate (continued)
HR
max
Saturation level of heart rate 240 bpm
HR
min
Threshold level of heart rate 25 bpm
ABP
bHR
ABP value at the center of baroreflex control
of heart rate function
52 mmHg
Respiratory Sinus Arrhythmia
K
ILV,s
Gain of respiratory Sinus Arrhythmia in
sympathetic branch
0.4 bpm /L
K
ILV,p
Gain of respiratory Sinus Arrhythmia in
parasympathetic branch
2.5 bpm /L
Sinoatrial node response characteristics
K
s
Gain in sinoatrial node to sympathetic
fluctuations
18 bpm/Hz
f
s
Corner frequency of sinoatrial nodal response
to sympathetic fluctuation
0.015 Hz
D
s
Time delay in sympathetic response 1.7 Sec
K
p
Gain in sinoatrial node to parasympathetic
fluctuations
3 bpm/Hz
F
p
Corner frequency of SA nodal response to
parasympathetic fluctuation
0.2 Hz
HR
intrinsic
Intrinsic heart rate 100 bmp
CV
HR
Coefficient of variation in heart rate 0.1 -
Frank-Starling Effect
SV
0
Baseline stroke volume 77 mL
RRI
0
Baseline R-R interval 0.83 sec
K
a
Sensitivity for afterload effect on stroke
volume
-0.1 mL/mmHg
51

Frank-Starling Effect (continued)
Q
m
Mitral flow late in diastole 15 mL/sec
CV
SV
Coefficient of variation in stroke volume 0.1 -
Respiration
F
ILV
Mean breathing frequency 0.3 Hz
V
ILV
Mean instantaneous lung volume 0.5 L
CV
FILV
Coefficient of variation in breathing
frequency
0.1 -
CV
VILV
Coefficient of variation in instantaneous lung
volume
0.1 -

52

Figure 3.3 Simulation Model A: Simulation model of cardiovascular system. Denoted in blue are continuous variables, green are
beat-to-beat variables, and red are key variables.
HR to ΔHR
s
HR
ΔHR
s
HR to ΔHR
p
HR
ΔHR
p
Baroreflex Control of HR
ABP to HR
HR
ABP
dt
d
K
m
−
MER
p
p
p
F 2 π
s
1
K
(s) H
+
=
s
s
s
F 2 π
s
1
K
(s) H
+
=
SA node
D
s
RC
t
e SBP P(t)
−
⋅ =
Windkessel
K
ILV,p
-K
ILV,s
RSA
HR
intrinsic
Noise
HR
ABP to TPR
TPR
ABP
D
pr
s T 1
K
(s) H
pr
pr
pr
+
=
α -adrenergic effector
Baroreflex Control of TPR
ILV
Saturation
SV= SV
0
+
K
a
·DBP
n-1
+
Q
m
·(RRI
n-1
- RRI
0
)
x
1
ILV
z
-1
1-beat Delay
x
1
Frank-Starling Effect
Noise
TPC
ABP
ABP
ILV
ABP
wk
ΔHR
s
ΔHR
p
HR
RRI
RRI
SBP
R RRI
n-1
DBP
n-1
DBP
PP
TPR TPC
C
1
SV
Noise
SV
53

3.2.1.1 Windkessel Model
A two-element Windkessel model is employed to characterize the overall
dynamic properties of the systemic arterial tree. The two elements consist of a resistor
and a capacitor connected in parallel. The resistor represents TPR and the capacitor
represents the arterial compliance. The ABP waveform thus takes a form of an
exponential decay as follows:

( ) ( )
RC
t
wk
n SBP t ABP
−
⋅ = e Equation 3.1
Figure 3.4 illustrates the generated ABP waveform using the above equation and
how the beat-to-beat values are defined.

Figure 3.4 ABP waveform generated by the Windkessel model and the corresponding
notations of beat-to-beat parameters

899.2 899.4 899.6 899.8 900 900.2 900.4 900.6 900.8
100
120
140
DBP(n+1)
DBP(n)
DBP(n-1)
PP(n+1)
PP(n)
PP(n-1)
SBP(n+1)
SBP(n)
SBP(n-1)
Blood Pressure
RRI(n-1) RRI(n)
ECG
54

SBP(n) denotes the systolic blood pressure value of the current beat, R denotes
TPR value determined at the end of the previous beat (n-1
th
beat), C denotes the arterial
compliance, which is assumed to be constant throughout the simulation, and t represents
the sample time point in the time course of the current RRI or beat duration. From this
equation, ABP decays at a rate determined by the time constant (R·C) until it reaches the
end of the RRI. The value that the ABP decays to is taken as the DBP value of the next
beat. PP(n) denotes the pulse pressure value of the current beat, which is the difference
between the SBP and DBP values of the current beat.
3.2.1.2 Mechanical Effect of Respiration on Blood Pressure (MER)
Changes in intrathoracic pressure induce alterations in venous return and the
filling of intrathoracic vessels and heart chambers, which subsequently lead to changes in
ABP. This mechanism is modeled as derivative of ILV and a negative gain, K
m
,
connected in series (Saul et al. 1991). Therefore, the fluctuations in ABP due to
respiration are closely related to respiratory flow rather than volume. The derivative is
implemented as a finite difference equation. Thus, the change in ABP due to mechanical
effect, denoted as ΔABP
ILV
, is defined as follows:
( ) ( ) ( ) ( ) 1 − − − = ∆ t ILV t ILV K t ABP
m ILV

Equation 3.2
3.2.1.3 Total Peripheral Resistance Baroreflex Control
The ABP is mapped to an inverse-sigmoidal-shaped gain to obtain the targeted
TPR. The targeted TPR is a value of TPR that arterial baroreflex would achieve if the
55

pressure were held constant (Madwed et al. 1989). The mapping function is represented
as follows:

( )
( ) ( )
( ) ( )
K
t ABP TPR
K
t ABP ABP
max min
target
central
bTPR
TPR TPR
t TPR
−
−
+








⋅ +
=
e
e
1

Equation 3.3
where TPR
min
represents the threshold level of targeted TPR as ABP becomes very high,
TPR
max
represents the saturation level of targeted TPR as ABP becomes very low,
ABP
bTPR
represents the inflection point or the central point of the inverse sigmoidal
function on the x-axis, and K is a parameter related to the slope at the inflection point. K
is related to the TPR baroreflex gain, K
bTPR
as follows:

bTPR
min max
K
TPR TPR
K
⋅
−
=
4

Equation 3.4
The targeted TPR is connected to the α-adrenergic effector mechanism, which is
modeled as a first-order lowpass filter, and a pure delay of 5 seconds in series. The
transfer function of the α-adrenergic effector mechanism is represented as follows:

( )
s T
K
s H
pr
pr
pr
+
=
1

Equation 3.5
where s denotes the complex variable in Laplace domain, H
pr
(s) represents the transfer
function of the α-adrenergic effector mechanism, K
pr
represents the gain of the smooth
muscle response, and T
pr
represents the time constant of the smooth muscle response.

56

3.2.1.4 Heart Rate Baroreflex Control
This part of the model is modified from Madwed et al. (1989). The ABP is
mapped to an inverse-sigmoidal-shaped gain to obtain the targeted HR. The assumption
is that the gain remains linear in the normal operating value of ABP and becomes lower
as ABP value deviated away from the normal range. The targeted HR is calculated from
the following equation:

( )
( ) ( )
( ) ( )
K
t ABP HR
K
t ABP ABP
max min
target
central
bHR
HR HR
t HR
−
−
+








⋅ +
=
e
e
1

Equation 3.6
where HR
min
represents the threshold level of targeted HR as ABP becomes very high,
HR
max
represents the saturation level of targeted HR as ABP becomes very low, ABP
bHR

represents the inflection point or the central point of the inverse sigmoidal function on the
x-axis, and K is a parameter related to the slope at the central point. K is related to the HR
baroreflex gain, K
bHR
as follows:

bHR
min max
K
HR HR
K
⋅
−
=
4

Equation 3.7
The targeted HR is then mapped to two piecewise linear functions, representing
the sympathetic and parasympathetic gains (Madwed et al. 1989). The function
representing sympathetic gain maps targeted HR value to a positive change in HR
whereas the other function representing parasympathetic gain maps targeted HR value to
a negative change in HR.

57

3.2.1.5 Respiratory Sinus Arrhythmia (RSA)
RSA represents the autonomically mediated coupling between respiration and HR
– HR increases during inspiration and decreases during expiration. Although RSA is
primarily governed by parasympathetic branch, the changes in HR due to RSA in this
model are divided into sympathetic and parasympathetic branches. ILV is translated to
changes in HR by two static gains with sympathetic gain being small and negative while
parasympathetic gain being larger and positive.
3.2.1.6 Sinoatrial (SA) Node Response Characteristics
This part of the model is modified from Saul et al. (1991). The transfer function
between instantaneous neural firing and HR is modeled as a single-pole lowpass filter.
The effect of vagal stimulation on HR is assumed to be instantaneous whereas the
sympathetic stimulation on HR is slower. Thus a pure delay is connected only to the
lowpass filter representing sympathetic response characteristic in series. Both transfer
functions are implemented as follows:

( )
s
s
s
F
s
K
s H
π 2
1 +
=
Equation 3.8

( )
p
p
p
F
s
K
s H
π 2
1 +
=

Equation 3.9
where s denotes the complex variable in Laplace domain; H
s
(s) and H
p
(s) denotes the
transfer functions between HR and sympathetic and parasympathetic activity,
58

respectively; K
s
and K
p
represents the gain in SA nodal response to sympathetic and
parasympathetic fluctuation, respectively; and F
s
and F
p
denotes the corner frequency of
SA nodal response to sympathetic and parasympathetic fluctuations, respectively.
3.2.1.7 Frank-Starling Effect
The control of SV is modified from the beat-to-beat model by Toska et al. (1996).
The SV is modeled as a combined effect of afterload and preload. The afterload is the
pressure the ventricle has to overcome in order to eject blood. Thus, higher pressure
results in greater afterload and subsequently results in smaller the SV. In this case, the
afterload is dependent on DBP of the previous beat and the afterload sensitivity. The
preload effect is modeled as being dependent on the mitral flow late in diastole and the
change in RRI of the previous beat from the baseline RRI. The longer the RRI increases
filling of the ventricle and thus increases the SV. The SV is calculated as follows:
( ) ( ) ( ) ( )
0 0
1 1 RRI n RRI Q n DBP K SV n SV
m a
− − ⋅ + − ⋅ + =
Equation 3.10
where SV(n) represents the SV of the current beat, SV
0
and RRI
0
represent baseline or the
expected value of SV and RRI, respectively; DBP(n–1) and RRI(n-1) represent the DBP
and RRI of the previous beat respectively; K
a
is an afterload sensitivity and assumed to
be constant; and Q
m
is the mitral flow late in diastole and assumed to be constant. The
computed SV is then divided by the arterial compliance, C, to generate PP.
( )
( )
C
n SV
n PP =
Equation 3.11

59

3.2.1.8 Random Noise Sources
Random noise is added to ILV, HR, SV, and TPC to simulate system noise. The
noise level is controlled by a variable called coefficient of variation (CV), which is
dependent on the mean (μ) and the standard deviation (σ) of the signal.

µ
σ
= CV
Equation 3.12
For ILV, both volume and breathing frequency are allowed to vary in a breath-to-
breath basis. For each of the cardiovascular variables, noise is added in a beat-to-beat
basis. All of these noise sources simulate system noise. This means that these random
effects propagate through this closed-loop system.
3.2.2 Simulation Procedure
For each simulation, the model starts with a heartbeat, generating a SV, and thus
PP, based on the initial values of DBP and RRI. SBP is then computed from an initial
value of DBP plus the newly computed PP. The newly computed SBP, RRI, along with
an initial value of TPR enter the Windkessel model to generate one full beat of ABP
waveform.
The sampling frequency of each simulation is set to 20 Hz. For each sampling
point (every 50 msec), all other continuous variables are computed. The output from the
Windkessel model, ABP
wk
(t) is combined with the change in ABP due to MER, ABP
ILV
,
to generate full ABP. ABP is then regulated by both TPR and HR baroreflexes. The HR
barereflex outputs from both sympathetic and parasympathetic branches are combined
with the changes in HR due to RSA. These changes then enter the SA node and the final
60

changes in HR due to sympathetic activities (ΔHR
s
) and parasympathetic activities
(ΔHR
p
) are obtained. Full HR is a sum of intrinsic HR, ΔHR
s
and ΔHR
p
. This concludes
the computations for continuous variables in one beat.
The last value of the continuous HR of the current beat is taken to represent the
value of the current beat in the beat-to-beat HR. The inverse of this current beat HR, i.e.
the current beat RRI, then determines the duration of the next heartbeat in the simulation.
The last value of the continuous ABP of the current beat is taken to represent the beat-to-
beat DBP of the next beat. Finally, the last value of the continuous TPR is taken to
represent the beat-to-beat TPR of the next heartbeat. With these three beat-to-beat initial
variables, the next heartbeat is ready to be generated. The whole process repeats until the
specified simulation time is reached.
ILV signal takes a form of sinusoidal wave and is predetermined before the start
of the simulation. Both the breathing frequency and volume per breath are allowed to
fluctuate randomly where the degree of random fluctuations is controlled by CV. White
noise is also added to the beat-to-beat HR, SV, and TPC to simulate noise in the system.
3.2.3 Simulated Interventions
Three types of interventions are introduced to the simulation model in order to
evaluate whether the simulation model is able to reproduce similar responses to what
have been observed in the experiments. The three interventions include: 1) an abrupt
increase in ABP, 2) an abrupt vasoconstriction, and 3) a deep inspiration or a sigh.
An abrupt increase in ABP is simulated by increasing the value of 3 consecutive
beats of DBP by 25, 50, and 25 mmHg such that
61

( ) ( )
( ) ( )
( ) ( ) 25 2 2
50 1 1
25
0
0
0
+ + = +
+ + = +
+ =
n DBP n DBP
n DBP n DBP
n DBP n DBP

Equation 3.13
where DBP
0
(n) is the supposed value of DBP calculated by the model if there were not
any intervention. The values of SBP are also modified in a similar fashion. These values
of DBP and SBP are modified before entering the Windkessel model.
An abrupt vasoconstriction is simulated by decreasing the value of TPC by 0.5
mL·mmHg
-1
·s
-1
such that
( ) ( ) 5 . 0
0
− = n TPC n TPC
Equation 3.14
where TPC
0
(n) is the supposed value of TPC calculated by the model if there were not
any intervention. The inverse of the modified TPC value, modified TPR, then enters the
Windkessel model.
Lastly, a sigh is simulated by increasing the value of one full breath of ILV by 3
times of the original value if there were not any intervention.
3.3 Minimal Modeling Methodology
3.3.1 Linear Model of Peripheral Vascular Conductance
One-Input Minimal Model
The one-input minimal model of peripheral vascular conductance (G
PV
), the
inverse of peripheral vascular resistance (R
PV
), relates fluctuations in MAP (ΔMAP) to
fluctuations in G
PV
(ΔG
PV
) through baroreflex control of peripheral vascular conductance
mechanism (BPC). This model is constructed under an assumption that respiration affects
62

G
PV
through the respiratory modulation effect on MAP. Further, the model assumes
linear and stationary relationship between MAP and G
PV
. The schematic of the model is
illustrated in Figure 3.5.

Figure 3.5 Schematic diagram of the one-input minimal model describing G
PV

fluctuations

The mathematical representation of the model is as follows:
( ) ( ) ( ) ( ) t T i t MAP i h t G
PV
M
i
BPC BPC PV σ
ε + − − ∆ ⋅ = ∆
∑
−
=
1
0

Equation 3.15
where h
BPC
represents the impulse response of the baroreflex control of peripheral vacular
conductance. h
BPC
quantifies the time-course of the changes in G
PV
resulting from an
abrupt increase in MAP of 1 mmHg. T
BPC
represents the latency associated with BPC
mechanism. This latency is assumed to be 3-8 beats (approximately 3-8 seconds) in
deBoer’s model (deBoer et al. 1987), 5 seconds in Madwed’s model (Madwed et al.
1989), and 2 seconds in van de Vooren’s model (van de Vooren et al. 2007). ε
G
PV

represents the dynamics of G
PV
that cannot be explained by the model.
Two-input Minimal Model
Unlike the one-input minimal model, this model assumes that respiration affects
G
PV
directly rather than through the modulation of MAP. Thus, it requires respiration as
h
BPC

Baroreflex Control of
Peripheral Vascular
Conductance
ε
G
PV

ΔG
PV
ΔMAP
63

the second input in addition to MAP as illustrated in Figure 3.6. The two-input minimal
model describes ΔG
PV
through two mechanisms: 1) BPC and 2) respiratory-peripheral
vascular conductance coupling (RPC). BPC relates ΔMAP to ΔG
PV
while RPC relates
changes in respiration (ΔILV) to ΔG
PV
. Again, both relationships are assumed to be linear
and stationary.

Figure 3.6 Schematic diagram of two-input minimal model describing G
PV
fluctuations
The mathematical representation of the model is as follows:

( ) ( ) ( )
( ) ( ) ( ) t T i t ILV i h
T i t MAP i h t G
PV
M
i
RPC RPC
M
i
BPC BPC PV
σ
ε + − − ∆ ⋅ +
− − ∆ ⋅ = ∆
∑
∑
−
=
−
=
1
0
1
0

Equation 3.16
where h
BPC
and h
RPC
represent the impulse responses of BPC and RPC. h
BPC
quantifies
the time-course of the changes in G
PV
resulting from an abrupt increase in MAP of 1
mmHg and h
RPC
quantifies the time-course of the changes in G
PV
following a very rapid
increase in ILV of 1 unit. T
BPC
and T
RPC
represent the latencies associated with BPC and
Respiratory-Peripheral
Vascular Conductance
Coupling
ΔILV h
RPC

Baroreflex Control of
Peripheral Vascular
Conductance
ε
G
PV

ΔG
PV
ΔMAP h
BPC

64

RPC mechanisms. The time delay in BPC mechanism is assumed to be the same as the
one-input model. For the latency in RPC dynamics, we observed relatively fast response
in our experimental data. Previous studies reported the time from the deep inspiration to
the peak response in vasoconstriction is 2-5 seconds (Bolton et al. 1936; Browse and
Hardwick 1969). In this study, the time delay in RPC, i.e. from the onset of the deep
inspiration to the onset of the vasoconstriction, is assumed to be 0-3 seconds. Lastly, ε
G
PV

represents the dynamics of G
PV
that cannot be explained by the model.
3.3.2 Nonlinear Model of Peripheral Vascular Conductance

Figure 3.7 Schematic diagram of the nonlinear model describing GPV fluctuations
ΔMAP
ε
G
PV

ΔG
PV

ΔILV
h
BPC

h
2BPC

h
BPC,RPC

h
2RPC

h
RPC

Respiratory-Peripheral
Vascular Conductance
Coupling
Baroreflex Control of
Peripheral Vascular
Conductance
Interaction Effect
of the 2 Inputs

65

While the linear models may be able to partially capture the dynamics of G
PV
,
they cannot provide information regarding nonlinear dynamics, which have been
commonly observed in various physiological systems, or the interaction between BPC
and RPC dynamics. In this study, we employ a two-input-one-output 2
nd
-order nonlinear
and stationary model to represent the relationships between ΔILV and ΔMAP to ΔG
PV
.
Figure 3.7 shows the configuration of the nonlinear model.
The mathematical representation of the model is as follows:
( ) ( ) ( )
( ) ( )
( ) ( ) ( )
( ) ( ) ( )
( ) ( ) ( )
( ) t
T j t ILV T i t MAP j i h
T j t ILV T i t ILV j i h
T j t MAP T i t MAP j i h
T i t ILV i h
T i t MAP i h t G
PV
RPC
M
i
M
j
BPC RPC BPC
RPC
M
i
M
j
RPC RPC
BPC
M
i
M
j
BPC BPC
M
i
RPC RPC
M
i
BPC BPC PV
σ
ε +
− − ∆ ⋅ − − ∆ ⋅ +
− − ∆ ⋅ − − ∆ ⋅ +
− − ∆ ⋅ − − ∆ ⋅ +
− − ∆ ⋅ +
− − ∆ ⋅ = ∆
∑ ∑
∑ ∑
∑ ∑
∑
∑
−
=
−
=
−
=
−
=
−
=
−
=
−
=
−
=
1
0
1
0
,
1
0
1
0
2
1
0
1
0
2
1
0
1
0
,
,
,

Equation 3.17
where h
BPC
and h
RPC
represent the linear dynamics of between the two system inputs and
the output; h
2BPC
and h
2RPC
represent the quadratic contributions of the system inputs;
h
BPC,RPC
is a cross-kernel which describes the 2
nd
-order interaction effect of the two inputs
on the output; T
BPC
and T
RPC
represent the latencies associated with BPC and RPC
mechanisms; and ε
G
PV
represents the extraneous influence on G
PV
that cannot be explained
by the model.
66

3.3.3 Model Estimation and Optimization
The minimal models proposed in this study are assumed to be stationary. The
characterization of both BPC and RPC dynamics can be accomplished by identifying
their corresponding system kernels relating the input(s) to the output of the system. In the
linear model case, the linear kernels are the impulse responses. Furthermore, due to an
inherent closed-loop structure of the model, it is necessary to impose causality constraints
in the model. This is achieved by imposing time delay in the model.
To estimate system kernels, we employ kernel expansion technique based on the
Volterra-Wiener theory of nonlinearity (Marmarelis 2004). This approach is chosen
because it can be applied on a relatively short and noise-contaminated data segment and
the system input requirement of being strictly white and Gaussian are relaxed(Marmarelis
1993). A system kernel can be represented as a weighted sum of a set of basis functions.
The Laguerre functions are often used as the orthonormal basis because of their built-in
exponential term, making them appropriate for modeling physiological systems
(Marmarelis 1993). In this study the Meixner basis functions (MBF), also known as
generalized Laguerre functions, are employed as they allow us to have control over the
rise time of the functions, making the approximation of kernels with slow start more
accurate (Asyali and Juusola 2005). Since each kernel can be expanded into a set of basis
functions, only the unknown expansion coefficients are required to be estimated, which,
greatly reduce the number of parameters to be estimated. As a consequence, an increase
in estimation accuracy can be obtained even when applied to relatively short data records
with the presence of noise (Marmarelis 1993).
67

Based on the nature of the dynamics of interest, the kernels are assumed to persist
for M sampling intervals where each sampling interval is 0.5 seconds. M represents the
memory of the system and is assumed to be 25 seconds (i.e. 50 samples). This implies
that the dynamics of BPC and RPC would become insignificant after 25 seconds. Using
this kernel expansion technique, each linear kernel such as in Equation 3.15-Equation
3.17 can be represented as

( ) ( ) ( )
∑
=
=
x
x
q
j
n
j x x
t B j c t h
1
) (
.
Equation 3.18
Each nonlinear self-kernel such as in Equation 3.17 can be represented as

( ) ( ) ( ) ( )
∑ ∑
= =
=
x x
x x
q
i
q
j
n
j
n
i xx xx
t B t B j i c t t h
1 1
2 1 2 1
) ( ) (
, , .
Equation 3.19
Lastly, the nonlinear cross-kernel in can be represented as

( ) ( ) ( ) ( )
∑ ∑
= =
=
x u
u x
q
i
q
j
n
j
n
i xu xu
t B t B j i c t t h
1 1
2 1 2 1
) ( ) (
, , .
Equation 3.20
h
x
represents the linear contribution of an input signal to the output. h
xx
represents the
quadratic contribution of an input signal to the output. h
xu
represents the effect on the
output the second-order nonlinear interaction of the two inputs. ) (
) (
t B
n
i
and ) (
) (
t B
n
j
are
the orthonormal sets of MBF with n
th
order of generalization. The larger the value of n,
the longer it takes for the MBF to reach its maximum value.
x
c correspond to the
expansion coefficients of the basis functions for the linear kernel.
xx
c and
xu
c
68

correspond to the expansion coefficients of the basis functions for the nonlinear self-
kernel and cross-kernel, respectively. q
x
and q
u
represent the Meixner function orders, the
total numbers of MBF, used in the expansion of the kernels.
Using the Equation 3.18-Equation 3.20, the three minimal models described in
section 3.3.1 and 3.3.2 can be rewritten, in order, as

( ) ( ) ( ) ( ) t t u j c t G
PV
BPC
j
q
j
BPC PV σ
ε + = ∆
∑
=1

Equation 3.21

( ) ( ) ( ) ( ) ( ) ( ) t t v j c t u j c t G
PV
RPC BPC
j
q
j
RPC j
q
j
BPC PV σ
ε + + = ∆
∑ ∑
= = 1 1

Equation 3.22

( ) ( ) ( ) ( ) ( ) ( )
( ) ( ) ( )
( ) ( ) ( )
( ) ( ) ( ) t v t u j i c
t v t v j i c
t u t u j i c
t t v j c t u j c t G
j i
q
i
q
j
RPC BPC
j i
q
i
q
j
RPC
j i
q
i
q
j
BPC
j
q
j
RPC j
q
j
BPC PV
BPC RPC
RPC RPC
BPC BPC
PV
RPC BPC
∑ ∑
∑ ∑
∑ ∑
∑ ∑
= =
= =
= =
= =
+
+
+
+ + = ∆
1 1
,
1 1
2
1 1
2
1 1
,
,
,
σ
ε

Equation 3.23
where u
j
(t) and v
j
(t) (or u
i
(t) and v
i
(t)) are the convolution of the basis functions ( ) t B
n
j
) (
(or ( ) t B
n
i
) (
) with the input ∆MAP and ∆ILV respectively:
69

( )
( )
( ) ( )
BPC
M
n
j j
T t MAP B t u − − ∆ =
∑
=
τ τ
τ 0

Equation 3.24
( )
( )
( ) ( )
RPC
M
n
j j
T t ILV B t v − − ∆ =
∑
=
τ τ
τ 0
.
Equation 3.25
The unknown coefficients c
BPC
, c
RPC
, c
2BPC
, c
2RPC
, and c
BPC,RPC
in Equation 3.21-
Equation 3.23 can now be estimated using least-squares minimization. This least-squares
minimization process was repeated for a range of values of delays (T
BPC
from 3.5 to 8
seconds and T
RPC
from 0 to 3 seconds), the order of generalization (n from 2 to 8), and
the Meixner function order (q
BPC
and q
RPC
from 2 to 6). For each combination of the
delays, order of generalization, and Meixner function order, the minimum description
length (MDL) (Rissanen 1982) is employed as a measure of the quality of data fitting and
can be computed as
( )
( )
M
M parameters of number total
J MDL
R
log
log
×
+ =
Equation 3.26
where
R
J is the variance of the residual errors between the observed and the predicted
output. The total number of parameters is the total number of the unknown coefficients to
be estimated. For the one-input and two-input linear models, the total numbers of
parameters are q
BPC
and q
BPC
+ q
RPC
, respectively. For the two-input 2
nd
-order model, the
total number of parameters is
70

( )
( )
.
RPC BPC
RPC RPC
BPC BPC
RPC BPC
q q
q q
q q
q q parameters of number total
⋅ +
+ ⋅
+
+ ⋅
+
+ =
2
1
2
1

Equation 3.27
The optimal set of parameters is selected based on the global search for minimal
MDL. It is worth noting that since M is fixed, MDL is mainly dependent on J
R
and the
total number of parameters. However, J
R
decreases with the higher Meixner function
order. Thus the global optimal set of parameters is the set with lowest Meixner function
order that yields the lowest variance of the residual errors. The unknown coefficients
c
BPC
, c
RPC
, c
2BPC
, c
2RPC
, and c
BPC,RPC
are computed using the optimal set of parameters.
The system kernels can then be estimated using Equation 3.18-Equation 3.20.
3.3.4 Physiological Interpretation of Nonlinear Dynamics
While the Volterra-Wiener approach to the system identification offers many
practical advantages, one of the limitations of this approach is how to interpret the
physiological significance of the estimated nonlinear kernels. In this section, we simulate
different scenarios that would aid us in understanding the underlying mechanisms
represented by the nonlinear kernels (Jo et al. 2007).
In most analyses of both the baroreflex control of HR and R
PV
, the baroreflex is
assumed to have linear gain over the pressure (input) range. However, many experiments
have shown that under the circumstances of the input becomes larger or smaller than the
normal operating range, the steady-state baroreflex response tends to reach saturation or
71

threshold level. To demonstrate that the 2
nd
-order BPC component reflects this
observation, representative h
BPC
and h
2BPC
are convolved with step functions of different
magnitude, ranging from 1 to 5 mmHg (Figure 3.8). We anticipate that the steady-state
response of the linear BPC component should increase linearly with each increase in the
step magnitude while the combined steady-state response of both linear and 2
nd
-order
BPC component should show saturation effect.

Figure 3.8 Illustration of how to obtain the simulated results from h
BPC
and h
2BPC

It has been previously observed that deep tidal breathing at normal breathing rate
enhances the modulatory effect on muscle sympathetic nerve activity (Seals et al. 1990).
This implies that the frequency response (gain) corresponding to higher V
T
would be
higher than that of lower V
T
. To demonstrate that the 2
nd
-order RPC component reflects
this observation, representative h
RPC
and h
2RPC
are convolved with a sinusoidal signal of
fixed frequency (0.5, 0.75, 1, 1.25, and 1.5 times the normal breathing frequency) and
Steady-state response
h
BPC
h
2BPC
Linear
Nonlinear
1 – 5 mmHg
Step functions with
different MAP magnitudes
72

fixed amplitude of V
T
(0.5, 1, 2, and 3 L) (Figure 3.9). The magnitude of the FFT of the
estimated G
PV
response at the sinusoidal input frequency is derived to represent the
frequency response of the RPC component. The frequency response is normalized to the
V
T
level. We anticipate that the normalized frequency response of the linear RPC
component would be the same for all V
T
levels; while the normalized frequency response
of the 2
nd
-order RPC component should be V
T
dependent.

Figure 3.9 Illustration of how to obtain the simulated results from h
RPC
and h
2RPC

Lastly, to demonstrate the physiological meaning of the 2
nd
-order cross-kernel
component, we simulate a case that mimics an experiment performed by Eckberg et al.
(1995). In that experiment, the neck suction (abrupt, brief pressure pulse to the neck) was
applied at different phases of the respiratory cycle, which resulted in larger increase in
sympathetic activity during expiration than during inspiration, showing the modulatory
effect of respiration on BPC. We thus convolve a representative h
BPC,RPC
with a semi-
sinusoidal tidal volume pattern and a positive MAP pulse signal triggered at different
time along the respiratory cycle (at an increment of 0.5 seconds). The semi-sinusoidal
Sinusoidal signal with different
V
T
amplitudes and frequencies
h
RPC
h
2RPC
y
RPC
y
2RPC
|Y
RPC
|
V
T
|Y
2RPC
|
V
T
Normalized
freq response
FFT
FFT
Linear
Nonlinear
73

waveform has amplitude of 1 L and a period (duration per breath) that is approximately
equal to the natural breathing frequency of the representative subject. The MAP pulse has
a pulse width of 0.5 seconds with a magnitude of 30 mmHg. We anticipate that the
resulting G
PV
following each MAP pulse would be dependent on the timing of the pulse
within the respiratory cycle.

Figure 3.10 Illustration of how to obtain the simulated results from h
BPC,RPC

3.3.5 Frequency Domain Representation of the Kernels
Another approach to further understand the physiological meaning of the
estimated kernels is to investigate these kernels from the frequency-domain perspective.
In brief, our goal is to determine a frequency response of a system such that we can
derive the gain of that system at the frequency range of interest. This section summarizes
the derivation of the frequency responses from the linear (1
st
-order) and nonlinear (2
nd
-
order) kernels as previously described by Marmarelis (Marmarelis et al. 1993; Chon et al.
V
T

Inspiration Expiration
MAP pulses simulating rapid BR stimulation
at different phases of breathing cycle
h
BPC,RPC
y
BPC,RPC
74

1994). Let’s assume that a system is exposed to an oscillatory input and the reason why
we assume this type of input is because, according to the Fourier analysis, any real
periodic signals (i.e. signals that occur in real, physical system) can be expressed as a
sum of analytical functions that take a form of
t j ω
e (where j is the imaginary number
( 1 − = j ), ω is the frequency in radian per time unit (ω = 2πf), and t is time). Therefore,
the most generalized form of the output can be obtained when such input is considered.
With these thoughts in mind, a frequency response of a system can be obtained simply by
taking the Fourier transform of the kernels representing that system. The magnitude of
the frequency response i.e. the amplitude of the output oscillations represents the system
gain. Note also that the system gain is dependent on the frequency.
Consider a system similar to the one described in section 3.3.2, the system output,
y(t), can be expressed as a sum of the output from each individual kernel:
( ) ( ) ( ) ( ) ( ) ( ) t y t y t y t y t y t y
xu uu xx u x
+ + + + =
Equation 3.28
where y
x
(t) and y
u
(t) represent the responses of the linear kernels whose inputs are x(t)
and u(t), respectively; y
xx
(t) and y
uu
(t) represent the responses of the second-order
nonlinear self-kernels whose inputs are x(t) and u(t), respectively; and y
xu
(t) represents
the response of the second-order nonlinear cross-kernel whose inputs are x(t) and u(t). Let
the two inputs be oscillatory inputs:
( ) ( )
t j t j
u x
b t u a t x
ω ω
e and e = =
Equation 3.29
75

where a and b represent the amplitudes of x(t) and u(t), respectively; and ω
x
and ω
u
are
the oscillation frequencies of x(t) and u(t), respectively. Then outputs of the linear
kernels, h
x
(t) and h
u
(t), are:

( ) ( )
( )
( )
( ) [ ] ( )
( ) ( )
t j
u u u
t j
x x x
t j
j
x
t j t j
x x
u
x x
x x x
bH t y
aH t h a
a h a a h t y
ω
ω ω
τ ω ω τ ω
ω
ω
τ τ τ τ
e
e e
d e e d e
0
=
= =
= =
−
∞
∞ −
−
∞
∫ ∫
F

Equation 3.30
where F represents the Fourier transform; and H
x
(ω) and H
u
(ω) represent the Fourier
transform of the linear kernels h
x
(t) and h
u
(t), respectively. From Equation 3.30, it means
that the output of these linear kernels would oscillate at the same frequency as their
respective inputs. The magnitudes of the frequency responses, denoted by |H
x
(ω)| and
|H
u
(ω)|, represent the linear gain of the system.
Similarly, the outputs of the second-order nonlinear kernels, h
xx
(t
1
,t
2
), h
uu
(t
1
,t
2
),
and h
xu
(t
1
,t
2
), to the oscillatory inputs, x(t) and u(t), are:

( ) ( )
( ) ( )
( )
( ) [ ]
( )
( ) ( )
( ) ( )
( )t j
u x xu xu
t j
u u uu uu
t j
x x xx
xx
t j
j j
xx
t j
t j t j
xx xx
u x
u
x
x
x x x
x x
abH t y
H b t y
H a
t t h a
h a
a a h t y
ω ω
ω
ω
ω
τ ω τ ω ω
τ ω τ ω
ω ω
ω ω
ω ω
τ τ τ τ
τ τ τ τ
+
∞
∞ −
− −
∞
∞ −
∞
− −
∞
=
=
=
=
=
=
∫ ∫
∫ ∫
e ,
e ,
e ,
, e
d d e e , e
d d e e ,
2 2
2 2
2 1
2 2
2 1 2 1
2 2
0
2 1
0
2 1
2 1
2 1
F

Equation 3.31
76

where H
xx
(ω,ω), H
uu
(ω,ω), and H
xu
(ω,ω) represent the 2D Fourier transform of the
second-order nonlinear kernels h
xx
(t
1
,t
2
), h
uu
(t
1
,t
2
), and h
xu
(t
1
,t
2
) respectively. From
Equation 3.31, the output of each nonlinear self-kernel will oscillate at twice the
frequency of its respective input. The output amplitude, given by the absolute value of the
Fourier transform of the kernel, also increases with the square of the input amplitude.
This means that output amplitude is dependent on the amplitude of the input i.e. the
stronger the input signal, the more pronounced the nonlinear behavior. The diagonal
elements of the Fourier transforms of the nonlinear self-kernels represent the nonlinear
response of a system at a given frequency. For the output of the nonlinear cross-kernel, it
will oscillate at the sum of the frequencies of its respective inputs; and the amplitude will
increase with the product of the two input amplitudes. The diagonal elements of the
Fourier transform of the cross-kernel provide the information about the effect of
interaction between the two inputs.
The derivation of the frequency response to the input containing multiple
frequencies can be achieved in a similar manner. Further details of the derivation can be
referred from Chon et al. (1994). In brief, for the nonlinear self-kernel case, there is an
additional term that appears in the output of the nonlinear self-kernel as a result of the
quadratic nonlinear interaction among the oscillations in each input. Similar to before, the
diagonal elements of the Fourier transform represent the nonlinear response at each of the
frequencies present in each input. If there is any interaction in the system, it would show
up in the off-diagonal elements of the nonlinear self-kernel Fourier transform. Lastly, the
off-diagonal elements of the cross-kernel Fourier transform represent various
77

combinations of the quadratic nonlinear interactions between the frequencies present in
the input signals.
3.3.6 Decorrelation of Model Inputs
In system identification, difficulty in identifying the model parameters can
potentially arise when there is collinearity. Collinearity occurs when the system inputs
(analogous to the predictor variables in regression analysis) are correlated to one another.
In physiological systems, it is known that blood pressure is modulated by respiration to a
certain degree (Toska and Eriksen 1993). We thus investigate whether decorrelation of
the system inputs would improve the estimation of the model parameters. To decouple
the correlation between ΔMAP and ΔILV, an autoregressive model with exogenous input
(ARX model) is used to model the coupling between the two signals with ΔILV being an
input and ΔMAP being an output (Chaicharn et al. 2009). The allowed ARX model order
ranges from 3 to 10 and the allow delay between ΔILV and ΔMAP ranges from 1.5 to 5
seconds. The optimal model order and delay are selected using MDL. The ARX model
prediction is the respiratory-correlated component of MAP, denoted as ΔMAP
rc
; while
the residuals of this ARX model is treated as the respiratory uncorrelated component of
MAP, denoted as ΔMAP
ru
. A 3-step process (Sangkatumvong 2011) is employed to
estimate the impulse responses of the 2-input linear model (see section 3.3.1). The
estimation of the impulse responses in each step undergoes the procedures described in
section 3.3.3.
78

Step 1: Estimate the temporary h
BPC
and h
RPC
(denoted as h
BPC,temp
and h
RPC,temp
,
respectively) from a 2-input-1-output model with ΔMAP
ru
and ΔILV being the inputs and
ΔG
PV
being the output (Figure 3.11).

Figure 3.11 Step 1 of the model estimation process using the decorrelated inputs
Step 2: Remove the blood pressure influence from ΔG
PV
such that it is only
influenced by respiratory effect (ΔG
PV,ILV
). Then estimate h
RPC
from a 1-input-1-output
model with ΔILV being an input and ΔG
PV,ILV
being an output (Figure 3.12).
( ) ( ) ( )
∑
−
=
− − ∆ ⋅ − ∆ = ∆
1
0
, , ,
M
i
temp BPC ru temp BPC PV ILV PV
T i t MAP i h t G σ
Equation 3.32

Figure 3.12 Step 2 of the model estimation process using the decorrelated inputs
Step 3: Remove the respiratory influence from ΔG
PV
to obtain the blood pressure
induced part of ΔG
PV
(ΔG
PV,MAP
). Then estimate h
BPC
from a 1-input-1-output model with
ΔMAP being an input and ΔG
PV,MAP
being an output (Figure 3.13).
ΔILV h
RPC,temp

ε
G
PV

ΔG
PV

ΔMAP
ru
h
BPC,temp

h
RPC

ε
G
PV

ΔG
PV,ILV
ΔILV
79

( ) ( ) ( )
∑
−
=
− − ∆ ⋅ − ∆ = ∆
1
0
,
M
i
RPC RPC PV MAP PV
T i t ILV i h t G σ
Equation 3.33

Figure 3.13 Step 3 of the model estimation process using the decorrelated inputs
The procedures described in this section are applied on the sickle cell dataset. The
normalized mean squared error (NMSE) between the observed data and the model
prediction of this method would be used to compare with the NMSE obtained from the
original method to determine whether the extra step of decorrelation of the inputs should
be incorporated in the final estimations of the BPC and RPC dynamics.
3.3.7 Residual Analysis
To evaluate the model performance, aside from being able to accurately estimate
the model components relating the inputs to the output, it is also essential to pay attention
to the residuals, the difference between the observed output and the predicted output,
representing the portion of the data that cannot be explained by the model. In this study,
the residual analysis consists of 3 parts: 1) the whiteness of the residuals, 2) the
correlation between the input and the output, and 3) the spectrum of the residuals. To test
the whiteness of the residuals, the normalized autocorrelation function of the residuals is
computed. If <5% of the autocorrelation function lies outside of the 95% confidence
interval, N 96 . 1 ± (where N is the number of samples), then the residuals are
h
BPC

ε
G
PV

ΔG
PV,MAP
ΔMAP
80

considered to be white. The whiteness of the residuals indicates that the fluctuations in
the output that are not due to random effects can be explained by the model. To test the
independence between the input and the output, the normalized cross-correlation function
is computed. Again, if the <5% of the cross-correlation function lies outside of the 95%
confidence interval, then the input and the output are consider to be independent of each
other. Otherwise, this indicates that the model cannot adequately describe how the input
relates to the part of the output. Lastly, the magnitude of the FFT of the residuals is
computed. By looking at the frequency components of the residuals, it can give us useful
information on possible mechanisms, based on our knowledge of the physiological
system under investigation, that are not sufficiently captured by the model and vice versa.
3.3.8 Estimation of Simulated Data
In this section, we investigate three aspects of the estimation. The first aspect is to
evaluate the accuracy of the estimated kernels of the 2-input linear model. The estimation
algorithm is applied on the “data” generated by the simulation model. The estimated
kernels are then compared with the “true” (simulated) kernels in the simulation model. To
test how the level of noise affects the accuracy of the estimation, the estimation algorithm
is applied on simulated data with different levels of system noise as well as measurement
noise. To explore the effect of the system noise, the coefficients of variation, controlling
beat-to-beat variability to SV, HR, and TPC, are increased from 0 to 0.1 (at an increment
of 0.02). We also investigate whether taking TPC, the output of the minimal model,
before and after adding the system noise would affect the estimation. Once simulated data
with different levels of system noise are obtained, different levels of the measurement
81

noise are added to the simulated TPC. The measurement noise is a white Gaussian noise
with zero mean. Its standard deviation is the same as the standard deviation of TPC
obtained from the default simulation model parameter configuration. The levels of the
measurement noise vary by changing the percentage of the standard deviation of the
measurement noise from 0 to 100% (at an increment of 25%) of the standard deviation.
Figure 3.14 shows the two configurations of how TPC are obtained from the simulation
model.

Figure 3.14 Configurations of how TPC, the output of the minimal model, was obtained
from the simulation model
The second aspect is to investigate how nonlinearity in the simulation model
affects the estimated kernels of the second-order nonlinear model. To induce more
nonlinearity to the simulated data, the linear operating range of the baroreflex control of
TPC
BPC

TPC
RPC

TPC before
system noise
TPC after
system noise
System Noise to TPC
Windkessel
Model
1
x
TPR
Measurement
Noise
Measurement
Noise
TPC, output of
Minimal model
TPC, output of
Minimal model
82

TPR is lowered. The estimation algorithm of the nonlinear model is then applied on the
simulated data from the original simulation model as well as the simulated data with the
modified baroreflex control of TPR. The estimated kernels, the nonlinear kernels in
particular, obtained from the original simulation model are then compared with the
estimated kernels from the simulation model with the modified baroreflex control of
TPR.
Lastly, we determine the relationship between the gains in the simulation model
and the system gains derived from the estimated linear kernels. This can be achieved by
deriving the system gains in both low-frequency (0.04-0.15 Hz) and high-frequency
(0.15-0.4 Hz) ranges from the estimated kernels. These derived system gains are then
compared with the respective gains in the simulation model.
83

Chapter 4. Results
This chapter presents the results of the analyses using the techniques described in
the previous chapter. The subject characteristics of both subject groups are first
displayed. The simulated data of the Simulation Model A under normal physiological
conditions as well as under interventions, are described. As previously mentioned,
vasoconstriction response is observed following a deep breath or a sigh. Simulation
Model A, which is developed based on previous literatures, shows that it cannot
reproduce such response to a sigh. We thus turn our attention to the respiratory effect on
peripheral vascular conductance (G
PV
= 1/R
PV
). The 1-input and the 2-input linear
minimal models are discussed in order to determine the effect of respiratory modulation
on G
PV
. The results suggest that respiration likely have direct modulatory effect on G
PV
.
This respiratory-coupling component was thus incorporated into the original simulation
model. The simulated data from this extended model are then displayed and compared
with the simulated data from the original model (Simulation Model A). The next section
shows the evaluation and validation results of the minimal modeling algorithm. This
includes the model prediction accuracy, the contribution of each estimated kernel to
describing the variability in G
PV
, as well as the ability of the minimal model to recover
the “true” characteristics of the system by feeding the “data” generated by the simulation
model to the minimal model, then compare the estimated impulse responses with the
“true” or the simulated impulse responses. Lastly, the results from applying minimal
84

modeling techniques on the collected data from both pediatric obese subjects and sickle
subjects are present.
4.1 Subject Characteristics
Obese Pediatrics Subjects
The subject characteristics of the selected obese pediatric subjects are displayed in
Table 4.1. Both groups were not statistically different in terms of age and BMI. Both
markers of the MetS + OSAS severity, OAHI and SI, in the metabolic group were
significantly higher than the control group. Figure 4.1 shows how the twenty selected
subjects clustered into two distinct groups in terms of the OSAS severity and insulin
sensitivity.
Table 4.1 Characteristics of the obese pediatric subjects in control, and metabolic
syndrome and obstructive sleep apnea syndrome groups
Ctrl (N = 10) MetS + OSAS (N = 10) P-Value
Gender (male/female) 7/3 10/0
Age (years) 12.4 ± 0.72 13.6 ± 0.59 0.140
BMI
†
(kg/m
2
) 30.4 ± 1.53 35.7 ± 2.41 0.078
OAHI
††
(events/hr) 1.56 ± 0.21 5.79 ± 1.15 <0.001
*

SI (×10
-4
min
-1
/ µU/ml) 6.67 ± 2.20 1.84 ± 0.17 0.021
*

Definition of abbreviations: Ctrl = control group; MetS + OSAS = metabolic syndrome and
obstructive sleep apnea syndrome group; BMI = body mass index; OAHI = obstructive apnea
hypopnea index; SI = insulin sensitivity.
Data show mean ± standard error (SE).
P-value from t-test (or Mann-Whitney rank sum test if failed normality test)
†
BMI ≥ 90
th
percentile for age and gender
††
Upper limit of normal in this age group = 1.5 events/hr
*
P < 0.05
85

Figure 4.1 Scatter plot of log(OAHI) and log(SI) of the control obese subjects (blue
circles) and metabolic syndrome and obstructive sleep apnea subjects (red crosses). The
dash lines show the median of log(OAHI) and the median of log(SI) from the total subject
pool.

Sickle Cell Subjects
The characteristics of the sickle cell subjects before and after blood transfusion
treatment are displayed in Table 4.2. After the blood transfusion treatment, the subjects
showed that the amount of hemoglobin in blood and the hematocrit (the volume
percentage of red blood cell in blood) significantly increased; while the reticulocyte count
(percentage of reticulocyte in the blood, giving a measure of how rapidly red blood cells
are produced and released by the bone marrow), the amount of sickle hemoglobin, and
HR significantly decreased. BMI, SBP, DBP, and SaO
2
before and after transfusion are
not statistically significantly different from each other.

-1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 3.5
-0.5
0
0.5
1
1.5
2
2.5
3
3.5
log(OAHI)
log(SI)
86

Table 4.2 Characteristics of the sickle cell subjects before and after transfusion
Pre-transfusion Post-transfusion P-Value
Gender (male/female) 7/7
Ethnicity (NonHis/His) 13/1
Age (years) 18.36 ± 1.56
Years on Transfusion 8.31 ± 1.33
BMI (kg/m
2
) 21.51 ± 1.08 21.48 ± 1.07 0.697
Hemoglobin (g/dl) 9.19 ± 0.20 11.59 ± 0.34 <0.001
*

Reticulocyte count (%) 12.45 ± 1.57 7.59 ± 1.09 <0.001
*

HCT (%) 27.14 ± 0.61 33.74 ± 0.99 <0.001
*

HbS (%) 36.36 ± 5.08 31.38 ± 4.09 0.014
*

HR (beats/min) 78.29 ± 1.97 69.71 ± 1.85 0.002
*

SBP (mmHg) 116.4 ± 2.63 115.1 ± 2.54 0.527
DBP (mmHg) 61.21 ± 2.16 60.21 ± 1.52 0.717
SaO
2
(%) 96 ± 0.01 97 ± 0.00 0.578
Definition of abbreviations: NonHis = non-Hispanic; His = Hispanic; BMI = body mass
index; HCT = hematocrit; HbS = sickle hemoglobin; HR = heart rate; SBP = systolic blood
pressure; DBP = diastolic blood pressure; SaO
2
= oxygen saturation by pulse oximetry.
Data show mean ± standard error (SE).
P-value from paired t-test (or Wilcoxon signed rank test if failed normality test)
*
P < 0.05

87

4.2 Simulation Model: Model A
Simulated data under normal physiological conditions (baseline) were generated
by the Simulation Model A (Figure 4.2, left panel). Table 4.3 shows the baseline values
of the key variables. RRI showed strong coupling with respiration through RSA
mechanism – inspiration led to shortening in RRI (i.e. increase in HR) and vice versa.
MAP oscillations followed the breathing frequency as respiration affected ABP through
the mechanical effect (direct effect) as well as through the RRI (indirect effect).
Table 4.3 Baseline values of the key variables generated by Simulation Model A
Variables Mean ± STD
ILV (L) 0.50 ± 0.35
HR (beats/min) 82.67 ± 10.64
SBP (mmHg) 117.46 ± 5.85
MAP (mmHg) 99.51 ± 5.11
DBP (mmHg) 83.97 ± 5.27
TPR (mmHg·s·mL
-1
) 1.00 ± 0.04
Definition of abbreviations: ILV = instantaneous lung volume; HR = heart rate; SBP =
systolic blood pressure; MAP = mean arterial pressure; DBP = diastolic blood pressure; TPC =
total peripheral resistance.

Lastly, TPC also showed coupling with respiration but this was likely due to the
modulation of respiration in ABP. The power spectra of the four signals (Figure 4.2, right
panel) all showed their peak frequency at 0.3 Hz, which was the breathing frequency.
88

However, the peak at the breathing frequency in MAP power spectrum was rather weak
compared to RRI and TPC. In addition, the low-frequency oscillations at around 0.08-
0.09 Hz could be observed in RRI, MAP as well as TPC power spectra.
4.2.1 Simulated Interventions
Intervention 1: An abrupt increase in ABP
An abrupt increase in ABP resulted in RRI lengthening in the following beat
through HR baroreflex mechanism. The lengthening of RRI then immediately caused a
drop in ABP as the Windkessel model was allowed more time for the blood pressure
pulse to decay. The drop in ABP led to the shortening of RRI, again through HR
baroreflex, which subsequently caused an immediate increase in MAP. These overshoot
and undershoot processes repeated for approximately 5 seconds before returning to the
baseline. TPC was not affected by the abrupt increase in ABP until 5 seconds later. The
abrupt increase in ABP caused TPC to rise, which immediately affected ABP, causing it
to increase. This increase in ABP was due to the drop in the time constant in the
Windkessel model, allowing shorter time for the blood pressure pulse to decay. The
oscillations in TPC followed the pattern of oscillations in MAP before returning to its
baseline approximately 10 seconds after the initiation of the abrupt increase in ABP. The
simulated data following this intervention is illustrated in Figure 4.3A.
Intervention 2: An abrupt vasoconstriction
An abrupt vasoconstriction immediately affected the time constant of the
Windkessel model, causing the ABP waveform slope (which is proportional to total
peripheral resistance) to decrease. This subsequently led to an increase in DBP and thus
89

an increase in ABP was observed. This increase in ABP then caused a lengthening of RRI
through HR baroreflex mechanism and the lengthening of RRI then caused blood
pressure to decrease. This oscillation in RRI repeated for approximately 5 seconds before
returning to its baseline. The increase in ABP caused by the abrupt vasoconstriction
triggered an increase in TPC approximately 5 seconds later through TPR baroreflex
responses. Again, the pattern of the changes in TPC after its increase followed the
fluctuations in ABP. TPC finally returned to its baseline after approximate 10 seconds
after the initiation of the abrupt vasoconstriction. The simulated data of an abrupt
vasoconstriction is illustrated in Figure 4.3B. Note that at time 151 second, the simulated
TPC, in fact, went down to 0.58 mL·mmHg
-1
·s
-1
. However, it was not shown in the figure
in order to allow the fluctuations in TPC following the induced vasoconstriction be more
visible.
Intervention 3: A sigh
The simulated data with a sigh is illustrated in Figure 4.3C. An increase in
inspiration caused by a sigh immediately resulted in the shortening of RRI (through the
RSA mechanism) and the small reduction in ABP (through the mechanical effect of
respiration). Shortly after, ABP increased due to the decrease in RRI (i.e. increase in HR)
as a result of an increase in cardiac output. An increase in ABP then caused a lengthening
in RRI through HR baroreflex as well as an increase in TPC through TPR baroreflex
approximately 5 seconds later. However, the Simulation Model A was not able to
reproduce the vasoconstriction response following a sigh as observed in the experimental
data.
90

Figure 4.2 Simulation Model A: simulated data under normal physiological conditions
0 0.1 0.2 0.3 0.4 0.5
0
1
2
P
ILV
(L
2
/Hz)
0 0.1 0.2 0.3 0.4 0.5
0
0.02
0.04
P
RRI
(s
2
/Hz)
0 0.1 0.2 0.3 0.4 0.5
0
100
200
P
MAP
(mmHg
2
/Hz)
0 0.1 0.2 0.3 0.4 0.5
0
1
2
x 10
-3 P
TPC
((mL.mmHg
-1
.s
-1
)
2
/Hz)
Frequency (Hz)
145 150 155 160 165 170 175 180
0
0.5
1
1.5
ILV (L)
145 150 155 160 165 170 175 180
0.5
1
RRI (s)
145 150 155 160 165 170 175 180
50
100
150
SBP, MAP, DBP (mmHg)
145 150 155 160 165 170 175 180
0.9
1
1.1
TPC (mL.mmHg
-1
.s
-1
)
Time (sec)
Beat-to-beat Time Series Power Spectra
91

Figure 4.3 Simulation Model A: simulated data with interventions: an abrupt increase in ABP (A), an abrupt vasoconstriction (B),
and a sigh (C)
145 150 155 160 165 170 175 180
0
1
2
3
ILV (L)
145 150 155 160 165 170 175 180
0.5
1
1.5
RRI (s)
145 150 155 160 165 170 175 180
50
100
150
SBP, MAP, DBP (mmHg)
145 150 155 160 165 170 175 180
0.8
1
TPC (mL.mmHg
-1
.s
-1
)
Time (sec)
145 150 155 160 165 170 175 180
0
1
2
3
ILV (L)
145 150 155 160 165 170 175 180
0.5
1
1.5
RRI (s)
145 150 155 160 165 170 175 180
50
100
150
SBP, MAP, DBP (mmHg)
145 150 155 160 165 170 175 180
0.8
1
TPC (mL.mmHg
-1
.s
-1
)
Time (sec)
145 150 155 160 165 170 175 180
0
1
2
3
ILV (L)
145 150 155 160 165 170 175 180
0.5
1
1.5
RRI (s)
145 150 155 160 165 170 175 180
50
100
150
SBP, MAP, DBP (mmHg)
145 150 155 160 165 170 175 180
0.8
1
TPC (mL.mmHg
-1
.s
-1
)
Time (sec)
A
An Abrupt Increase in ABP An Abrupt Vasoconstriction A Sigh
B C
92

4.3 Respiratory Effect on Peripheral Vascular Conductance
4.3.1 Indirect vs. Direct Respiratory Effect on G
PV

Fluctuation in MAP, ILV and PATampN, reflecting fluctuation in G
PV
, were fed
into the linear one-input and two-input minimal models in order to investigate the effect
of respiration on G
PV
. In this analysis, data from sickle cell subjects were utilized. We
would focus on how well each model was able to predict the output. Normalized mean
squared error (NMSE) was employed to determine the extent to which the model
prediction fits the observed data. Two-way repeated measures ANOVA showed that the
grouped NMSE of the two-input model was overall smaller than that of the one-input
model (p<0.001) as shown Figure 4.4.

Figure 4.4 Grouped bar chart of 1-input-model and 2-input-model NMSE (mean ± SE)
across transfusion treatment groups. * Denotes p<0.001 for Holm-Sidak post hoc test.
Overall, the two-input model was able to capture more dynamics of G
PV

compared to the one-input model. Furthermore, the breathing frequency peak in the
Treatment
Pre-transfusion Post-transfusion
NMSE
0.5
0.6
0.7
0.8
0.9
1.0
1.1
1-input Linear Model
2-input Linear Model
* *
93

power spectrum of the residuals of the second model was also greatly diminished (Figure
4.5). These results suggest that it is important to include the direct effect of respiration
along with baroreflex-mediated effects as the key mechanisms that determine natural
variations in G
PV
.

Figure 4.5 Representative power spectra of the ILV, MAP, PATamp, the one-input linear
model residuals (dash-dotted magenta line), and the two-input linear model residuals
(dashed green lines).
4.3.2 Respiratory Input Signal
Using the linear two-input minimal model, we investigated which of the 3
respiratory signals would best capture the dynamics of the respiratory-peripheral vascular
conductance coupling (RPC). ILV, breath-to-breath tidal volume (V
T
) and breath-to-
breath expiration time (T
E
) were the candidates of this investigation. This analysis was
applied on the sickle cell data. Again, NMSE was employed as a measure of how well the
fluctuations in the G
PV
could be explained by each of the respiratory input. Two-way
0 0.1 0.2 0.3 0.4
0
0.5
1
P
ILV
0 0.1 0.2 0.3 0.4
0
20
40
P
MAP
0 0.1 0.2 0.3 0.4
0
2
4
Frequency (Hz)
P
PATamp
& P
resid

PATamp
1-input Resid
2-input Resid
94

repeated measures ANOVA showed that there was a difference in grouped NMSE
derived from each of the respiratory inputs (p<0.001).

Figure 4.6 Grouped bar chart of NMSE of the two-input linear model using ILV, V
T
and
T
E
as the respiratory inputs (mean ± SE) across transfusion treatment groups. * Denotes
p<0.05 for Holm-Sidak post hoc test.
Figure 4.6 shows the effect of each of the respiratory inputs on the NMSE. The
NMSE of the model with ILV and V
T
as the respiratory input were significantly lower
than when T
E
was the input in both pre-transfusion (p<0.004 and p<0.002, respective)
and post-transfusion (p=0.001 and p<0.001, respectively). There was no significant
difference in NMSE between the model using ILV and the model using V
T
as the
respiratory input. These results suggest that the fluctuations in G
PV
are better described
by ILV and V
T
than T
E
. In this study, we decided to use ILV as the respiratory input for
the rest of our analyses. This is because ILV is a continuous signal and thus it should be
able to show a smoother modulatory effect of respiration on the G
PV
while still contains
the information about the magnitude of the tidal volume like V
T
.
Treatment
Pre-transfusion Post-transfusion
NMSE
0.5
0.6
0.7
0.8
0.9
1.0
ILV
Vt
Te *
*
*
*
95

4.3.3 Respiratory Coupling in Non-sigh vs. Sigh Data
Next, we investigated whether the RPC impulse response (h
RPC
) estimated from
data segments containing no sighs would be different from the h
RPC
estimated from data
segments containing sighs. The analysis employed the two-input model, which was
applied on data from 7 sickle cell disease subjects before receiving blood transfusion. For
each subject, the non-sigh segment was obtained from the baseline recording where there
was not any spontaneous sighs while the sigh segment was obtained from the first
coached sigh to the end of the recording. Figure 4.7 shows the average h
RPC
estimated
from the non-sigh (left) and the sigh segments (right).

Figure 4.7 Estimated h
RPC
from non-sigh and sigh segments of 7 sickle cell subjects before
receiving blood transfusion
Two-way repeated measures ANOVA with the two factors being 1) the type of
the data segments (non-sigh vs. sigh) and 2) the time course of the estimated impulse
response (1-10 second with 1-second interval) was applied to test whether the respiratory
coupling was different when there was no sigh. Note that only the first 10 seconds of the
0 5 10 15 20
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
x 10
-3
Time (sec)
h
RPC
0 5 10 15 20
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
x 10
-3
Time (sec)
h
RPC
Non-Sigh Sigh
96

impulse response were tested because this was where most of the dynamics occurred. The
statistical result showed that the estimated h
RPC
from non-sigh and sigh data segments
were not different from each other (p=0.344). However, dynamics of the h
RPC
estimated
from both types of segments changed over time as reflected by the significant result in
the time factor (p<0.001). Post-hoc test showed that the 1
st
second of the non-sigh h
RPC

was different from h
RPC
at 4-8 seconds and the 2
nd
second was different from 4-6 seconds
(p<0.05). For the h
RPC
of the sigh-segment, its 1
st
second was different from 4-5 seconds
(p<0.05). Figure 4.8 shows the first 10 seconds of the estimated h
RPC
from both non-sigh
and sigh segments.
Time (sec)
1 2 3 4 5 6 7 8 9 10
h
RPC
-0.002
-0.001
0.000
0.001
0.002
Non-sigh
Sigh
**
**
**
*
*
†
†

Figure 4.8 Grouped bar chart of the first ten seconds of the h
RPC
estimated from non-sigh
and sigh segments of the sickle cell subjects (mean ± SE).
* Denotes significant difference compared to non-sigh h
RPC
at 1 second (p<0.05).
** Denotes significant difference compared to non-sigh h
RPC
at 1 and 2 second (p<0.05).
†
Denotes significant difference compared to non-sigh h
RPC
at 1 second (p<0.05).

97

4.4 Extended Simulation Model: Model B
Based on the results in section 4.3, respiration likely has an effect on G
PV
directly
rather than only through the modulatory effect on blood pressure. Therefore, in this
section, we would incorporate the RPC component into the Simulation Model A. This
extended simulation model would be denoted as the Simulation Model B. The RPC
component to be incorporated into the simulation model was estimated from data
collected from the obese pediatric subjects during supine posture.
4.4.1 Incorporating RPC into the Simulation Model
The RPC estimated from the experimental data was incorporated into the
simulation model such that it would be more representative of the actual cardiovascular
system. In this study, we searched for the h
RPC
of the subject whose NMSE of the two-
input linear model decreased the most from its respective one-input linear model.
This representative h
RPC
is connected to a gain, K
RPC
, to produce the sigh-
vasoconstriction reflex block. The input to this block is ILV and the output is the change
in TPC (ΔTPC
RPC
). ILV is convolved with h
RPC
to generate ΔTPC
RPC
. This ΔTPC
RRC
is
added to the output of the TPR baroreflex control block (TPC
BPC
) to make up the total
TPC. Thus, TPC is now influenced directly not only by ABP but also by ILV. Figure 4.9
shows the schematic diagram of the simulation model after incorporating the sigh-
vasoconstriction mechanism. We also allow the possibility of the RPC mechanism to be
nonlinear by adding a switch block that would determine the magnitude of K
RPC

depending on the tidal volume. Currently, the threshold in the switch block is set such
98

that K
RPC
remains constant throughout the simulation regardless of the tidal volume i.e.
the relationship between the vasoconstriction response and the tidal volume is assumed to
be linear. Under normal physiological conditions, K
RPC
is assumed to be equal to 1.
Figure 4.10 shows simulated data (left panel) and their corresponding power
spectra (right panel) under normal physiological conditions by Simulation Model B.
Table 4.4 shows the baseline values of the key variables. Similar to the Simulation Model
A, this model demonstrated strong coupling between RRI and respiration. This was
verified by a distinct peak at the breathing frequency in the power spectrum of the RRI
(Figure 4.10, right panel). Another distinct oscillations in the RRI occurred at 0.1 Hz.
MAP oscillations, to a degree, followed the breathing frequency but the oscillations
around 0.1 Hz was much stronger as shown in its power spectrum. TPC also showed
coupling with respiration as well as the oscillations at 0.1 Hz. In comparison with the
Simulation Model A, the HR generated by this model was lower by approximately 10
beats/min (i.e. increase in RRI by 0.1 seconds). The average blood pressure values (SBP,
MAP, and DBP) in both models were comparable. The most obvious change in the
Simulation Model B was an increase in the baseline TPR (i.e. decrease in TPC). The
source of this increase was the result of the addition of the sigh-vasoconstriction
mechanism to the model. Also, the variability of the TPR generated by this model was
also slightly higher than the TPR generated by the Simulation Model A.

99

Table 4.4 Baseline values of the key variables generated by Simulation Model B
Variables Mean ± STD
ILV (L) 0.50 ± 0.36
HR (beats/min) 71.53 ± 11.18
SBP (mmHg) 119.33 ± 6.19
MAP (mmHg) 100.85 ± 5.02
DBP (mmHg) 85.19 ± 5.65
TPR (mmHg·s·mL
-1
) 1.17 ± 0.05
Definition of abbreviations: ILV = instantaneous lung volume; HR = heart rate; SBP =
systolic blood pressure; MAP = mean arterial pressure; DBP = diastolic blood pressure; TPC =
total peripheral resistance.

4.4.2 Simulated Interventions
Intervention 1: An abrupt increase in ABP
In general, both Simulation Model A and B produced similar response to an
abrupt increase in ABP. The RRI increased as a result of the abrupt rise in ABP thru HR
baroreflex, which immediately caused the ABP to drop as the Windkessel model was
allowed more time for the blood pressure pulse to decay. Approximately 5 seconds later,
the TPR baroreflex caused the TPC to increase as a response to the abrupt increase in
ABP. The slight difference between the responses generated by each simulation model
was mainly due to reduction in the overall TPC (i.e. increase in TPR). Since the TPR was
higher, the time constant in the Windkessel model also increased, thus decreased the rate
of decay in the pressure pulse. This means that, as a result of a lengthening in the RRI,
the ABP would not decrease as much as it would do in the Simulation Model A. This
100

then led the RRI to increase even further. Therefore, the increase in the RRI due to the
abrupt change in ABP appeared to last longer in the Simulation Model B compared to the
Simulation Model A. The simulated data following this intervention is illustrated in
Figure 4.11A.
Intervention 2: An abrupt vasoconstriction
Again, the responses to an abrupt vasoconstriction in both Simulation Model A
and B were overall similar. An abrupt increase in TPR led to an increase in the time
constant of the Windkessel model, resulting in an increase in DBP. Thus, an immediate
increase in ABP was observed. The HR baroreflex then caused the lengthening in the
RRI in the following beat due to the increase in ABP while the TPR baroreflex caused the
TPC to slightly increase approximately 5 seconds later. Furthermore, due to the overall
TPC in this model was lower than the Simulation Model A, the response to the abrupt
vasoconstriction appeared to last longer because of the lengthening in the RRI. The
simulated data of an abrupt vasoconstriction is illustrated in Figure 4.11B. Note that at
time 151 second, the simulated TPC drop, in fact, went down to 0.38 mL·mmHg
-1
·s
-1
.
However, it was not shown in the figure in order to allow the fluctuations in TPC
following the induced vasoconstriction be more visible.
Intervention 3: A sigh
Figure 4.11C illustrates the simulated data during a sigh generated by the
Simulation Model B. Similar to the response produced by the Simulation Model A, the
shortening of the RRI occurred when the sigh was initiated, which then caused ABP to
increase in the following beat. The main difference in the response to the sigh produced
101

by this model compared to the previous model was the vasoconstriction response.
Approximately 2 seconds after the sigh, the TPC began to decrease until it reached the
nadir approximately 4 seconds later. As a result of the decrease in TPC, both ABP and
RRI increased. The increase in the ABP was likely due to the feedforward effect of the
vasoconstriction to the blood pressure – decreased in the Windkessel model’s time
constant led to lower rate of decay in the pressure pulse, which caused an increase in
ABP. The increase in the RRI was likely the byproduct of the vasoconstriction –
vasoconstriction led to an increase in ABP, which then lengthened the RRI through HR
baroreflex mechanism. The increase in TPC that would have happened in response to the
increase in ABP (at 156 second) like in the Simulation Model A was not evident in the
Simulation Model B. This was likely because the vasodilation response due to the TPR
baroreflex was combined with the stronger sigh-vasoconstriction response. Thus, only the
vasoconstrictive response was observed following a sigh.
Figure 4.12 displays the response following a sigh in the experimental data and
the simulated data generated by the Simulation Model A and B. Clearly, without an
addition of the sigh-vasoconstriction mechanism (Simulation Model A), the reduction in
the TPC as observed in the experimental data could not be reproduced. In terms of the
cardiovascular variability, the power spectra of the simulated data from both Simulation
Model A and B were reasonably comparable to the power spectra of the experimental
data (Figure 4.13). To summarize, with the Simulation Model B operating in a closed-
loop fashion, it demonstrated the ability to simulate the patterns of cardiovascular
variability observed in experimental data, including the sigh-vasoconstriction response.
102

4.4.3 Continuous vs. Intermittent Sigh-Vasoconstriction Mechanism
We investigated whether the RPC mechanism operated at all time (continuous)
but became apparent when a sigh occurred or it was only triggered during deep
inspiration (intermittent). To simulate the second hypothesis, the RPC component in the
Simulation Model B was configured such that it would be activated when the tidal
volume of the current breath exceeded 2 liters. Otherwise, TPC would only be governed
by the baroreflex control of TPR. Figure 4.14 shows the representative experimental data
during a sigh and the simulated data during a sigh with different RPC states (continuous
vs. intermittent). At the initiation of the sigh in the intermittent RPC case, RRI began to
drop (i.e. increase in HR). This immediately caused ABP to increase. Approximately 2
seconds after the initiation of the sigh, the TPC began to decrease until it reached the
nadir approximately 4 seconds later. As a result of the decrease in TPC, ABP and RRI
both increased. Ten seconds after the initiation of the sigh, TPC increased as a result of
the previous increase in ABP through TPR baroreflex. The main differences between the
continuous RPC and intermittent RPC cases were 1) more pronounced overshoot in TPC
after the vasoconstriction and 2) the overall high average TPC, which was similar to the
TPC value of the Simulation Model A.

103

Figure 4.9 Simulation Model B: Simulation model after incorporating the sigh-vasoconstriction mechanism (lower right corner).
Denoted in blue are continuous variables, green are beat-to-beat variables, and red are key variables
HR to ΔHR
s
HR
ΔHR
s
HR to ΔHR
p
HR
ΔHR
p
Baroreflex Control of HR
ABP to HR
HR
ABP
dt
d
K
m
−
MER
p
p
p
F 2 π
s
1
K
(s) H
+
=
s
s
s
F 2 π
s
1
K
(s) H
+
=
SA node
D
s
RC
t
e SBP P(t)
−
⋅ =
Windkessel
K
ILV,p
-K
ILV,s
RSA
HR
intrinsic
ABP to TPR
TPR
ABP
D
pr
s T 1
K
(s) H
pr
pr
pr
+
=
α -adrenergic effector
Baroreflex Control of TPR
Zero-order Hold
ILV
Saturation
SV= SV
0
+
K
a
·DBP
n-1
+
Q
m
·(RRI
n-1
- RRI
0
)
x
1
ILV
z
-1
1-beat Delay
x
1
C
1
Frank-Starling Effect
Noise
TPC
ABP
TPC
ABP
ILV
ABP
wk
ΔHR
s
ΔHR
p
HR
RRI
RRI
SBP
R RRI
n-1
DBP
n-1
DBP
SV
PP
Vt
Sigh Reflex
h(t)
t
Vt
Sigh Reflex Gain
threshold
K
RPC
Sigh-vasoconstriction Reflex
Vt
ILV
TPC
RPC
TPC
BPC
TPR
Noise
HR
Noise
SV
x
1
104

Figure 4.10 Simulation Model B: simulated data under normal physiological conditions
0 0.1 0.2 0.3 0.4 0.5
0
1
2
P
ILV
(L
2
/Hz)
0 0.1 0.2 0.3 0.4 0.5
0
0.05
0.1
P
RRI
(s
2
/Hz)
0 0.1 0.2 0.3 0.4 0.5
0
100
200
P
MAP
(mmHg
2
/Hz)
0 0.1 0.2 0.3 0.4 0.5
0
1
2
x 10
-3 P
TPC
((mL.mmHg
-1
.s
-1
)
2
/Hz)
Frequency (Hz)
145 150 155 160 165 170 175 180
0
0.5
1
1.5
ILV (L)
145 150 155 160 165 170 175 180
0.5
1
1.5
RRI (s)
145 150 155 160 165 170 175 180
50
100
150
SBP, MAP, DBP (mmHg)
145 150 155 160 165 170 175 180
0.8
1
TPC (mL.mmHg
-1
.s
-1
)
Time (sec)
Beat-to-beat Time Series Power Spectra
105

Figure 4.11 Simulation Model B: simulated data with interventions: an abrupt increase in ABP (A), an abrupt vasoconstriction (B),
and a sigh (C)
145 150 155 160 165 170 175 180
0
1
2
3
ILV (L)
145 150 155 160 165 170 175 180
0.5
1
1.5
RRI (s)
145 150 155 160 165 170 175 180
50
100
150
SBP, MAP, DBP (mmHg)
145 150 155 160 165 170 175 180
0.6
0.8
1
TPC (mL.mmHg
-1
.s
-1
)
Time (sec)
145 150 155 160 165 170 175 180
0
1
2
3
ILV (L)
145 150 155 160 165 170 175 180
0.5
1
1.5
RRI (s)
145 150 155 160 165 170 175 180
50
100
150
SBP, MAP, DBP (mmHg)
145 150 155 160 165 170 175 180
0.6
0.8
1
TPC (mL.mmHg
-1
.s
-1
)
Time (sec)
145 150 155 160 165 170 175 180
0
1
2
3
ILV (L)
145 150 155 160 165 170 175 180
0.5
1
1.5
RRI (s)
145 150 155 160 165 170 175 180
50
100
150
SBP, MAP, DBP (mmHg)
145 150 155 160 165 170 175 180
0.6
0.8
1
TPC (mL.mmHg
-1
.s
-1
)
Time (sec)
An Abrupt Increase in ABP An Abrupt Vasoconstriction A Sigh
A B C
106

Figure 4.12 Data during a sigh from the experimental data (A), the Simulation Model A (B), and the Simulation Model B (C).
Note that the y-scale of the plot in the last row in panel A is different from panel B and C.
Experiment Simulation Model A Simulation Model B
A B C
0 20 40 60
0
2
ILV (L)
0 20 40 60
0.5
1
1.5
RRI (s)
0 20 40 60
80
100
120
MAP (mmHg)
0 20 40 60
0.5
1
PATamp (au)
Time (sec)
0 20 40 60
0
2
ILV (L)
0 20 40 60
0.5
1
1.5
RRI (s)
0 20 40 60
80
100
120
MAP (mmHg)
0 20 40 60
0.8
1
1.2
TPC (mL.mmHg
-1
.s
-1
)
Time (sec)
0 20 40 60
0
2
ILV (L)
0 20 40 60
0.5
1
1.5
RRI (s)
0 20 40 60
80
100
120
MAP (mmHg)
0 20 40 60
0.6
0.8
1
TPC (mL.mmHg
-1
.s
-1
)
Time (sec)
107

Figure 4.13 Normalized power spectra of data under normal physiological conditions from the experimental data (A), the
Simulation Model A (B), and the Simulation Model B (C)
0 0.1 0.2 0.3 0.4 0.5
0
0.5
1
P
ILV
0 0.1 0.2 0.3 0.4 0.5
0
0.5
1
P
RRI
0 0.1 0.2 0.3 0.4 0.5
0
0.5
1
P
MAP
0 0.1 0.2 0.3 0.4 0.5
0
0.5
1
P
PATamp
Frequency (Hz)
0 0.1 0.2 0.3 0.4 0.5
0
0.5
1
P
ILV
0 0.1 0.2 0.3 0.4 0.5
0
0.5
1
P
RRI
0 0.1 0.2 0.3 0.4 0.5
0
0.5
1
P
MAP
0 0.1 0.2 0.3 0.4 0.5
0
0.5
1
P
TPC
Frequency (Hz)
0 0.1 0.2 0.3 0.4 0.5
0
0.5
1
P
ILV
0 0.1 0.2 0.3 0.4 0.5
0
0.5
1
P
RRI
0 0.1 0.2 0.3 0.4 0.5
0
0.5
1
P
MAP
0 0.1 0.2 0.3 0.4 0.5
0
0.5
1
P
TPC
Frequency (Hz)
Experiment Simulation Model A Simulation Model B
A B C
108

Figure 4.14 Data during a sigh from the experimental data (A), the Simulation Model B with RPC operating continuously (B),
and the Simulation Model B with RPC operating intermittently (C).
Note that the y-scale of the plot in the last row in panel A is different from panel B and C.
0 20 40 60
0
1
2
3
ILV (L)
0 20 40 60
0.5
1
1.5
RRI (s)
0 20 40 60
80
100
120
MAP (mmHg)
0 20 40 60
0.4
0.6
0.8
1
1.2
PATamp (au)
Time (sec)
0 20 40 60
0
1
2
3
ILV (L)
0 20 40 60
0.5
1
1.5
RRI (s)
0 20 40 60
80
100
120
MAP (mmHg)
0 20 40 60
0.6
0.8
1
TPC (mL.mmHg
-1
.s
-1
)
Time (sec)
0 20 40 60
0
1
2
3
ILV (L)
0 20 40 60
0.5
1
1.5
RRI (s)
0 20 40 60
80
100
120
MAP (mmHg)
0 20 40 60
0.6
0.8
1
TPC (mL.mmHg
-1
.s
-1
)
Time (sec)
Experiment Continuous RPC Intermittent RPC
A B C
109

4.5 Minimal Model Evaluation and Validation
4.5.1 Model Prediction Accuracy
Obese Pediatric Subjects
The accuracy of each model was evaluated by computing the NMSE. The average
NMSEs of the linear and second-order nonlinear model applied on obese pediatric
subjects in supine and standing are displayed in Figure 4.15. Two-way repeated measures
ANOVA showed that the NMSE of the second-order nonlinear model was significantly
lower that the NMSE of the linear model (p<0.001). Post hoc test also showed this
reduction within each posture.

Figure 4.15 Grouped bar chart of linear-model and 2
nd
-order-nonlinear-model NMSE
(mean ± SE) across postures. * Denotes p<0.001 for Holm-Sidak post hoc test.
The residual analysis showed that the residuals of both models were not white and
there also existed the correlation between the residuals and the inputs. However, the
nonlinear model showed that it was able to capture the dynamics in the low frequency
and the breathing frequency ranges. Figure 4.16 displays the power spectra of the
Postures
Supine Stand
NMSE
0.0
0.2
0.4
0.6
0.8
1.0
1.2
Linear
2
nd
-order Nonlinear
*
*
110

PATamp signal, the residuals of the linear model, and the residuals of the nonlinear
model. The top row shows the power spectra of the subjects with the greatest difference
between the NMSE from the linear model and the NMSE of the nonlinear model in
supine and standing. The second row shows the power spectra of the subjects with
moderate NMSE reduction in supine and standing. Lastly, the third row shows the power
spectra of the subjects with the smallest NMSE reduction in both postures.

Figure 4.16 Representative power spectra of the PATamp signal (solid blue lines), the
linear model residuals (dashed green lines) and the second-order nonlinear model residuals
(dotted red lines), for supine and standing postures. From top to bottom: the spectra from
subjects with the largest, moderate and smallest NMSE reduction.
0 0.1 0.2 0.3 0.4
0
2
4
6
8
Supine
0 0.1 0.2 0.3 0.4
0
2
4
6
8
Stand
0 0.1 0.2 0.3 0.4
0
2
4
6
8
0 0.1 0.2 0.3 0.4
0
2
4
6
8
0 0.1 0.2 0.3 0.4
0
2
4
6
8
Frequency (Hz)
0 0.1 0.2 0.3 0.4
0
2
4
6
8
Frequency (Hz)
Largest NMSE
Reduction
Moderate NMSE
Reduction
Smallest NMSE
Reduction
111

Sickle Cell Subjects
Similar to the case of the obese pediatric subjects, there was a reduction in NMSE
when the second-order nonlinear model was applied compared to the linear model
(p<0.001). The post hoc test showed that the reduction in NMSE was significant in both
blood transfusion treatment groups (p<0.001). The average NMSE from both models and
both treatment groups are shown in Figure 4.17.

Figure 4.17 Grouped bar chart of linear-model and 2
nd
-order-nonlinear-model NMSE
(mean ± SE) across treatment groups. * Denotes p<0.001 for Holm-Sidak post hoc test.
The auto-correlation of the residuals of both model types showed that the
residuals were not white. The residuals and both inputs also still showed dependency on
each other. The power spectra of the residuals of both models showed that the nonlinear
model was better at capturing the dynamics of G
PV
at the low frequency range. Figure
4.18 shows the power spectra of the G
PV
signal, the residuals of the linear model, and the
residuals of the nonlinear model. The top row shows the power spectra of the subjects
with the greatest difference between the NMSE from the linear model and the NMSE of
Treatment
Pre-transfusion Post-transfusion
NMSE
0.0
0.2
0.4
0.6
0.8
1.0
Linear
2nd-order Nonlinear *
*
112

the nonlinear model in pre- and post-transfusion while the bottom row shows the power
spectra of the subjects with the smallest NMSE reduction in both treatment groups

Figure 4.18 Representative power spectra of the G
PV
signal (solid blue lines), the linear
model residuals (dashed green lines) and the second-order nonlinear model residuals
(dotted red lines), for pre- and post-transfusion. The top and the bottom rows: the spectra
from subjects with the largest and smallest NMSE reduction.
4.5.2 Decorrelation of the Model Inputs
The decorrelation of the model inputs technique described in section 3.3.6 was
applied on the sickle cell data to test whether it would improve the accuracy of the model
prediction. Figure 4.19 displays the average NMSE of the two-input linear model without
and with the input decorrelation. The NMSE of the model with the inputs being
decorrelated was slightly lower than the original model. Two-way repeated measures
ANOVA confirmed that the decorrelation of the input technique was not able to improve
0 0.1 0.2 0.3 0.4
0
1
2
3
4
x 10
-5
Pre-tranfusion
0 0.1 0.2 0.3 0.4
0
1
2
3
4
x 10
-5
Post-transfusion
0 0.1 0.2 0.3 0.4
0
0.2
0.4
0.6
0.8
1
x 10
-4
Frequency (Hz)
0 0.1 0.2 0.3 0.4
0
0.2
0.4
0.6
0.8
1
x 10
-4
Frequency (Hz)
Largest NMSE
Reduction
Smallest NMSE
Reduction
113

the performance of the model significantly in both blood transfusion treatment groups
(p=0.204).

Figure 4.19 Grouped bar chart of NMSE of the linear model and linear model with input
decorrelation (mean ± SE) across treatment groups.
4.5.3 Nonlinear Model Kernel Contributions to G
PV
Variability
For further insight into the mechanisms behind the variability in the G
PV
, we
investigated how much each kernel in the second-order nonlinear model contributed to
describing the variability of the model output. The ratio of the variance of model
prediction made by each kernel and the variance of the model output was employed as a
quantitative measure of the kernel contribution. Higher ratio means that the kernel could
describe more dynamics of the model output.
Figure 4.20 displays the contribution of each kernel on the variability of PATamp
collected from obese pediatric subjects. Among all the kernels in the model, the nonlinear
cross-kernel h
BPC,RPC
had the highest contribution to the variability in PATamp during
supine. Linear and nonlinear RPC, h
RPC
and h
2RPC
, provided the second highest
Treatment
Pre-transfusion Post-transfusion
NMSE
0.5
0.6
0.7
0.8
0.9
Original
Input Decorrelation
114

contribution. The nonlinear self-kernel h
2BPC
had the lowest contribution to the variability
in PATamp. In standing posture, the contribution mostly came from h
BPC,RPC
and h
2RPC
.
Both linear terms, h
BPC
and h
RPC
, had approximately equal contribution of the PATamp
variability. Again, the nonlinear self-kernel h
2BPC
had the lowest contribution to the
variability in PATamp.

Figure 4.20 Grouped bar chart of the ratio of the variance of the model prediction
by each kernel in the 2
nd
-order nonlinear model and the variance of PATamp (mean ± SE)
across postures (supine and standing).
* Denotes significant difference compared to full model (p<0.001).
** Denotes significant difference compared to full model and h
BPCRPC
(p<0.001).
Figure 4.21 displays the contribution of each kernel on the variability of G
PV

collected from sickle cell subjects before and after blood transfusion treatment. The
contributions from each kernel in both treatment groups were similar. Among all the
kernels in the model, the linear RPC term, h
RPC
, had the highest contribution to the
variability in G
PV
. The kernel with the second highest contribution was the cross-kernel
term, h
BPC,RPC
. The third highest contribution which contributed almost equally the same
Postures
Supine Stand
Prediction Variance / Data Variance
0.0
0.1
0.2
0.3
0.4
0.5
Full
h
BPC
h
RPC
h
2BPC
h
2RPC
h
BPCRPC
*
*
**
*
*
*
*
*
*
*
115

amount as the cross-kernel was the nonlinear self-kernel RPC term, h
2RPC
. The fourth
highest contribution was the nonlinear self-kernel BPC term, h
2BPC
. Lastly, the linear
BPC term, h
BPC
, had the lowest contribution to the fluctuations in G
PV
.

Figure 4.21 Grouped bar chart of the ratio of the variance of the model prediction
by each kernel in the 2
nd
-order nonlinear model and the variance of G
PV
(mean ± SE)
across treatment groups.
** Denotes significant difference compared to full model and h
RPC
(p<0.05).
4.5.4 Estimation of Simulated Data
Effect of System Noise and Measurement Noise
Figure 4.22 and Figure 4.23 show the estimated h
BPC
from the simulated data with
different levels of system noise and measurement noise obtained from the Simulation
Model B. When TPC was taken before adding the system noise (Figure 4.22), the
estimation algorithm was able to accurately estimate h
BPC
in terms of the delay of 5
seconds and the shape of the impulse response, with an exception to the zero system
noise case (the first column). With increasing level of measurement noise, the estimated
Treatment
Pre-transfusion Post-transfusion
Prediction Variance / Data Variance
0.0
0.1
0.2
0.3
0.4
0.5
0.6
Full
h
BPC

h
RPC

h
2BPC

h
2RPC
h
BPCRPC
**
**
**
**
116

h
BPC
, which should take a form of exponential decay, became more blunted. However, the
delay, the shape and the magnitude were still closely resembled the simulated (“true”)
impulse response. When TPC was taken after adding the system noise (Figure 4.23), the
estimated h
BPC
started from a negative value followed by an overshoot and an undershoot
within first 5 seconds. The following overshoot, however, moderately resembled the
shape of the simulated h
BPC
(thin black traces). As the level of the system noise increased,
the initial oscillation in the first 5 seconds became larger. At the highest system noise
level (CV = 0.1), the positive peak at 5 seconds became smaller than the peak of the
simulated impulse response overshoot at approximately 15 seconds. The addition of the
measurement noise had unnoticeable effect on the estimated h
BPC
.
Figure 4.24 and Figure 4.25 show the estimated h
RPC
from the simulated data with
different levels of system noise and measurement noise obtained from the Simulation
Model B. When TPC was taken before adding the system noise (Figure 4.24), the
estimated h
RPC
were recovered reasonably well in terms of capturing the negative
dynamics in TPC. The estimated delay of h
RPC
was 1 second longer than the delay in the
simulated impulse response. The estimation compensated this by estimating faster
decreasing rate. The drop in the estimated h
RPC
was also not as smooth as the simulated
h
RPC
and the magnitude of the negative peak was slightly larger. When the level of
system noise and the measurement noise were both high, the estimated h
RPC
became
smoother and the estimation of the delay was also correct. When TPC was taken after
adding the system noise (Figure 4.25), the estimated h
RPC
at CV = 0.02 – 0.04 (second
and third columns) showed some oscillations. For cases with higher system noise, the
117

estimated h
RPC
was smooth but the initial part of the estimation was not accurate as it
started with a positive peak. This was likely a compensatory effect for the initial negative
peak in h
BPC
. The addition of the measurement noise had little effect on the estimated
h
RPC
. For all cases where TPC was taken after the addition of system noise, the
estimation of the delay was 1 second longer than the simulated delay and the negative
peak magnitude was slightly larger than the peak of the simulated impulse response.
Effect of Continuous vs. Intermittent Sigh-Vasoconstriction
The estimated h
BPC
and h
RPC
from the simulated data with different operating RPC
states are shown in Figure 4.26. The estimated h
BPC
from the simulated data without RPC
and continuous RPC were almost identical to the simulated h
BPC
but the estimated delays
in both cases were 0.5 seconds longer than the simulated impulse response. The estimated
h
BPC
from the simulated data with intermittent RPC, however, showed an initial small
negative peak before becoming positive at 5 seconds. The positive peak magnitude was
slightly higher than the simulated h
BPC
. Following the positive peak was an undershoot
with brief oscillations. The dynamics of h
BPC
returned to zero after approximately 23
seconds. The estimated h
RPC
from the simulated data without RPC was small and
different from the simulated h
RPC
. It showed brief oscillations that lasted for
approximated 5 seconds. Its negative peak magnitude was around 5 times smaller than
the simulated h
RPC
. In the continuous RPC case, as expected, the estimated h
RPC
was
similar to the simulated h
RPC
in terms of both the time delay, the shape (negative
response) and the duration. However, its shape was not as smooth as the simulated
impulse response. The h
RPC
estimated from the simulated data with intermittent RPC
118

showed large negative response, approximately 2 times larger than the simulated h
RPC
.
However, the duration of the negative response was shorter than the simulated h
RPC
by
approximately 4 seconds. The negative response was followed by a small overshoot and
brief oscillations before returning to zero at 20 seconds.
Effect of Increased TPR Baroreflex Nonlinearity
The second-order nonlinear model estimation algorithm was applied on the
simulated data generated by the Simulation Model B. Both simulated data from the
original model and the model with increased TPR baroreflex nonlinearity were obtained
from the case where TPC was taken after the addition of system noise. Figure 4.27 shows
the estimated kernels from both simulated data sets. For the original model, the estimated
h
BPC
did not contain the initial negative part as the one estimated using the linear model
algorithm (Figure 4.23). The estimated h
RPC
was negative and smooth, similar to the
simulated h
RPC
but with larger magnitude. For the model with increased TPR baroreflex,
the estimated h
BPC
was different from the simulated one. It had the initial negative value
and the second overshoot did not resemble the shape of the simulated h
BPC
as in the
estimated h
BPC
using the linear model estimation algorithm (Figure 4.23, 2
nd
-6
th
columns).
The nonlinear BPC kernel, h
2BPC
, had small magnitude when estimated from the
original simulation model. However, the magnitude of h
2BPC
was much larger when it
was estimated from the simulation model with increased nonlinearity. The diagonal
elements of the h
2BPC
during the first 5 seconds of the estimated kernel was positive, then
followed by a small undershoot before returning to zero after 17 seconds. The nonlinear
RPC kernel, h
2RPC
, estimated from the model with increased nonlinearity in TPR
119

baroreflex control also displayed larger magnitude compared to the kernel estimated from
the original simulation model. The h
2RPC
estimated from the original simulation model
was positive and the reason for this was likely to compensate for the overestimated h
RPC
.
The h
2RPC
estimated from the simulation model with increased nonlinearity was negative.
Its initial negative part was followed by an overshoot, which was then followed by an
even larger undershoot. The dynamics returned to zero baseline after 17 seconds. Lastly,
the nonlinear cross-kernel, h
BPC,RPC
, estimated from the original simulation model was
almost symmetrical and negative. Like in the aforementioned estimated nonlinear
kernels, the h
BPC,RPC
estimated from the simulation model with increased nonlinearity had
larger magnitude. The h
BPC,RPC
was asymmetric with one side being more positive and the
other side being more negative.
Relationship between Simulated Gains and Estimated System Gains
Figure 4.28 shows how the system gains derived from the estimated kernels
changed with different levels of the simulated gains. Both low- and high-frequency BPC
gains were linearly correlated with the simulated TPR baroreflex gain (K
b,TPR
). As the
K
b,TPR
increased, BPC gains in both frequency ranges also increased. The low-frequency
BPC gain was more sensitive to the change the K
b,TPR
as reflected by steeper slope
compared to the high-frequency BPC gain. Similar to the BPC gains, the RPC gains in
both frequency ranges increased with increasing simulated gain of the sigh-
vasoconstriction reflex (K
RPC
). The low-frequency RPC gain showed linear correlation
with K
RPC
. However, the high-frequency RPC gains showed saturation effect as the
simulated gain became lower.
120

Figure 4.22 Estimated h
BPC
from the simulated data with different levels of system noise and measurement noise when TPC was
taken before adding system noise. Thin black tracing shows the simulated (“true”) h
BPC
.
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
Time (sec)
0 5 10 15 20 25
-2
0
2
x 10
-3
Time (sec)
0 5 10 15 20 25
-2
0
2
x 10
-3
Time (sec)
0 5 10 15 20 25
-2
0
2
x 10
-3
Time (sec)
0 5 10 15 20 25
-2
0
2
x 10
-3
Time (sec)
0 5 10 15 20 25
-2
0
2
x 10
-3
Time (sec)
Increasing System Noise: CV = 0 – 0.1
Increasing Measurement Noise: 0 – 100%
121

Figure 4.23 Estimated h
BPC
from the simulated data with different levels of system noise and measurement noise when TPC was
taken after adding system noise. Thin black tracing shows the simulated (“true”) h
BPC
.
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
0 5 10 15 20 25
-2
0
2
x 10
-3
Time (sec)
0 5 10 15 20 25
-2
0
2
x 10
-3
Time (sec)
0 5 10 15 20 25
-2
0
2
x 10
-3
Time (sec)
0 5 10 15 20 25
-2
0
2
x 10
-3
Time (sec)
0 5 10 15 20 25
-2
0
2
x 10
-3
Time (sec)
0 5 10 15 20 25
-2
0
2
x 10
-3
Time (sec)
Increasing System Noise: CV = 0 – 0.1
Increasing Measurement Noise: 0 – 100%
122

Figure 4.24 Estimated h
RPC
from the simulated data with different levels of system noise and measurement noise when TPC was
taken before adding system noise. Thin black tracing shows the simulated (“true”) h
RPC
.
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
Time (sec)
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
Time (sec)
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
Time (sec)
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
Time (sec)
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
Time (sec)
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
Time (sec)
Increasing System Noise: CV = 0 – 0.1
Increasing Measurement Noise: 0 – 100%
123

Figure 4.25 Estimated h
RPC
from the simulated data with different levels of system noise and measurement noise when TPC was
taken after adding system noise. Thin black tracing shows the simulated (“true”) h
RPC
.
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
Time (sec)
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
Time (sec)
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
Time (sec)
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
Time (sec)
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
Time (sec)
0 5 10 15 20 25
-0.06
-0.04
-0.02
0
0.02
Time (sec)
Increasing System Noise: CV = 0 – 0.1
Increasing Measurement Noise: 0 – 100%
124

Figure 4.26 Estimated h
BPC
and h
RPC
from simulated data without operating RPC (left), with continuously operating RPC (middle),
and RPC operating only with the presence of a sigh (right), in comparison to the simulated impulse responses (thin black lines)

0 10 20 30
-5
0
5
10
x 10
-4
0 10 20 30
-0.06
-0.04
-0.02
0
0.02
0 10 20 30
-5
0
5
10
x 10
-4
0 10 20 30
-0.06
-0.04
-0.02
0
0.02
0 10 20 30
-5
0
5
10
x 10
-4
0 10 20 30
-0.06
-0.04
-0.02
0
0.02
Without RPC Continuous RPC Intermittent RPC
h BPC h RPC
125

Figure 4.27 Estimated kernels (2
nd
-order nonlinear model) from the original Simulation
Model B “data” (left) and from “data” generated by Simulation Model B with increased
TPR baroreflex nonlinearity (right).
Original Increased TPR Baroreflex Nonlinearity
0 10 20 30
-0.04
-0.02
0
Time (sec)
h
RPC
0 10 20 30
-0.04
-0.02
0
Time (sec)
h
RPC
0 10 20 30
-3
-2
-1
0
1
x 10
-3

h
BPC
0 10 20 30
-4
-2
0
2
x 10
-3

h
BPC
126

Figure 4.28 Relationship between simulated gain and system gain derived from the estimated linear kernels at low- and high-
frequency ranges. Top row: TPR baroreflex gain (K
b,TPR
) vs. system gains of h
BPC
. Bottom row: sigh-vasoconstriction reflex gain
(K
RPC
) vs. system gains. of h
RPC
. Note that the ordinate scales for the low-frequency RPC gains are larger than the corresponding
scales for the high-frequency RPC gains.
0 0.5 1 1.5
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
x 10
-4
K
b,TPR
|H
BPC
|
LF
0 0.5 1 1.5
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
x 10
-4
K
b,TPR
|H
BPC
|
HF
0 0.5 1 1.5
0
0.2
0.4
0.6
0.8
1
x 10
-3
K
RPC
|H
RPC
|
HF
0 0.5 1 1.5
-2
0
2
4
6
8
x 10
-3
K
RPC
|H
RPC
|
LF
127

4.6 Effect of Metabolic Syndrome and OSAS on G
PV
Control
4.6.1 Estimated Kernels
Linear Kernels
The estimated linear BPC and RPC kernels from both linear and nonlinear models
are shown in Figure 4.29 and Figure 4.30, respectively. The estimated h
BPC
of the linear
model in both subject groups (control and MetS + OSAS) and both postures (supine and
standing) were larger than the h
BPC
estimated from the nonlinear model. They all showed
the initial negative peak followed by an overshoot at approximately 5-6 seconds. In
terms of postures, h
BPC
estimated from the linear model (h
BPC,Linear
) during supine had
larger positive peak compared to standing. During supine, the h
BPC,Linear
of the both
subject groups had comparable magnitude of the positive peak. However, during
standing, the positive peak of the h
BPC,Linear
of the MetS + OSAS group was smaller than
the control. For the h
BPC
estimated from the nonlinear model (h
BPC,Nonlinear
), the positive
peaks of the kernels in both subject groups were larger during supine compared to
standing. The positive peak of h
BPC,Nonlinear
in control was larger than the MetS + OSAS
during supine. However, this observation was reversed in standing posture.
Similar to h
BPC
, the estimated linear RPC kennels of the linear model (h
RPC,Linear
)
in subject groups and both postures were larger than the h
RPC
estimated from the
nonlinear model (h
RPC,Nonlinear
). All estimated h
RPC
were overall negative. The negative
peak of h
RPC
occurred between 0-2 seconds. The overall dynamics of h
RPC
were slower
than h
RPC
. The magnitude of the h
RPC,Linear
were comparable in both subject groups but
128

there were more oscillations in supine compared to standing . For the h
RPC,Nonlinear
during
supine posture, the negative response was followed by larger overshoot in both subject
groups. However, the h
RPC,Nonlinear
during standing did not have the overshoot following
the negative peak.
Second-order Nonlinear Kernels
The average second-order BPC kernels of both subject groups and postures are
shown in Figure 4.31. Overall, the dynamics of the h
2BPC
were slower than the linear BPC
kernel. h
2BPC
during supine in both subject groups started with negative peak while the
h
2BPC
during standing started with positive peak. In control group, the magnitude of the
h
2BPC
during supine was smaller than during standing. However, the dynamics of the
h
2BPC
during supine lasted longer compared to standing. Out of the four conditions, the
h
2BPC
of the MetS + OSAS group during supine had substantially larger magnitude.
Figure 4.32 shows the average second-order RPC kernels of both subject groups
and postures. All h
2RPC
showed negative response. The h
2RPC
of the MetS + OSAS group
during supine had the largest negative peak and the longest dynamics. In terms of
postures, the h
2RPC
of both subject groups during supine had larger magnitude than
standing.
The average second-order cross-kernels of both subject groups and postures are
shown in Figure 4.33. The h
BPC,RPC
showed a combination of both positive and negative
dynamics consistent with the dynamics found in the self-kernels. The magnitudes of the
h
BPC,RPC
were larger during supine compared to standing. Again, the MetS + OSAS group
during supine showed the largest kernel magnitude, in particular, the positive peak.
129

Figure 4.29 Average linear BPC kernels estimated from the linear (blue solid tracings) and
nonlinear model (red dashed tracings), corresponding to supine and standing postures in
control and MetS + OSAS groups.

Figure 4.30 Average linear RPC kernels estimated from the linear (blue solid tracings) and
nonlinear model (red dashed tracings), corresponding to supine and standing postures in
control and MetS + OSAS groups.
0 10 20 30
-0.02
0
0.02
h
BPC
Control
0 10 20 30
-0.02
0
0.02
h
BPC
Metabolic Syndrome
0 10 20 30
-0.02
0
0.02
h
BPC
Time (sec)
0 10 20 30
-0.02
0
0.02
h
BPC
Time (sec)
Supine Stand
MetS + OSAS
0 10 20 30
-0.4
-0.2
0
0.2
h
RPC
Control
0 10 20 30
-0.4
-0.2
0
0.2
h
RPC
Metabolic Syndrome
0 10 20 30
-0.4
-0.2
0
0.2
h
RPC
Time (sec)
0 10 20 30
-0.4
-0.2
0
0.2
h
RPC
Time (sec)
Supine Stand
MetS + OSAS
130

Figure 4.31 Average 2
nd
-order BPC kernels, corresponding to supine and standing postures
in control and MetS + OSAS groups.

Figure 4.32 Average 2
nd
-order RPC kernels, corresponding to supine and standing postures
in control and MetS + OSAS groups.
Supine Stand
MetS + OSAS
Supine Stand
MetS + OSAS
131

Figure 4.33 Average 2
nd
-order cross-kernels, corresponding to supine and standing
postures in control and MetS + OSAS groups.
4.6.2 Physiological Interpretation of Nonlinear Kernels
The simulation results for the physiological interpretation of the nonlinear kernels
are displayed in this section. Figure 4.34 shows the step-responses of the linear and the
combined linear and nonlinear components of BPC for different MAP step magnitudes.
The linear BPC step-responses increased linearly as the MAP step magnitude increased.
The combined linear and nonlinear BPC step-responses also increased with the increasing
MAP step magnitude but as the step magnitude became higher, the differential increase
became smaller.
Figure 4.35 shows the normalized frequency response of the linear and the
combined linear and nonlinear RPC components to sinusoidal signals with different
breathing frequencies and V
T
levels. The frequency responses of the linear RPC
Supine Stand
MetS + OSAS
132

component to different levels of the tidal volume all lined up perfectly on top of each
other, reflecting that the linear RPC responses were independent of the tidal volume
levels. However, the frequency responses of the combined linear and nonlinear RPC
showed dependency on the tidal volume levels – as the tidal volume increased, the
response became higher. Also, as the breathing frequency became lower than the natural
breathing frequency, the frequency response became higher. On the other hand, as the
breathing frequency became higher than the natural breathing frequency, the response
became smaller.
Figure 4.36 shows the simulation of the neck suction at different phases of the
breathing cycle (Eckberg 1995). The response to the MAP pulse was larger when the
MAP pulse was triggered during expiration. The peak response in the control subject
occurred just before midway through expiratory phase while the peak response in the
metabolic subject occurred right at the onset of the expiratory phase.

Figure 4.34 The linear (blue circles) and linear+nonlinear (red squares) BPC steady-state
response as a function of MAP step magnitude in control (left) and MetS + OSAS (right)
representative subjects.
1 2 3 4 5
1
2
3
4
5
6
7
x 10
-3
BPC response
MAP Step Amplitude (mmHg)
1 2 3 4 5
0.01
0.02
0.03
0.04
0.05
0.06
0.07
BPC response
MAP Step Amplitude (mmHg)
Control Metabolic Syndrome MetS + OSAS
133

Figure 4.35 The normalized frequency responses of the linear (top) and linear+nonlinear
(bottom) RPC component for different values of tidal volume and breathing frequency in
control (left) and MetS + OSAS (right) representative subjects.

Figure 4.36 Peak responses of the BPC-RPC cross-kernel to a 30 mmHg MAP pulse train
triggered at different times during the inspiratory/expiratory cycles in control (left) and
MetS + OSAS (right) representative subjects.

0.1 0.2 0.3 0.4 0.5
0
2
4
6
Linear RPC
0.1 0.2 0.3 0.4 0.5
-50
0
50
100
150
200
250
Linear + 2nd-order RPC
Breathing Frequency (Hz)
0.1 0.2 0.3 0.4 0.5
0
2
4
6
Linear RPC
0.1 0.2 0.3 0.4 0.5
-50
0
50
100
150
200
250
Linear + 2nd-order RPC
Breathing Frequency (Hz)
Control Metabolic Syndrome
0.5 - 3.0 L 0.5 - 3.0 L
0.5 L
1.0 L
2.0 L
3.0 L
0.5 L
1.0 L
2.0 L
3.0 L
MetS + OSAS
0 1 2 3 4
0
0.2
0.4
0.6
0.8
1
1.2
Vt
0 1 2 3 4
3
4
5
6
7
Peak Response
Time (sec)
Inspiration Expiration
0 1 2 3 4 5
0
0.2
0.4
0.6
0.8
1
1.2
Vt
0 1 2 3 4 5
3
4
5
6
7
Peak Response
Time (sec)
Inspiration Expiration
Control Metabolic Syndrome MetS + OSAS
134

4.6.3 Effect of Orthostatic Stress on System Gains
The overall effect of orthostatic stress on obese pediatric subjects (N = 49) caused
both linear and nonlinear gains to decrease. For the BPC gains derived from the linear
minimal model, paired t-test showed that both the low- and high-frequency BPC gains
decreased from supine to standing (p<0.001 and p<0.001, respectively). Similarly, the
RPC gains derived from the linear minimal model became smaller in standing compared
to supine (p=0.008 and p<0.001, respectively). As for the gains derived from the second-
order nonlinear minimal model, the gains also decreased due to orthostatic stress (Figure
4.37). The low- and high-frequency linear BPC and RPC gains showed significant
reduction from supine to standing (p<0.001 in all 4 cases). The low- and high-frequency
nonlinear BPC gain decreased from supine to standing (p<0.001 and p=0.022,
respectively). For the nonlinear RPC gains, only the high-frequency gain significantly
decreased from supine to standing (p<0.001). Lastly, both low- and high-frequency gains
of the cross-kernel showed significant reduction from supine to standing (p<0.001 and
p=0.001, respectively).

135

Figure 4.37 Nonlinear-model gains from supine to standing in both low- and high-
frequency ranges in obese pediatric subjects (N=49)

|H
BPC
|
Frequency Range
LF HF
0.000
0.001
0.002
0.003
0.004
0.005
|H
RPC
|
Frequency Range
LF HF
0.00
0.05
0.10
0.15
0.20
0.25
0.30
Supine
Stand
|H
2BPC
|
Frequency Range
LF HF
0
5e-6
1e-5
2e-5
2e-5
|H
2RPC
|
Frequency Range
LF HF
0.00
0.01
0.02
0.03
|H
BPC,RPC
|
Frequency Range
LF HF
0.0000
0.0004
0.0008
0.0012
** p < 0.001
* p = 0.001
†
p = 0.022
**
**
**
**
**
**
**
*
†
136

To investigate how the MetS + OSAS affected the autonomic reactivity to
orthostatic stress, the percentage changes in the BPC and RPC gains in both low- and
high-frequency ranges from supine to standing were computed. Figure 4.38 shows the
linear-model BPC and RPC gains. There was no statistical difference in the BPC or RPC
reactivity between the control and the MetS + OSAS groups. However, the gains in
supine postures were in general higher than the gains in standing. The nonlinear-model
BPC and RPC reactivity is shown in Figure 4.39. For the nonlinear-model gains, most of
them showed negative change from supine to standing. Mann-Whitney rank sum test
showed that the reduction in the gains of the 2
nd
-order BPC in the low-frequency range
from supine to standing in control subjects was smaller than in metabolic subjects
(p=0.031).

Figure 4.38 Linear-model gains from supine to standing in both low- and high-frequency
ranges (median and IQR) across subject groups.
Supine Stand
10
-4
10
-3
10
-2
Supine Stand
10
-4
10
-3
10
-2
Supine Stand
10
-5
10
-4
10
-3
10
-2
Supine Stand
10
-5
10
-4
10
-3
10
-2
Supine Stand
10
-2
10
-1
Supine Stand
10
-2
10
-1
Supine Stand
10
-2
10
-1
Supine Stand
10
-2
10
-1
|H
BPC
|
|H
RPC
|
Ctrl MetS Ctrl MetS
Low-Frequency Gains High-Frequency Gains
MetS + OSA MetS + OSAS
137

Figure 4.39 Nonlinear-model gains from supine to standing in both low- and high-
frequency ranges (median and IQR) across subject groups.

Supine Stand
10
-3
10
-2
Supine Stand
10
-3
10
-2
Supine Stand
10
-5
10
-4
10
-3
10
-2
Supine Stand
10
-5
10
-4
10
-3
10
-2
Supine Stand
10
-2
10
-1
Supine Stand
10
-2
10
-1
Supine Stand
10
-4
10
-3
10
-2
10
-1
Supine Stand
10
-4
10
-3
10
-2
10
-1
Supine Stand
10
-6
10
-5
Supine Stand
10
-6
10
-5
Supine Stand
10
-9
10
-8
10
-7
10
-6
10
-5
Supine Stand
10
-9
10
-8
10
-7
10
-6
10
-5
Supine Stand
10
-3
10
-2
10
-1
Supine Stand
10
-3
10
-2
10
-1
Supine Stand
10
-5
10
-4
10
-3
10
-2
Supine Stand
10
-5
10
-4
10
-3
10
-2
Supine Stand
10
-4
10
-3
Supine Stand
10
-4
10
-3
Supine Stand
10
-7
10
-6
10
-5
10
-4
Supine Stand
10
-7
10
-6
10
-5
10
-4
|H
BPC,RPC
|
|H
BPC
|
|H
RPC
|
|H
2BPC
|
|H
2RPC
|
Ctrl MetS Ctrl MetS
Low-Frequency Gains High-Frequency Gains
MetS + OSA MetS + OSAS
138

4.7 Effect of Blood Transfusion on G
PV
Control
4.7.1 Estimated Kernels
Linear Kernels
Figure 4.40 shows the average linear BPC kernels estimated from the linear and
second-order nonlinear model in pre- and post-transfusion treatment groups. In the pre-
transfusion group, the linear-model BPC kernel was larger than the nonlinear-BPC
kernel. Both of them started with an initial positive peak. The linear-model BPC kernel
showed more oscillations compared to the nonlinear-model one. In the post-transfusion
group, the nonlinear-model BPC kernel was substantially larger than the linear-model
kernel and started with a negative peak. However, the linear-model BPC kernel in post-
transfusion group started with a positive peak. The linear-model BPC kernel in pre-
transfusion group had larger positive peak compared to post-transfusion. On a contrary,
the positive peak of the nonlinear-model BPC kernel of pre-transfusion was smaller than
post-transfusion.

Figure 4.40 Average linear BPC kernels estimated from the linear (blue solid tracings) and
nonlinear model (red dashed tracings), corresponding to pre- and post-transfusion.
0 10 20 30
-2
-1
0
1
2
x 10
-5
h
BPC
Time (sec)
Pre-transfusion
0 10 20 30
-2
-1
0
1
2
x 10
-5
h
BPC
Time (sec)
Post-transfusion
139

Figure 4.41 displays the average linear RPC kernels estimated from the linear and
second-order nonlinear model in pre- and post-transfusion treatment groups. All
estimated linear RPC showed negative response. In the pre-transfusion group, both the
linear- and nonlinear-model RPC kernels were comparable in terms of both the
magnitude and the dynamics. In the post-transfusion group, the linear-model RPC kernel
had larger negative peak compared to the nonlinear-model one. Both of them also had
similar dynamics. Both linear- and nonlinear-model RPC kernels in pre-transfusion group
was smaller than their respective counterparts in the post-transfusion group. However, the
change in the kernel magnitude was larger between the linear-model kernels from both
treatment groups.

Figure 4.41 Average linear RPC kernels estimated from the linear (blue solid tracings) and
nonlinear model (red dashed tracings), corresponding to pre- and post-transfusion.
Second-order Nonlinear Kernels
The average second-order BPC kernels from both pre- and post-transfusion
groups are shown in Figure 4.42. The kernel magnitude of the pre-transfusion group was
evidently smaller compared to the post-transfusion group. The pre-transfusion h
2BPC

started with a small negative peak followed by an overshoot. The opposite occurred to the
0 10 20 30
-8
-6
-4
-2
0
x 10
-4
h
RPC
Time (sec)
Pre-transfusion
0 10 20 30
-8
-6
-4
-2
0
x 10
-4
h
RPC
Time (sec)
Post-transfusion
140

post-transfusion h
2BPC
– the kernel started with a positive peak followed by longer and
larger oscillations. The dynamics in the second-order BPC kernels in both groups
decayed to zero after approximately 20 seconds.
The average second-order RPC kernels from both pre- and post-transfusion
groups are shown in Figure 4.43. Both the pre- and post-transfusion kernels started with
negative peaks. However, the post-transfusion h
2RPC
has larger negative response. The
diagonal elements in the pre-transfusion h
2RPC
were positive but the off-diagonal elements
were negative. The opposite effect was observed in the post-transfusion h
2RPC
. The
diagonal elements were negative between 0 to 12 seconds while the off-diagonal
elements were positive.
The average second-order cross-kernels from both pre- and post-transfusion
groups are shown in Figure 4.44. The h
BPC,RPC
showed a combination of both positive and
negative dynamics consistent with the dynamics found in the self-kernels. The kernel of
the post-transfusion group had larger dynamics compared to the pre-transfusion group.
The dynamics corresponding to the RPC in the post-transfusion kernel lasted longer
compared to the pre-transfusion group while the duration of the dynamics corresponding
to the BPC in both treatment groups were comparable.

141

Figure 4.42 Average 2
nd
-order BPC kernels, corresponding to pre- and post-transfusion.

Figure 4.43 Average 2
nd
-order RPC kernels, corresponding to pre- and post-transfusion.

Figure 4.44 Average 2
nd
-order cross-kernels, corresponding to pre- and post-transfusion.
142

4.7.2 System Gains
To determine the effect of blood transfusion treatment on the control of G
PV
, the
percent change in the BPC and RPC gains from pre-transfusion to post-transfusion were
calculated. Figure 4.45 shows the linear-model BPC and RPC gains in both low- and
high-frequency ranges in both pre- and post-transfusion groups. Overall, the low-
frequency BPC and RPC gains increased from pre- to post-transfusion; while the high-
frequency BPC and RPC gains slightly decreased from pre- to post-transfusion. One-
sample t-test indicated that there was an increase in the low-frequency BPC gain from
pre- to post-transfusion but it was not able to reach the statistical significance (p=0.091).

Figure 4.45 Linear-model gains from pre-transfusion to post-transfusion in both low- and
high-frequency ranges (median and IQR).
Pre Post
10
-6
10
-5
Pre Post
10
-7
10
-6
10
-5
Pre Post
10
-5
10
-4
10
-3
Pre Post
10
-6
10
-5
10
-4
|H
BPC
|
|H
RPC
|
Low-Frequency Gains High-Frequency Gains
143

The nonlinear-model BPC and RPC gains from pre- to post-transfusion are
displayed in Figure 4.46. There was no statistical significant change from pre- to post-
transfusion in the gains derived from the linear components. One-sample signed rank test
showed that there was an increase in the second-order BPC gain in both the low- and
high-frequency ranges from pre- to post-transfusion (p=0.017 and p=0.009, respectively).
The second-order RPC gain in the high-frequency range also showed an increase from
pre- to post-transfusion but the change did not achieve the significant level (p=0.078).
The high-frequency gain derived from the second-order cross-kernel also increased from
pre- to post-transfusion (p=0.009). The low-frequency gain of the cross-kernel also
increased after transfusion but the increase did not achieve the significant level
(p=0.078).

144

Figure 4.46 Nonlinear-model gains from pre-transfusion to post-transfusion in both low-
and high-frequency ranges (median and IQR).
Pre Post
10
-6
10
-5
10
-4
Pre Post
10
-8
10
-7
10
-6
10
-5
Pre Post
10
-6
10
-5
10
-4
10
-3
Pre Post
10
-7
10
-6
10
-5
10
-4
Pre Post
10
-9
10
-8
10
-7
10
-6
Pre Post
10
-11
10
-10
10
-9
10
-8
10
-7
Pre Post
10
-7
10
-6
10
-5
10
-4
Pre Post
10
-9
10
-8
10
-7
10
-6
Pre Post
10
-7
10
-6
10
-5
Pre Post
10
-9
10
-8
10
-7
|H
BPC,RPC
|
|H
BPC
|
|H
RPC
|
|H
2BPC
|
|H
2RPC
|
Low-Frequency Gains High-Frequency Gains
145

Chapter 5. Discussion
Using the minimal modeling approach to quantify the sigh-vasoconstriction
response observed in the experimental data, we were able to identify the respiratory
coupling on G
PV
component that existed not only with but also without the presence of a
sigh. We demonstrated that the simulation model of cardiovascular variability developed
based on previously published models could not simulate the vasoconstriction response
following a sigh. This was possible only when the respiratory-G
PV
coupling mechanism,
which was estimated using the minimal model, was incorporated into the simulation
model. We also showed, using the simulation model, that taking the respiratory
modulation effect on G
PV
into account was essential for accurate estimation of the
dynamics of the baroreflex control of TPR.
We investigated the possibility of nonlinear dynamics in both the baroreflex
control, BPC, and the respiratory coupling, RPC, mechanisms and found that both
mechanisms exhibited nonlinear behavior. Including the nonlinear dynamics as well as
the interaction in addition to just the linear dynamics significantly decreased the model
prediction error. Most of the contributions to the variability in G
PV
originated from the
second-order nonlinear respiratory coupling term and the interaction effect between ABP
and respiration. We found that in obese pediatric subjects, the system gains derived from
the estimated BPC and RPC kernels decreased from supine to standing. The reduction in
the gains suggests that, the sympathetic modulation of G
PV
decreased, even though it is
well known that orthostatic stress increases sympathetic tone. There was a greater
146

reduction in the BPC nonlinear behavior in subjects with more severe degrees of
metabolic syndrome and obstructive sleep apnea than the non-severe ones. In the sickle
cell group, the nonlinear behavior in BPC as well as the interaction effect between ABP
and respiration increased after blood transfusion treatment.
To validate these findings, we tested the minimal model on the simulated “data”
under various simulation gains and increased nonlinearity. We found that the minimal
model was able to detect the increased nonlinear behavior generated in the simulation
model. Furthermore, the gains derived from the kernels estimated from the simulated data
showed correlation with the gains in the simulation model. By employing structured and
minimal modeling approaches in tandem, we developed an extended model of blood
pressure variability that reproduces the respiratory modulation effect on G
PV.
We have
also demonstrated that it is necessary to take this respiratory modulation effect into
account in order to achieve accurate TPR baroreflex estimation. The results also suggest
that the system gains derived from the estimated kernels could be employed as potentially
useful biomarkers of sympathetic nervous system functions. In the last section of this
chapter, we will discuss the limitations of the study.
5.1 Respiratory Effect on Peripheral Vascular Conductance
To date, numerous studies have extensively studied respiratory modulation on
heart rate but not many have investigated the respiratory modulation of R
PV
. From the
experimental data, we observed modulation of respiration on amplitude of PAT signal as
well as PU signal, which were employed as the surrogate measure of G
PV
, the inverse of
R
PV
. This observation was confirmed by the power spectra of both surrogate G
PV
signals.
147

Their spectra showed not only a peak in the low-frequency range, but also a peak at the
breathing frequency. Furthermore, previous studies have shown that muscle sympathetic
nerve activity fluctuates with a strong respiratory periodicity (Eckberg 1995).
Considering these observations, we turned our attention to further explore the effect of
respiration on G
PV
. Two hypotheses were proposed: 1) respiration affects G
PV
through its
modulation on ABP and 2) respiration modulates G
PV
directly. Using the minimal
modeling approach to characterize the relationship between ABP and G
PV
, we found that
most of the dynamics in G
PV
were not captured by the one-input model (first assumption
that respiration affects G
PV
through its modulation on ABP) as reflected in high
prediction error. Power spectral analysis was performed on the model residuals and we
found that power spectra of the residuals in a number of cases contained a visible peak at
the breathing frequency. This suggests that the respiratory component observed in G
PV

was not entirely governed by respiratory modulation on ABP. Based on the second
hypothesis, we developed another model that, in addition to ABP, also allowed a direct
effect of respiration on G
PV
(two-input model). We found that the prediction error was
significantly reduced and the breathing frequency peak in the power spectrum of the
residuals of the two-input model was also visibly attenuated or disappeared. These
findings support our hypothesis that respiration modulates G
PV
directly.
In addition to the modulatory effect of respiration, we also observed the sigh-
vasoconstriction response in our experimental data, consistent with what was also
previously reported by other studies (Bolton et al. 1936; Browse and Hardwick 1969;
Baron et al. 1996). This sigh-vasoconstriction response exists not only in adults but also
148

in infants (Inwald et al. 1996; Galland et al. 2000). Moreover, the vasoconstriction
follows not only a spontaneous sigh but also a voluntary deep breath as reported by
(Bolton et al. 1936) and by us in the coached sighs experiments . Though most of these
studies described the sigh-vasoconstriction as a reflex that was triggered only with the
presence of a sigh, we speculated that it could be a mechanism that operated
continuously; however, its response only became noticeable when the drive was large
enough such as when taking a deep breath. Both of the estimated respiratory coupling
components, h
RPC
, from the non-sigh and coached sigh segments showed negative
response with reasonably comparable in magnitude and dynamics (Figure 4.7). This
reveals that there is a mechanism similar to the sigh-vasoconstriction mechanism even
without the presence of a sigh.
Since the effect of respiratory modulation on G
PV
only became visible when there
was a sigh, thus earning the name “sigh-vasoconstriction reflex”, we investigated which
feature of respiratory signal would best capture this modulatory effect. A sigh is
characterized by large inspiratory volume and long expiration time (due to large
inspiratory volume). Three respiratory signals: instantaneous lung volume (ILV), breath-
to-breath tidal volume (V
T
) and breath-to-breath expiration time (T
E
) were investigated.
We found that fluctuations in G
PV
were better described by ILV and V
T
compared to T
E
.
This suggests that the lung volume might play a more significant role in modulating G
PV
.
Another study by Baron et al. (1996) also found that the respiratory modulation on blood
flow could be detected when there was high tidal volume and slow respiratory rate. Other
studies investigated the modulatory effect of respiration on peroneal muscle sympathetic
149

nerve activities (MSNA). Seals et al. (1990) reported that there was an increase in the
modulatory effect on MSNA during deep, slow frequency breathing but sustained
changes in breathing depth or pattern did not change the frequency of MSNA discharge.
The following study (Seals et al. 1993) confirmed the findings of the previous study that
V
T
and breathing frequency affected the within-breath variation of MSNA. Additionally,
they found that the potentiation of the within-breath pattern of the MNSA outflow at
increased tidal volume was not related to changes intravascular or intrathoracic pressures
due to baroreflex, or influenced by marked variations in the respiratory motor drive; but it
was instead dependent on the pulmonary vagal lung inflation feedback (Seals et al. 1993).
Although the respiratory effect was reported to be more correlated with V
T
, studies
showed that there was higher MSNA discharge during expiration (Eckberg 1995). This
led us to speculate that T
E
could be an appropriate respiratory signal. However, we found
that T
E
had the worst performance out of the three respiratory signals being tested in
terms of capturing G
PV
variability. Perhaps the fact that we found the increas in MSNA
firing during expiration was because of lower lung volume during expiration. Seals et al.
(1993) confirmed that at rest, approximately 70% of the MSNA discharge occurred
during the lower volume phase of the respiratory cycle. Therefore, we concluded that
ILV should be employed as a respiratory input for the minimal model since it contained
information about the tidal volume, the breathing frequency, as well as the inspiration
and expiration durations.

150

5.2 Simulation Models
In this study, we presented a closed-loop model of short-term cardiovascular
variability by employing a combination of structured and minimal modeling approaches.
The original simulation model (Simulation Model A) was able to generate the variability
of HR, ABP and TPR that were comparable to the actual measurements under normal
physiological conditions. The variability in HR was governed by the baroreflex control of
heart rate, the sinoatrial node, and the respiratory sinus arrhythmia mechanisms. It was
evident from the simulated data that the HR was tightly coupled with respiration. The
power spectrum of HR displayed dominant peaks at low-frequency range around 0.1 Hz
and at breathing frequency around 0.3 Hz as observed in the spectra from the
experimental data (Figure 4.13). The variability in ABP came from the baroreflexes
control of HR and TPR, the mechanical effect of respiration, as well as the influences
from the preload and afterload that affected PP. We demonstrated that the power
spectrum of MAP had a dominant peak at approximately 0.1 Hz. There was also
respiratory modulation in MAP but the strength of the modulation at the breathing
frequency was much weaker than the 0.1 Hz component. In the Simulation Model A, the
variation in TPR arose from the baroreflex control of TPR. Similar to the power spectra
of HR and ABP, the TPR spectrum showed dominant peak at 0.1 Hz and another smaller
but evident peak at the breathing frequency. It is known that there are two main debating
theories behind the origin of the 0.1-Hz oscillations. The first one was that these
oscillations are central in origin while the second theory suggests that they reflect the
resonance due to the baroreflex feedback. Since our simulation models (both A and B)
151

were able to produce this 0.1-Hz rhythm without the need of an additional oscillator, this
further supported the second hypothesis. Similar findings were also found in previous
studies (deBoer et al. 1987; Julien 2006; van de Vooren et al. 2007). However, this
finding still does not rule out the possibility that such oscillations could be intrinsic and
originate from the central nervous system.
Despite the fact that this developed simulation model (Simulation Model A) was
able to reproduce expected responses in a physiological system, it failed to reproduce the
vasoconstrictive response following a deep inspiration or a sigh, the kind of response
which had been observed in previous studies (Bolton et al. 1936; Browse and Hardwick
1969; Aso et al. 1997; Sangkatumvong et al. 2011) as well as in our experimental data.
Using the minimal modeling approach, the RPC component was identified. We
demonstrated that only after this component was incorporated into the simulation model
(denoted as the Simulation Model B) and the simulation model was allowed to operate in
a closed-loop manner that the sigh-vasoconstriction response could be accurately
predicted. To our knowledge, this is the first simulation model that has incorporated the
direct modulatory effect of respiration on TPR and is able to reproduce the sigh-
vasoconstriction response.
One should note that the simulation models developed in this study omit other
regulatory influences that may contribute to short-term cardiovascular variability such as
cardiopulmonary reflexes, chemoreflexes, and local vascular factors. The simulation
models also assumed respiratory input as an external input to the system, i.e. independent
of the cardiovascular effect, when in reality both the respiratory and cardiovascular
152

systems are integrated and closely interact with each other (Cheng et al. 2010). However,
in spite of this model’s highly simplified structure compared to the much more complex
physiological system, it was able reproduce expected responses as observed in a
physiological system, in particular the vasoconstriction response following a sigh. The
utility of the Simulation Model B is that it provides us with an improved understanding of
the underlying mechanisms that lead to oscillations at different frequencies. It can
potentially simulate cardiovascular variability under different conditions by adjusting the
appropriate simulation model parameters accordingly.
5.3 Minimal Models
The spontaneous variabilities observed in measurements of cardiovascular
parameters such as HR and ABP reflect the dynamic behavior of regulatory mechanisms
as they interact with perturbations that originate from external sources or other organ
systems (Parati et al. 2006). In this study, we aimed to characterize the mechanisms that
introduce variability in G
PV
. The fluctuations in G
PV
are generally assumed to be
governed autonomically by the baroreflex control of total peripheral resistance.
Combined with fluctuations in heart rate and stroke volume, these fluctuations in G
PV

generate blood pressure variability. However, we showed earlier that ABP alone cannot
sufficiently describe G
PV
dynamics. Respiratory modulation also plays a prominent role
in producing the G
PV
fluctuations. With the RPC as the additional mechanism, we
demonstrated that the “extended” minimal model with two inputs significantly improved
the ability of the model to capture the dynamics in G
PV
compared to the one-input model.
In addition, by having respiration as a separate arm of the dynamics, it enables us to
153

dissociate the confounding effects of respiration from blood pressure that contribute to
the variability in G
PV
. This point is demonstrated in Figure 5.1. The impulse response of
BPC, h
BPC
, was estimated from the data generated by the Simulation Model B. Note that
this simulation model has two mechanisms governing the changes in TPR: the baroreflex
control and the respiratory modulation. The more accurate estimation of h
BPC
could be
achieved only when the two-input model was applied, otherwise, the estimated h
BPC
also
captured the negative response, which was actually the response due to the respiratory
modulation on TPR.

Figure 5.1 Estimated h
BPC
from data generated by the Simulation Model B. Solid blue
trace shows h
BPC
estimated using MAP as the only input to the 1-input minimal model.
Dashed red trace shows h
BPC
estimated using MAP and ILV as the inputs to the 2-input
minimal model. Thin black trace shows the simulated (“true”) h
BPC
.

In this study, our focus was to examine the mechanisms responsible for short-term
variability of G
PV
. For this reason, the minimal model was assumed to be stationary
0 5 10 15 20 25
-1.5
-1
-0.5
0
0.5
1
1.5
x 10
-3
Time (sec)

1-input Model
2-input Model
Simulated
154

despite the fact that the actual physiological systems are commonly known to be non-
stationary. According to the guidelines about heart rate variability analysis, the Task
Force of The European Society of Cardiology and The North American Society of Pacing
and Electrophysiology (1996) suggested that, for a short-term recording, a 5-minute
segment was preferred. In order to minimize the effect of non-stationarity, we therefore
utilized relative short recordings but with sufficient duration for the systems to be
characterized accurately. Also, the minimal model proposed in this study was designed to
investigate the closed-loop relationships among ABP, ILV and G
PV
. Normally, to
accurately characterize the mechanisms of interest from a system that continuously
operates in a closed-loop manner, one needs to “open the loop” to isolate the component
from the rest of the system. “Opening of the loop” is traditionally achieved by the use of
pharmacological or surgical procedures. Since the measurements made in this study were
obtained when the subjects were in their normal closed-loop physiological condition, the
measured signals contained both the feedback and feedforward dynamics. To isolate each
types of dynamics without having to physically or pharmacologically “open the loop”, we
made use of the fact that there are time delays in each dynamic arm of the model, and
imposed causality constraints on the model structure. By doing so, we essentially
computationally opened the loop (Belozeroff et al. 2002).
We first tried to explain the variability in G
PV
by assuming linear relationships
between MAP and G
PV
through BPC, and ILV and G
PV
through RPC. The dynamics of
the BPC impulse response, h
BPC
, were not entirely consistent with our understanding of
the underlying physiology. In general, one would expect that with an increase in blood
155

pressure, the vascular conductance should also increase (i.e. vasodilation) in order to
dampen the effect of blood pressure. However, the estimated h
BPC
often showed the initial
negative peak but later followed by a positive response that usually occurred at around 5
seconds. One explanation for this finding was that the imposed time delay in the BPC
dynamics was not long enough. As a result, the estimated h
BPC
also captured the strong
feedforward response (vasoconstriction leads to increase in blood pressure) as well as the
feedback response (baroreflex control of peripheral conductance). On the other hand, the
dynamics of the RPC impulse response were as expected – vasoconstriction (decrease in
G
PV
) occurred as the ILV increased. We then proceeded to employing a second-order
nonlinear minimal model to see whether the nonlinear characteristics were involved in
the modulation of G
PV
variability. Despite these possible nonlinear contributions, the
nonlinear model also allowed the characterization of a potential multiplicative interaction
of the two inputs of the system. The nonlinear characteristics were confirmed by
comparing the prediction NMSE of the linear and nonlinear models. The prediction
NMSE significantly reduced when the second-order dynamics were included in the
model. It was observed that the contribution of the second-order kernels were higher for
the RPC self-kernel (h
2RPC
) and the interaction term (h
BPC,RPC
). This indicated that the
nonlinear mechanisms involving ILV interaction had higher contribution to the variability
in G
PV
. The power spectra of the residuals obtained from the nonlinear model showed
that the nonlinear model was able to capture the variability of G
PV
in the low frequency
as well as the breathing frequency range better than the linear model.
156

While there are numerous studies on the baroreflex control of HR dynamics,
fewer have investigated the identification of baroreflex control of TPR (Mukkamala et al.
2003; Aljuri and Cohen 2004; Mukkamala et al. 2006; Chen et al. 2008). This group
proposed a mathematical model of the linear dynamic properties of the arterial and the
cardiopulmonary TPR baroreflexes. Mukkamala et al. (2003) proposed two identification
algorithms: direct and indirect identification. The direct identification estimated arterial
TPR baroreflex by relating fluctuations in ABP to fluctuations in TPR where TPR was
computed from ABP and CO measured by Doppler ultrasound technique described in
Eriksen and Walloe (1990). The dynamics of the cardiopulmonary TPR baroreflex was
estimated from fluctuations in SV, a surrogate measure of right atrial transmural pressure
(RATP) and was computed from the measured CO, and fluctuations in TPR. The indirect
identification, on the other hand, assumed two mechanisms governing the fluctuations in
ABP. The first mechanism related CO to ABP. It assumed that an increase in CO would
cause an increase in ABP through systemic arterial tree, which then activate the arterial
TPR baroreflex as well as the systemic arterial tree arc in order to decrease TPR
(Mukkamala et al. 2003). The second mechanism related SV to ABP. It assumed that SV
fluctuations influenced ABP fluctuations through cardiopulmonary TPR baroreflex as
well as the inverse heart-lung unit and systemic arterial tree (Mukkamala et al. 2003).
Thus, an increase in SV would cause an increase in RATP, which would then activate the
cardiopulmonary TPR baroreflex and systemic arterial tree arc in order to increase TPR
and maintain ABP (Mukkamala et al. 2003). They found that the identification using
indirect identification algorithm was more reliable than the direct identification. Chen et
157

al. (2008) extended the model by Mukkamala et al. (2003; 2006) and they showed that
the impulse response of the arterial baroreflex control of TPR could be estimated reliably.
They showed that the group average estimated arterial TPR baroreflex impulse response
from the control condition exhibited negative feedback dynamics, which was consistent
with the known baroreflex physiology. However, their estimated impulse response was
not associated with a time delay. This finding runs counter to current understanding of
the physiology (Sagawa 1983). Examination of our estimated h
BPC
, demonstrates that the
positive peak, which we speculate to represent the feedback response from the baroreflex,
occurs at approximately 5 seconds. However, most of our estimated h
BPC
contained the
initial negative peak.
In terms of the nonlinear TPR baroreflex behavior, Raymundo et al. (1989)
investigated the arterial and cardiopulmonary baroreflexes in awake dogs and did not find
evidence of nonlinear baroreflex behavior over a wider range than the spontaneous beat-
to-beat fluctuations. However, we found evidence of the nonlinear dynamics of BPC. The
discrepancy in these findings could be due to the species difference. Previous studies
reported the evidence of the second-order nonlinear couplings in the baroreflex control of
HR (Chon et al. 1997; Jo et al. 2007). To our knowledge, we are not aware of other
studies that investigated the identification of the respiratory-peripheral vascular
conductance/resistance coupling. However, other studies have explored the nonlinear
behavior of the coupling between respiration and heart rate variability (Chon et al. 1996;
Jo et al. 2007). Chon et al. (1996) estimated the second-order kernel relating ILV to HR
and found that following the administration of propranolol (sympathetic blockade), the
158

magnitude of the nonlinear kernel became smaller and the reduction in the magnitude
became more prominent with atropine administration (parasympathetic blockade). Jo et
al. (2007) reported the increase in the nonlinear kernel magnitude from wakefulness to
sleep. In our study, we found that the magnitude of the second-order kernel of RPC
decreased from supine to standing. In terms of the cross-kernel, (Chon et al. 1996) found
that the magnitude of the cross-kernel, relating the interaction of ILV and SBP to HR,
decreased when propranolol or atropine was administered. In our case, the magnitude of
the h
BPC,RPC
decreased from supine to standing.
The methods available for noninvasive quantitative assessment of autonomic
function include, but are not limited to, spectral analysis, structured modeling and
minimal modeling. The spectral analysis method offers the ease of computation of the
spectral indices, which undeniably provides important information about the autonomic
functions (Parati et al. 1995). The issue with the spectral analysis is that it only considers
information about the output of the system (Sleight et al. 1995) such as fluctuations in
ABP or in this case would be fluctuations in G
PV
. However, it provides little information
about the underlying mechanisms that contribute to the variability in the output. The
structured models are constructed based on our understanding of the underlying
mechanisms in the system. Therefore, it can offer straightforward physiological
interpretation of the results. However, the danger of this method arises when the assumed
structure is incorrect and the fact that it requires prior knowledge about the system, it is
not readily applicable when we try to identify a system that has never been studied. The
minimal modeling approach allows the characterization of the interactions among
159

biosignals and thus can provide information about the underlying mechanisms.
Furthermore, it also enables us to dissociate the confounding effects of one input from
other inputs that contribute to the variability in the output. In addition, the identification
of the system dynamics is based solely on the input-output data and thus does not require
prior knowledge about the system (Chon et al. 1994). In this study, the estimation of the
kernels representing the dynamics of BPC and RPC was achieved by employing the basis
function expansion technique. Each kernel was constructed by a weighted sum of the
basis functions. The important advantage of this method was that it substantially reduced
the number of unknown parameters to be estimated compared to other methods such as
the autoregressive model. As a result, an increase in estimation accuracy can be achieved
even when applied to relatively short data records with the presence of noise (Marmarelis
1993).
5.4 Kernel Estimation Errors
To validate the accuracy of the kernel estimation, we applied the minimal model
estimation to the simulated “data” produced by Simulation Model B (section 4.5.4). The
estimated impulse responses were then compared with the impulse responses derived
from the associated components of the simulation model. We found that higher system
noise improved the estimation and the estimation was robust to the effect of measurement
noise (assumed to be white Gaussian noise). In addition, the accuracy of the h
BPC

estimation also depended on whether the TPC was taken before or after the addition of
the system noise (Figure 3.14). When the TPC before the addition of the system noise
was employed as the output of the minimal model, the estimation of hBPC was almost
160

identical to the simulated impulse response. On the other hand, when the TPC after the
addition of the system noise was employed as the output, the estimation showed the
initial negative response. However, there was still a positive peak at 5 seconds with the
shape and dynamics that are reasonably comparable to the simulated impulse response.
The negative peak in the estimated h
BPC
could be attributed to the strong feedforward
effect (vasoconstriction leads to an increase in ABP). In the simulation model, the TPR
that went into the Windkessel, which subsequently generated ABP, was a combination
between the calculated TPR and the random system noise. Therefore, this feedforward
effect became amplified when the TPC after the addition of the system was employed as
the output of the minimal model. On a contrary, the TPR baroreflex feedback effect was
amplified when the TPC before the addition of the system was employed as the output of
the minimal model and thus resulted in accurate estimation of h
BPC
. Unfortunately, we do
not have the luxury of separating the feedback and the feedforward effect when taking
measurements from a physiological system without “opening the loop” (i.e. similar to the
case of using TPC after the addition of system noise). With this in mind, this could
explain why our estimation of h
BPC
sometimes showed the initial negative response. Next,
to validate whether the second-order nonlinear minimal model would be able to capture
the nonlinear behavior in the system, the nonlinearity in the TPR baroreflex was
increased by narrowing the linear operating range. We demonstrated that the second-
order nonlinear minimal model was able to capture the increase in nonlinear behavior in
the simulation model as all of the nonlinear kernels showed much more dynamics and
larger amplitudes compared to the kernel estimated from the base simulation model.
161

5.5 Physiological Interpretation of Nonlinear Dynamics
The Volterra-Wiener approach to the system identification offers many practical
advantages. In the linear case, the kernels represent the impulse responses (Marmarelis
and Marmarelis 1978). However, the physiological interpretation of the higher-order
kernels is more complicated. Adopting an approach similar to that introduced by Jo et al.
(2007), we simulated different scenarios using the estimated kernels from different
subject groups to provide more insight into nonlinear behavior represented by the
nonlinear kernels. The simulations of the BPC steady-state response derived from both
the combination of linear and nonlinear kernels demonstrated the saturation effect in the
MAP-TPC relation and the rate of saturation was dependent on the dynamics and the
magnitude of the second-order BPC kernel. The simulations of the RPC response show
that the normalized RPC frequency response was dependent on tidal volume. That is the
larger the magnitude of the second-order RPC kernel, the higher the dependency of the
RPC frequency response upon the tidal volume. Earlier studies showed similar nonlinear
behavior in the RSA (Hirsch and Bishop 1981; Eckberg 1995). Figure 5.2 showed the
tidal volume had significant contribution to RSA response. As the tidal volume increased,
the RSA response also increased.
162

Figure 5.2 Changes in R-R interval at different breathing frequencies and tidal volumes
(Hirsch and Bishop 1981; Eckberg 1995)

The simulations of the interaction between BPC and RPC showed the respiratory
modulation of the baroreflex gain. Both control and MetS + OSAS subjects exhibited
higher BPC peak response during expiration compared to during inspiration. Previous
study by (Eckberg et al. 1985) showed that neck suction applied during expiration caused
a larger increase in the sympathetic activity (Figure 5.3). Another study (Eckberg and
Orshan 1977) investigated the changes in P-P interval in response to as neck suction at
different phase of respiration and found that the response was larger during expiration
compared to during inspiration (Figure 5.4).

163

Figure 5.3 Changes in muscle sympathetic nerve activity (left) after delivering neck
suction at different times in the respiratory cycle (Eckberg 1995).

Figure 5.4 Changes in P-P interval after delivering neck suction at different times in the
respiratory cycle (Eckberg 1995)
164

5.6 System Gains and Autonomic Functions
Many objective tests have been developed for the assessment of autonomic
function. Standard clinical autonomic stress tests include the Valsalva maneuver, heart
rate and blood pressure response due to orthostatic stress, heart rate response to deep
breathing, and blood pressure response to sustained handgrip (Ewing et al. 1985). While
these standard tests are useful, they offer little information about the underlying
mechanisms causing such changes in the test results. Moreover, some of these tests
require patient cooperation, and cannot be applied under certain conditions, such as
during sleep. In the present study, we introduced the BPC and RPC system gains derived
from the minimal model of G
PV
variability as potential autonomic function biomarkers.
We found that in obese pediatric subjects, both BPC and RPC gains of both linear and
second-order nonlinear minimal models decreased in response to orthostatic stress
(Figure 4.38 and Figure 4.39). The reduction in the system gains due to postural change
could be explained by the figure below (Figure 5.5).
In the steady-state during supine posture, the operating point is represented by the
blue dot on the black curve with the slope at that operating point being the gain of the
system in supine posture. Upon standing, blood pressure decreases which subsequently
leads to an increase in sympathetic activity, thus causing the operating point to move up
along the black curve. However, in order to return the blood pressure to the same level as
before, compensatory mechanisms shift the curve to the right (orange curve) and the new
operating point in standing posture would now be at the red dot where the slope or the
gain is lower than it was during supine posture. It could be postulated that the reduction
165

Figure 5.5 Shift in TPR baroreflex during postural change
in the linear gain, reflecting the baroreflex sensitivity, could be “explained” by the shift in
the TPR baroreflex curve. It should be noted that the compensatory shift of the TPR
baroreflex curve is to the right, rather than upwards, since there is likely a ceiling to SNA
activity, especially in subjects that already have high sympathetic tone.
The nonlinear modulation of TPR (and by implication, SNA activity) appears to
follow a trend similar to the linear gains following change in posture from supine to
standing: in standing, the nonlinear gains are diminished. This is reflected by the orange
curve being less skewed (i.e. more linear). The reduction in nonlinear behavior is also
evident from the magnitudes of the time-domain kernels themselves: the kernel
magnitudes in standing posture are much smaller than in supine posture (Figure 4.31 and
Figure 4.32). The interaction, represented by the cross-kernel, also shows less dynamics
in standing compared to supine posture (Figure 4.33). From the obese pediatric subject
Autonomic Input
(ABP or Respiration)
SNA or SVR
Stand: High symp tone
Supine: Low symp tone
Supine
Stand
166

group, two subgroups were selected: the first group had less severe MetS + OSAS,
denoted as the control group; while the second group had more severe MetS + OSAS
indicated by their insulin sensitivity and number of obstructive apnea/hypopnea events
per hour (Figure 4.1). We found that the attenuated nonlinear behavior of BPC was more
pronounced in the MetS + OSAS group than the control group.
In the case of the system gains derived from the sickle cell subjects before and
after blood transfusion. We found that the gains showed the tendency to increase after
blood transfusion. Similar reasoning to the change in posture case could also be applied
here. The pre-transfusion gains would be analogous to the standing posture while the
post-transfusion gains would be analogous to the supine posture (Figure 5.5). Before
receiving blood transfusion treatment, the TPR baroreflex operating point was at the less
steep slope (i.e. lower gain) and exhibited small nonlinear behavior. After the transfusion,
the slope became steeper (i.e. higher gain) and exhibited larger nonlinear behavior (i.e.
more skewed curve). In addtion to the increased nonlinear behavior in BPC after blood
transfusion treatment, the interaction between the blood pressure and respiration also
became more pronouced as shown by the increase in the high-frequency gain derived
from the cross-kernel (Figure 4.46, last row).
To the extent of our knowledge, no other studies have explored the gains of the
TPR baroreflex and the coupling of respiration and peripheral vascular resistance
mechanisms via the minimal modeling approach. However, there have been other studies
on heart rate variability. Belozeroff et al. (2003) investigated the linear gains of the
arterial baroreflex control of heart rate (ABR) and the respiratory-cardiac coupling (RCC)
167

in healthy and untreated obstructive sleep apnea syndrome (OSAS) subjects. Although
they found no significant difffernce in the ABR gains between supine and standing, there
was a significant reduction in the RCC gains in standing compared to supine. They also
found that both ABR and RCC gains were significantly lower in OSAS subjects than in
control subjects. The decrease in the ABR gain in OSAS subjects was attributed to an
elevated sympathetic tone; while the decrease in the RCC gain was attributed to
impairment of the parasympathetic control of heart rate (Belozeroff et al. 2003). Another
similar study, but in male obese pediatric subjects with varying degrees of OSAS
severity, also found similar results that there was a reduction in the RCC gain in standing
versus supine postures, but no statistically signicant difference in the ABR gain between
supine and standing (Oliveira et al. 2009; Oliveira 2011). A study on sickle cell subjects
and the effect of chronic blood transfusion treatment showed that there was no significant
change in the linear ABR or RCC gains after the transfusion treatment (Sangkatumvong
2011).
To validate whether the derived system gains were related to the actual gains in
the system, we estimated the impulse responses from simulated data produced by
Simulation Model B over a wide range of values for the TPR baroreflex gains and RPC
gains. We found that a strong correlation between the derived system gains and the gains
in the simulation model (Figure 4.28). Based on both the experimental evidence as well
as our simulation results, we conclude that the system gains derived from the estimated
kernels of our minimal model of TPR variability are potentially useful as noninvasive
168

biomarkers with sufficient sensitivity to detect the changes in autonomic functions in the
obese pediatric population as well as patients with sickle cell disease.
5.7 Limitations of the Study
5.7.1 Simulation Model
In terms of the simulation model, although the simulation model developed in this
study was able to simulate the responses similar to what would be expected in a
physiological system under normal conditions, it still lacked a number of major
regulatory mechanisms that also contribute to the blood pressure control. Other
regulatory systems involved in short-term regulation of blood pressure include, but not
limited to, hormonal regulation, renal regulation, microvascular circulation, and, in
particular, the cardiopulmonary reflexes as well as the chemoreflexes. Furthermore, all
elements in this simulation model were either static linear or static nonlinear. However,
physiological systems are well-known for being non-stationary. An example of a non-
stationary mechanism that was not included in this model would be arterial baroreflex
resetting. In addition, to ensure that the simulation model could also operate under
extreme conditions such as hemorrhage or exercise, threshold and saturation level would
be required to ensure that the computations would not run into numerical errors.
5.7.2 Minimal Models
An important limitation in our minimal modeling approach was the assumption of
stationarity, while physiological systems are inherently non-stationary in behavior. The
next point that should be noted is pertinent to the model inputs. In principle, to achieve
169

good estimation of the kernels, the model inputs should be independent of each other
(Chon et al. 1996). However, ABP, respiration, TPR and other cardiovascular variables
such as HR each plays a part of a larger and complex closed-loop control system and so
they are in some way dependent on each other. In addition, an ideal input to the system
should contain many frequency components to be able to broadly stimulate the system
and thus allow more accurate estimation of the system. Although the inputs employed in
this study (MAP and ILV) contain certain degrees of variability but these are the
spontaneously occurring variability in a normal physiological system. Thus, the
frequency content in these inputs may not be broad enough to stimulate the system
sufficiently, and there is clearly some correlation between the inputs. To overcome both
issues, we could introduce a random breathing frequency protocol to broaden the spectral
component of the cardiovascular variability (Berger et al. 1989). Regarding the blood
pressure signal to be employed in the estimation model, DBP could be one of the
promising factors as reported by Eckberg (1995) that muscle sympathetic traffic became
higher as DBP decreased (Figure 5.6). However, the blood pressure measurements in this
study were derived from peripheral sites. The peripheral waveform is known to be
distorted by the pressure wave transmission and reflection in the arterial tree. Therefore,
MAP was used instead as the integration of the blood pressure contour over each beat
would result in some degree of cancellation of the waveform distortions (Elstad et al.
2001).
Lastly, we found that the residuals of both the linear and the second-order
nonlinear model were usually not white. This suggests that there are still other
170

mechanisms or higher-order of nonlinearity that contribute to the variability in R
PV
but
were not included in the minimal models proposed in this study. Other studies have
shown that the cardiopulmonary baroreflex also plays a great role in regulating R
PV

(Mukkamala et al. 2003; Hughson et al. 2004).

Figure 5.6 Muscle sympathetic traffic measured at different diastolic pressure levels
during one breath (Eckberg 1995)
5.7.3 Surrogate Measures of TPR
To assess the sympathetic activation in response to baroreflex stimulation or
respiratory modulation, it would be ideal if we could obtain the blood pressure or
respiration (as an input) and the SNA measurement (as an output) directly, continuously,
171

as well as noninvasively. To date, the method considered to be a gold standard of
measuring sympathetic nervous activity has been the radiotracer technology which
measures the norepinephrine spillover (Lambert et al. 1998; Malpas 2010). However, this
technique only provides single and regional measurement and thus would not be
applicable for the modeling purpose which requires continuous measurement. Another
widely used technique for measuring the sympathetic nerve activity is the direct
measurement of the MSNA from the peroneal nerve. This method offers continuous
recording of SNA and was found to have correlation with cardiac output, and thus also
TPR, in males (Charkoudian et al. 2005; Malpas 2010). In addition, it had also been
reported to have good correlation with the norepinephrine spillover at baseline (Wallin et
al. 1992). However, the correlation between MSNA and TPR was not found in females
and the relationship with the norepinephrine spillover does not hold for all conditions
(Malpas 2010). Other disadvantages of this technique include its invasiveness and the
requirement of technical expertise to perform the MSNA recording, making it not readily
applicable for clinical usage. For these reasons, we proposed that instead of measuring
the SNA itself, we could employ the vasoconstriction response measurement since the
vasoconstriction response reflects changes in TPR and changes in TPR are known to be
largely sympathetically mediated (Malpas 2002).
One drawback of using TPR as an index of vascular resistance is that TPR is a
highly global measure – sympathetic stimulation can lead to vasoconstriction or
vasodilation in various arterial/capillary beds, depending on the types of receptors (alpha
versus beta) present in those regional beds. In this study, we utilized PAT and laser
172

Doppler flowmetry as the surrogate measures of peripheral vascular resistance, R
PV.
By
the very nature of these measurements, they provide information about the vascular
resistance in the associated local regions. Previous studies have demonstrated that
peripheral vasoconstriction response could be measured noninvasively by PAT (Schnall
et al. 1999; O'Donnell et al. 2002; Kuvin et al. 2003; Rubinshtein et al. 2010) and laser
Doppler flowmetry (Inwald et al. 1996; Aso et al. 1997; Allen et al. 2002; Mayrovitz and
Groseclose 2002). PAT, however, does not directly measure TPR but rather the changes
in volume at the fingertip (Kuvin et al. 2003). The measurement of vasoconstriction
response using PAT is considered to be a regional response. Furthermore, whether the
relationship between R
PV
and amplitude of PAT is linear or not is yet to be determined.
That is to say the change in the PAT amplitude due to an increase in R
PV
by one unit may
not be the same as the change in the amplitude of PAT due to a decrease in R
PV
of one
unit. In addition, PAT measurements are given in arbitrary units, making it inconvenient
to compare across subjects. PAT signals collected from different postures were
normalized to each subject’s own baseline in supine. Therefore, the information that can
be derived from PAT would be the measure of reactivity, such as how its gain changes
with posture, rather than absolute gain, such as baroreflex sensitivity in each posture.
While PAT is considered to be a regional measurement, the laser Doppler flowmetry
measures the microvascular perfusion on the skin (Oberg et al. 1984), and thus focuses on
an even more local response compared to PAT. Similar to PAT, microvascular perfusion
measured by the laser Doppler flowmetry is also given in arbitrary units. Lastly, we
should also keep in mind that while both PAT and laser Doppler flowmetry provide
173

surrogate measures of R
PV
but not SNA itself, since vasoconstriction is an effector
response to sympathetic outflow. Consequently, it is possible for reductions in PAT
amplitude or microvascular perfusion to not always reflect low sympathetic gain, but
rather, these changes could also be due to diminished effector gain.
Nexfin (BMEYE, Amsterdam, The Netherlands) offers a measurement of TPR
(or systemic vascular resistance, SVR, as denoted by Nexfin) as one of its output. Nexfin
is a noninvasive continuous blood pressure monitoring device. It estimates cardiac output
from the finger arterial pressure using pulse contour analysis (Wesseling et al. 1993;
Truijen et al. 2012) and, therefore, TPR is estimated from cardiac output and blood
pressure. Although this Nexfin parameter offers an alternative measure of TPR, it is not a
direct measure of TPR but rather, an estimate based on the decay phase of the arterial
pressure pulse measured at the finger.

174

Chapter 6. Future Work and Conclusions
6.1 Future Work
6.1.1 Closed-loop Model of Cardiovascular Variability
The minimal model of R
PV
variability proposed in this study could be extended to
incorporate the other parts of the cardiovascular regulation such as the heart rate
baroreflex control, the respiratory-cardiac coupling, the mechanical effect of respiration
on blood pressure (Belozeroff et al. 2002; Khoo 2008; Chaicharn et al. 2009) into the
larger closed-loop model of cardiovascular variability as presented in Figure 1.2. This
closed-loop model would give a more comprehensive picture of the cardiovascular
control system and subsequently let us gain more insight into the dynamics of the
regulation of this system in individual subjects.
6.1.2 Linear Time-varying Model
In addition to the fact that most physiological systems are generally non-
stationary in nature, a time-varying model would also offer the ability to track responses
induced by the experimental interventions such as cold face stimulation. Using similar
methodology employed by Chaicharn et al. (2009), the linear minimal model of G
PV

variability could be extended to incorporate time-varying BPC and RPC dynamics as
follows
175

( ) ( ) ( )
( ) ( ) ( ) t T i t ILV i t h
T i t MAP i t h t G
PV
M
i
RPC RPC
M
i
BPC BPC PV
σ
ε + − − ∆ ⋅ +
− − ∆ ⋅ = ∆
∑
∑
−
=
−
=
1
0
1
0
,
,

Equation 6.1
where h
BPC
(t,i) and h
RPC
(t,i) are the time-varying impulse responses. Similar to the time-
invariant model, the kernel expansion technique (section 3.3.3) can be applied such that
Equation 6.1 can be rewritten as follows

( ) ( ) ( ) ( ) ( ) ( ) t t v t c t u t c t G
PV
RPC BPC
j
q
j
RPC
j j
q
j
BPC
j PV σ
ε + + = ∆
∑ ∑
= = 1 1

Equation 6.2
where ) (t c
BPC
j
and ) (t c
RPC
j
are the time-varying expansion coefficients to be estimated.
These coefficients can be estimated using an autoregressive model with an adaptive
parameter estimation algorithm, such as recursive least-squares (Ljung 1999).
The forgetting factor approach would be employed for the recursive estimation.
The forgetting factor (λ) reflects the memory of the adaptive filter and can vary from 0 to
1. When λ = 1 all data before the present time are used in the estimation of the current
parameter. Small λ implies that the most recent data points are weighted more heavily.
With this technique, we can adjust λ to allow the model prediction to be less sensitive to
the remote past input. λ is chosen such that the cost function J
w
is minimized (Chaicharn
et al. 2009).
) (
2
1
0
i t e J
t
i
i
w
− =
∑
−
=
λ
Equation 6.3
where e is the error between the observed data and the prediction.
176

The linear time-varying model would be useful for tracking responses that are
time-varying in nature such as response to cold face stimulation, or the transition between
the changes in posture from supine to standing using tilt table; or tracking responses of
the systems that are known to be changing with time even without any experimental
interventions.
6.1.3 Principal Dynamic Modes
The Volterra-Wiener approach has shown to be useful in modeling the
physiological systems. However, as the model order becomes higher, the estimated
kernels also increase in dimensionality, making the estimation as well as the
physiological interpretation more difficult (Marmarelis et al. 2012). To overcome this
limitation, the concept of principal dynamic modes (PDM) – that the dynamics of a
system can be represented by a minimum set of basis functions that are unique to that
system (Marmarelis 1997) – was introduced. PDM approach has been applied in various
applications, for example, analysis of renal autoregulation in rats (Marmarelis et al.
1999), contributions of autonomic nervous system to heart rate variability (Zhong et al.
2004), and cerebral flow autoregulation (Marmarelis et al. 2012). Using the PDM
technique, the Volterra-Wiener model can be represented in a more compact form. At the
same time, PDM analysis retains the essential dynamic characteristics of the estimated
system while rejecting kernel contributions that could be due to the corrupting effects of
noise. Furthermore, PDM analysis, in principle, also allows the incorporation of high-
order (>2) without making the system in question non-identifiable due to over-
parameterization.
177

6.1.4 Multi-input Minimal Model of R
PV

Peripheral blood flow, a surrogate measure of R
PV
, is known to be controlled by
both extrinsic and intrinsic pathways (Levy et al. 2007). The extrinsic pathway is mainly
mediated by sympathetic nerve signaling from the central nervous system and influences
peripheral blood flow by regulating vascular resistance through arterial baroreflex. This
part of the mechanism is already included in the current minimal model. The intrinsic
pathway, however, is controlled by the local effect of nitric oxide, which acts to dilate the
vessels. By incorporating this intrinsic control of the peripheral blood flow, we may be
able to dissociate the effect of the sympathetic control from the effect of the endothelial
local control of the R
PV
. The model structure was proposed by Sangkatumvong (2011).
Another mechanism that also deserves further investigation is the cardiopulmonary
baroreflex. Other studies reported that it also plays a significant role in the regulation of
R
PV
(Mukkamala et al. 2003; Hughson et al. 2004). Besides the mechanisms involved in
the regulation of R
PV
, we can also introduce an external input such as pain or other
stressors, which are used as the stimuli in the experiments into the model.
6.1.5 Reduction in Parameter Estimate Variability
In the minimal modeling applied in this work, one set of estimated kernel
coefficients was determined from each dataset. The estimated parameters did not contain
any error bounds or confidence limits. We believe it would be useful to introduce a
bootstrapping method that can provide a rough indication of the error associated with
each estimated model parameter. We have performed preliminary trials with an iterative
method, as detailed in this section.
178

Figure 6.1 Procedure to generate surrogate data using AAFT method
The Amplitude-adjusted Fourier transform (AAFT) is a method for generating
surrogate data (Theiler et al. 1992), which still retains the distribution as well as
approximately the power spectrum of the original signal (Dolan and Spano 2001). The
Original signal, x(t)
Gaussian random signal, y(t),
with the same variance as x(t)
Amplitude-adjusted y(t), y
AA
(t)
Surrogate y(t), y
surr
(t)
Magnitude of y
AA
(t), |Y
AA
(f)|
Random Phase
Surrogate x(t), x
surr
(t)
|X(f)|
|X
surr
(f)|
Magnitude of x(t) is preserved.
F
F
-1
F
179

steps taken to generate a set of surrogate data using the AAFT method are illustrated in
Figure 6.1. In brief, the Fourier transform is first applied to the original data. To preserve
its power, only the phase component of the original data is randomized then combined
with the original magnitude component. Applying inverse Fourier transform, the
surrogate data is obtained. By applying the AAFT method to the model residuals, the
surrogate residuals can then be added to the model prediction to generate a “new” output
time-series and the new estimation begins again. Thus, this technique allows us to
artificially “repeat” the experiments. Using the estimation from all iterations, the
confidence interval of the estimation can be constructed, give a sense of how reliable the
estimation is. Figure 6.2 shows the estimated BPC and RPC impulse responses from a
sickle cell disease subject. Each grey lines show the estimation with an application of
AAFT method (25 iterations).

Figure 6.2 Estimated h
BPC
and h
RPC
of a sickle cell disease subject. Thick red and blue
tracing are the estimated impulse responses using the actual data. Thin grey lines are the
estimations using the “new” output data, which is generated using the AAFT method.
0 5 10 15 20 25
-2
-1
0
1
2
3
4
x 10
-4
h
BPC
Time (sec)
0 5 10 15 20 25
-15
-10
-5
0
x 10
-4
h
RPC
Time (sec)
180

6.2 Conclusions
In the present study, our goal is to develop a computational model of R
PV

regulation, and to use this model to provide information about sympathetic modulation of
the peripheral vasculature. SNA plays a major role in autonomic regulation as such
assessment of SNA has been the research focus since the knowledge of sympathetic tone
provides information regarding the underlying autonomic physiology as well as the
clinical state of the subject being tested. To date, many techniques have been developed
to quantify the SNA but they are oftentimes invasive, require technical expertise as well
as can be costly, making them not readily applicable for clinical usage. The regulation of
TPR is known to be largely sympathetically mediated and the changes in TPR are
reflected in the vasoconstriction response. For these reasons, another goal of this study is
to employ a noninvasive measurement of TPR in the computational model, in which
sympathetic modulation can be inferred.
The variability in TPR is generally attributed to baroreflex control, although there
is abundant evidence of respiratory modulation of sympathetic outflow. Additionally,
many studies have reported that deep breaths, or sighs, can lead to peripheral
vasoconstriction and yet the respiratory modulation of TPR has been little studied. Thus
we would also seek to identify how respiration plays a role in the modulation of TPR –
whether it comes from the effect of respiratory modulation on ABP or from the direct
modulation of respiration. To investigate the two main mechanisms responsible for the
variability in TPR: baroreflex control and the respiratory modulation of TPR, we adopted
PAT and laser Doppler flowmetry as the surrogate measures of R
PV
. Both PAT and laser
181

Doppler flowmetry have been reported to be able to detect peripheral vasoconstriction
response in various applications. PAT measures the volume changes due to blood
pressure changes at the finger tip. Vasoconstriction would lead to attenuation in the pulse
amplitude in the PAT signal and thus the amplitude of the PAT signal can be employed
as an indicator of vasoconstriction response. Laser Doppler flowmetry measures the skin
microvascular perfusion. Vasoconstriction would be reflected as a reduction in the blood
flow and perfusion.
By means of structured and minimal modeling approaches, we found that
respiration likely affects R
PV
through a direct modulation and this direct respiratory
modulation likely occurs spontaneously with or without the presence of a sigh, in spite of
the fact that its effect only becomes visible when there is a sigh. We demonstrated using a
simulation model, which consists of main components in the cardiovascular system
including the modulation of respiration on blood pressure, that the sigh-vasoconstriction
could not be reproduced unless an additional mechanism of direct respiratory modulation
on TPR was incorporated into the simulation model. Applying these models to the
physiological measurements obtained from obese pediatric subjects with varying degrees
of severity in metabolic syndrome and OSAS during orthostatic stress revealed that the
linear gains of both the baroreflex control of peripheral vascular conductance (BPC) and
the respiratory-peripheral vascular conductance coupling (RPC) significantly decreased
from supine to standing. These results suggest that their sympathetic modulation of R
PV

diminished during standing despite the increase in sympathetic tone.
182

We further investigated the nonlinear dynamics and the interaction effect involved
in the modulation of R
PV
through BPC and RPC mechanisms in both obese pediatric
subjects exposed to orthostatic stress as well as sickle cell disease subjects before and
after blood transfusion treatment. In both subject groups, the results reveal that both BPC
and RPC exhibited nonlinear behavior. There also existed the interaction between the
dynamics of BPC and RPC. By including the nonlinear dynamics and the interaction
effect, in addition to the linear dynamics, the model prediction error was significantly
decreased. Mechanisms with large contribution to R
PV
variability were the second-order
nonlinear RPC dynamics and the interaction effect between ABP and respiration terms.
In the obese pediatric group, the nonlinear behavior became smaller when exposed to
orthostatic stress. Subjects with more severe degrees of metabolic syndrome and OSAS
showed larger reduction in nonlinear behavior of the baroreflex control of TPR. In the
sickle cell disease group, blood transfusion treatment led to an increase in nonlinear TPR
baroreflex gain as well as the interaction between BPC and RPC dynamics.
In conclusion, in the present study, we have developed a minimal model of R
PV

regulation using the noninvasive PAT and laser Doppler flowmetry measurements.
Through a combination of the structured and minimal modeling, we have demonstrated
the direct influence of respiration on R
PV
but the mechanism behind this modulation
remains to be investigated. Taking this respiratory modulation effect into account is
essential for achieving accurate TPR baroreflex estimation. Finally, the system gains
derived from the estimated kernels may constitute potentially useful biomarkers of
sympathetic nervous system function.
183

References
Akselrod S., Gordon D., Madwed J. B., Snidman N. C., Shannon D. C. and Cohen R. J.
(1985). Hemodynamic regulation: investigation by spectral analysis. The
American journal of physiology 249, H867-875.
Aljuri N. and Cohen R. J. (2004). Theoretical considerations in the dynamic closed-loop
baroreflex and autoregulatory control of total peripheral resistance. American
journal of physiology. Heart and circulatory physiology 287, H2252-2273.
Allen J., Frame J. R. and Murray A. (2002). Microvascular blood flow and skin
temperature changes in the fingers following a deep nspiratory gasp.
Physiological measurement 23, 365-373.
Almond N. E. and Wheatley A. M. (1992). Measurement of hepatic perfusion in rats by
laser Doppler flowmetry. The American journal of physiology 262, G203-209.
Arvidsson D., Svensson H. and Haglund U. (1988). Laser-Doppler flowmetry for
estimating liver blood flow. The American journal of physiology 254, G471-476.
Aso Y., Inukai T. and Takemura Y. (1997). Evaluation of skin vasomotor reflexes in
response to deep inspiration in diabetic patients by laser Doppler flowmetry. A
new approach to the diagnosis of diabetic peripheral autonomic neuropathy.
Diabetes care 20, 1324-1328.
Asyali M. H. and Juusola M. (2005). Use of meixner functions in estimation of Volterra
kernels of nonlinear systems with delay. IEEE transactions on bio-medical
engineering 52, 229-237.
Baron R., Habler H. J., Heckmann K. and Porschke H. (1996). Respiratory modulation of
blood flow in normal and sympathectomized skin in humans. J Auton Nerv Syst
60, 147-153.
Baselli G., Cerutti S., Civardi S., Malliani A. and Pagani M. (1988). Cardiovascular
variability signals: towards the identification of a closed-loop model of the neural
control mechanisms. IEEE transactions on bio-medical engineering 35, 1033-
1046.
Belozeroff V., Berry R. B. and Khoo M. C. (2003). Model-based assessment of
autonomic control in obstructive sleep apnea syndrome. Sleep 26, 65-73.

184

Belozeroff V., Berry R. B., Sassoon C. S. and Khoo M. C. (2002). Effects of CPAP
therapy on cardiovascular variability in obstructive sleep apnea: a closed-loop
analysis. American journal of physiology. Heart and circulatory physiology 282,
H110-121.
Berger R. D., Saul J. P. and Cohen R. J. (1989). Assessment of autonomic response by
broad-band respiration. IEEE transactions on bio-medical engineering 36, 1061-
1065.
Berntson G. G., Bigger J. T., Jr., Eckberg D. L., Grossman P., Kaufmann P. G., Malik
M., Nagaraja H. N., Porges S. W., Saul J. P., Stone P. H. and van der Molen M.
W. (1997). Heart rate variability: origins, methods, and interpretive caveats.
Psychophysiology 34, 623-648.
Bogert L. W., Harms M. P., Pott F., Secher N. H., Wesseling K. H. and van Lieshout J. J.
(2004). Reconstruction of brachial pressure from finger arterial pressure during
orthostasis. J Hypertens 22, 1873-1880.
Bolton B., Carmichael E. A. and Sturup G. (1936). Vaso-constriction following deep
inspiration. J Physiol 86, 83-94.
Broch O., Renner J., Gruenewald M., Meybohm P., Schottler J., Caliebe A., Steinfath M.,
Malbrain M. and Bein B. (2012). A comparison of the Nexfin(R) and
transcardiopulmonary thermodilution to estimate cardiac output during coronary
artery surgery. Anaesthesia 67, 377-383.
Browse N. L. and Hardwick P. J. (1969). The deep-breath-venoconstriction reflex. Clin
Sci 37, 125-135.
Chaicharn J., Lin Z., Chen M. L., Ward S. L., Keens T. and Khoo M. C. (2009). Model-
based assessment of cardiovascular autonomic control in children with obstructive
sleep apnea. Sleep 32, 927-938.
Chalacheva P. and Khoo M. C. (2013). An extended model of blood pressure variability:
Incorporating the respiratory modulation of vascular resistance. Conf Proc IEEE
Eng Med Biol Soc 2013, 3825-3828.
Charkoudian N., Joyner M. J., Johnson C. P., Eisenach J. H., Dietz N. M. and Wallin B.
G. (2005). Balance between cardiac output and sympathetic nerve activity in
resting humans: role in arterial pressure regulation. J Physiol 568, 315-321.

185

Chen X., Kim J. K., Sala-Mercado J. A., Hammond R. L., Elahi R. I., Scislo T. J., Swamy
G., O'Leary D. S. and Mukkamala R. (2008). Estimation of the total peripheral
resistance baroreflex impulse response from spontaneous hemodynamic
variability. American journal of physiology. Heart and circulatory physiology
294, H293-301.
Cheng L., Ivanova O., Fan H. H. and Khoo M. C. (2010). An integrative model of
respiratory and cardiovascular control in sleep-disordered breathing. Respiratory
physiology & neurobiology 174, 4-28.
Chon K. H., Chen Y. M., Marmarelis V. Z., Marsh D. J. and Holstein-Rathlou N. H.
(1994). Detection of interactions between myogenic and TGF mechanisms using
nonlinear analysis. The American journal of physiology 267, F160-173.
Chon K. H., Mukkamala R., Toska K., Mullen T. J., Armoundas A. A. and Cohen R. J.
(1997). Linear and nonlinear system identification of autonomic heart-rate
modulation. IEEE Eng Med Biol Mag 16, 96-105.
Chon K. H., Mullen T. J. and Cohen R. J. (1996). A dual-input nonlinear system analysis
of autonomic modulation of heart rate. IEEE transactions on bio-medical
engineering 43, 530-544.
deBoer R. W., Karemaker J. M. and Strackee J. (1987). Hemodynamic fluctuations and
baroreflex sensitivity in humans: a beat-to-beat model. The American journal of
physiology 253, H680-689.
Dolan K. T. and Spano M. L. (2001). Surrogate for nonlinear time series analysis. Phys
Rev E Stat Nonlin Soft Matter Phys 64, 046128.
Eckberg D. L. (1983). Human sinus arrhythmia as an index of vagal cardiac outflow. J
Appl Physiol 54, 961-966.
Eckberg D. L. (1995). Respiratory Sinus Arrhythmia and Other Human Cardiovascular
Neural Periodicities. In Regulation of Breathing, 2nd edn, ed. J. A. Dempsey and
A. I. Pack, pp. 669-740. Marcel Dekker, New York.
Eckberg D. L., Nerhed C. and Wallin B. G. (1985). Respiratory modulation of muscle
sympathetic and vagal cardiac outflow in man. J Physiol 365, 181-196.
Eckberg D. L. and Orshan C. R. (1977). Respiratory and baroreceptor reflex interactions
in man. The Journal of clinical investigation 59, 780-785.

186

Eeftinck Schattenkerk D. W., van Lieshout J. J., van den Meiracker A. H., Wesseling K.
R., Blanc S., Wieling W., van Montfrans G. A., Settels J. J., Wesseling K. H. and
Westerhof B. E. (2009). Nexfin noninvasive continuous blood pressure validated
against Riva-Rocci/Korotkoff. Am J Hypertens 22, 378-383.
Elstad M., Toska K., Chon K. H., Raeder E. A. and Cohen R. J. (2001). Respiratory sinus
arrhythmia: opposite effects on systolic and mean arterial pressure in supine
humans. J Physiol 536, 251-259.
Eriksen M. and Walloe L. (1990). Improved method for cardiac output determination in
man using ultrasound Doppler technique. Med Biol Eng Comput 28, 555-560.
Esler M., Jennings G., Korner P., Willett I., Dudley F., Hasking G., Anderson W. and
Lambert G. (1988). Assessment of human sympathetic nervous system activity
from measurements of norepinephrine turnover. Hypertension 11, 3-20.
Eun H. C. (1995). Evaluation of skin blood flow by laser Doppler flowmetry. Clin
Dermatol 13, 337-347.
Ewing D. J., Martyn C. N., Young R. J. and Clarke B. F. (1985). The value of
cardiovascular autonomic function tests: 10 years experience in diabetes. Diabetes
care 8, 491-498.
Fukushima H., Kelstrup J., Fukushima S., Umemoto T. and Sagawa H. (1983). Isolation,
partial purification and preliminary characterization of a bacteriocin from
Streptococcus mutans Rm-10. Antonie Van Leeuwenhoek 49, 41-50.
Galland B. C., Taylor B. J., Bolton D. P. and Sayers R. M. (2000). Vasoconstriction
following spontaneous sighs and head-up tilts in infants sleeping prone and
supine. Early Hum Dev 58, 119-132.
Goldstein D. S., Chang P. C., Eisenhofer G., Miletich R., Finn R., Bacher J., Kirk K. L.,
Bacharach S. and Kopin I. J. (1990). Positron emission tomographic imaging of
cardiac sympathetic innervation and function. Circulation 81, 1606-1621.
Goor D. A., Sheffy J., Schnall R. P., Arditti A., Caspi A., Bragdon E. E. and Sheps D. S.
(2004). Peripheral arterial tonometry: a diagnostic method for detection of
myocardial ischemia induced during mental stress tests: a pilot study. Clinical
cardiology 27, 137-141.
Grassi G. and Esler M. (1999). How to assess sympathetic activity in humans. J
Hypertens 17, 719-734.
187

Grote L., Zou D., Kraiczi H. and Hedner J. (2003). Finger plethysmography--a method
for monitoring finger blood flow during sleep disordered breathing. Respiratory
physiology & neurobiology 136, 141-152.
Hannemann M. M., Liddell W. G., Shore A. C., Clark P. M. and Tooke J. E. (2002).
Vascular function in women with previous gestational diabetes mellitus. J Vasc
Res 39, 311-319.
Henriksen L. and Sejrsen P. (1976). Local reflex in microcirculation in human cutaneous
tissue. Acta physiologica Scandinavica 98, 227-231.
Hirsch J. A. and Bishop B. (1981). Respiratory sinus arrhythmia in humans: how
breathing pattern modulates heart rate. The American journal of physiology 241,
H620-629.
Hughson R., O'Leary D., Shoemaker J. K., Lin D., Topor Z., Edwards M. and Tulppo M.
(2004). Searching for the Vascular Component of the Arterial Baroreflex.
Cardiovascular Engineering: An International Journal 4, 155-162.
Imholz B. P., Wieling W., van Montfrans G. A. and Wesseling K. H. (1998). Fifteen
years experience with finger arterial pressure monitoring: assessment of the
technology. Cardiovascular research 38, 605-616.
Inwald D., Hathorn M. K. and Costeloe K. (1996). The deep breath vasoconstriction
reflex--a new tool for autonomic assessment in infancy? Early Hum Dev 45, 55-
61.
Jepsen H. and Gaehtgens P. (1993). Postural vascular response in human skin: passive
and active reactions to alteration of transmural pressure. The American journal of
physiology 265, H949-958.
Jo J. A., Blasi A., Valladares E. M., Juarez R., Baydur A. and Khoo M. C. (2007). A
nonlinear model of cardiac autonomic control in obstructive sleep apnea
syndrome. Ann Biomed Eng 35, 1425-1443.
Julien C. (2006). The enigma of Mayer waves: Facts and models. Cardiovascular
research 70, 12-21.
Kandel E. R., Schwartz J. H. and Jessell T. M. (2000). Principles of neural science. New
York, McGraw-Hill, Health Professions Division.
Kanters J. K., Hojgaard M. V., Agner E. and Holstein-Rathlou N. H. (1996). Short- and
long-term variations in non-linear dynamics of heart rate variability.
Cardiovascular research 31, 400-409.
188

Katona P. G. and Jih F. (1975). Respiratory sinus arrhythmia: noninvasive measure of
parasympathetic cardiac control. J Appl Physiol 39, 801-805.
Khoo M. C. (2008). Modeling of autonomic control in sleep-disordered breathing.
Cardiovasc Eng 8, 30-41.
Kuvin J. T., Patel A. R., Sliney K. A., Pandian N. G., Sheffy J., Schnall R. P., Karas R.
H. and Udelson J. E. (2003). Assessment of peripheral vascular endothelial
function with finger arterial pulse wave amplitude. Am Heart J 146, 168-174.
La Rovere M. T., Mortara A., D. P. G. and Bernardi L. (1995). Baroreflex sensitivity and
heart rate variability in the assessment of the autonomic status. In Heart Rate
Variability, ed. M. Malik and A. J. Camm, pp. 189-205. NY Futura Publishing
Company, Inc., New York.
Lambert G. W., Kaye D. M., Thompson J. M., Turner A. G., Cox H. S., Vaz M., Jennings
G. L., Wallin B. G. and Esler M. D. (1998). Internal jugular venous spillover of
noradrenaline and metabolites and their association with sympathetic nervous
activity. Acta physiologica Scandinavica 163, 155-163.
Lanfranchi P. A. and Somers V. K. (2002). Arterial baroreflex function and
cardiovascular variability: interactions and implications. Am J Physiol Regul
Integr Comp Physiol 283, R815-826.
Levy M. N., Pappano A. J. and Berne R. M. (2007). Cardiovascular physiology.
Philadelphia, PA, Mosby Elsevier.
Lindsberg P. J., O'Neill J. T., Paakkari I. A., Hallenbeck J. M. and Feuerstein G. (1989).
Validation of laser-Doppler flowmetry in measurement of spinal cord blood flow.
The American journal of physiology 257, H674-680.
Ljung L. (1999). System identification : theory for the user. Upper Saddle River, NJ,
Prentice Hall PTR.
Madwed J. B., Albrecht P., Mark R. G. and Cohen R. J. (1989). Low-frequency
oscillations in arterial pressure and heart rate: a simple computer model. The
American journal of physiology 256, H1573-1579.
Malpas S. C. (2002). Neural influences on cardiovascular variability: possibilities and
pitfalls. American journal of physiology. Heart and circulatory physiology 282,
H6-20.
Malpas S. C. (2010). Sympathetic nervous system overactivity and its role in the
development of cardiovascular disease. Physiol Rev 90, 513-557.
189

Malpas S. C., Leonard B. L., Guild S. J., Ringwood J. V., Navakatikyan M., Austin P. C.,
Head G. A. and Burgess D. E. (2001). The sympathetic nervous system's role in
regulating blood pressure variability. IEEE Eng Med Biol Mag 20, 17-24.
Marmarelis P. Z. and Marmarelis V. Z. (1978). The White-Noise Method in System
Identification. In Analysis of Physiological Systems: The White-Noise Approach,
pp. 131-180. Plenum Press.
Marmarelis V., Shin D. and Zhang R. (2012). Linear and nonlinear modeling of cerebral
flow autoregulation using principal dynamic modes. Open Biomed Eng J 6, 42-55.
Marmarelis V. Z. (1993). Identification of nonlinear biological systems using Laguerre
expansions of kernels. Ann Biomed Eng 21, 573-589.
Marmarelis V. Z. (1997). Modeling methodology for nonlinear physiological systems.
Ann Biomed Eng 25, 239-251.
Marmarelis V. Z. (2004). Nonlinear Dynamic Modeling of Physiological Systems.
Hoboken, NJ, Wiley-Interscience.
Marmarelis V. Z., Chon K. H., Chen Y. M., Marsh D. J. and Holstein-Rathlou N. H.
(1993). Nonlinear analysis of renal autoregulation under broadband forcing
conditions. Ann Biomed Eng 21, 591-603.
Marmarelis V. Z., Chon K. H., Holstein-Rathlou N. H. and Marsh D. J. (1999). Nonlinear
analysis of renal autoregulation in rats using principal dynamic modes. Ann
Biomed Eng 27, 23-31.
Mayrovitz H. N. and Groseclose E. E. (2002). Inspiration-induced vascular responses in
finger dorsum skin. Microvasc Res 63, 227-232.
Mukkamala R., Kim J. K., Li Y., Sala-Mercado J., Hammond R. L., Scislo T. J. and
O'Leary D. S. (2006). Estimation of arterial and cardiopulmonary total peripheral
resistance baroreflex gain values: validation by chronic arterial baroreceptor
denervation. American journal of physiology. Heart and circulatory physiology
290, H1830-1836.
Mukkamala R., Mathias J. M., Mullen T. J., Cohen R. J. and Freeman R. (1999). System
identification of closed-loop cardiovascular control mechanisms: diabetic
autonomic neuropathy. The American journal of physiology 276, R905-912.
Mukkamala R., Toska K. and Cohen R. J. (2003). Noninvasive identification of the total
peripheral resistance baroreflex. American journal of physiology. Heart and
circulatory physiology 284, H947-959.
190

Mullen T. J., Appel M. L., Mukkamala R., Mathias J. M. and Cohen R. J. (1997). System
identification of closed-loop cardiovascular control: effects of posture and
autonomic blockade. The American journal of physiology 272, H448-461.
Neff R. A., Wang J., Baxi S., Evans C. and Mendelowitz D. (2003). Respiratory sinus
arrhythmia: endogenous activation of nicotinic receptors mediates respiratory
modulation of brainstem cardioinhibitory parasympathetic neurons. Circulation
research 93, 565-572.
Nilsson G. E. (1984). Signal processor for laser Doppler tissue flowmeters. Med Biol Eng
Comput 22, 343-348.
O'Donnell C. P., Allan L., Atkinson P. and Schwartz A. R. (2002). The effect of upper
airway obstruction and arousal on peripheral arterial tonometry in obstructive
sleep apnea. Am J Respir Crit Care Med 166, 965-971.
Oberg P. A., Tenland T. and Nilsson G. E. (1984). Laser-Doppler flowmetry--a non-
invasive and continuous method for blood flow evaluation in microvascular
studies. Acta Med Scand Suppl 687, 17-24.
Oliveira F. M. (2011). Autonomic and Metabolic Effects of Obstructive Sleep Apnea in
Childhood Obesity. Doctor of Philosophy Dissertation, University of Southern
California.
Oliveira M. M., da Silva N., Timoteo A. T., Feliciano J., Silva S., Xavier R., Rocha I.,
Silva-Carvalho L. and Ferreira R. (2009). Alterations in autonomic response
head-up tilt testing in paroxysmal atrial fibrillation patients: a wavelet analysis.
Revista portuguesa de cardiologia : orgao oficial da Sociedade Portuguesa de
Cardiologia = Portuguese journal of cardiology : an official journal of the
Portuguese Society of Cardiology 28, 243-257.
Parati G., Mancia G., Di Rienzo M. and Castiglioni P. (2006). Point: cardiovascular
variability is/is not an index of autonomic control of circulation. J Appl Physiol
101, 676-678; discussion 681-672.
Parati G., Pomidossi G., Ramirez A., Cesana B. and Mancia G. (1985). Variability of the
haemodynamic responses to laboratory tests employed in assessment of neural
cardiovascular regulation in man. Clinical science 69, 533-540.
Parati G., Saul J. P., Di Rienzo M. and Mancia G. (1995). Spectral analysis of blood
pressure and heart rate variability in evaluating cardiovascular regulation. A
critical appraisal. Hypertension 25, 1276-1286.
Preiss G. and Polosa C. (1974). Patterns of sympathetic neuron activity associated with
Mayer waves. The American journal of physiology 226, 724-730.
191

Rajan V., Varghese B., van Leeuwen T. G. and Steenbergen W. (2009). Review of
methodological developments in laser Doppler flowmetry. Lasers Med Sci 24,
269-283.
Rauh R., Posfay A. and Muck-Weymann M. (2003). Quantification of inspiratory-
induced vasoconstrictive episodes: a comparison of laser Doppler fluxmetry and
photoplethysmography. Clin Physiol Funct Imaging 23, 344-348.
Raymundo H., Scher A. M., O'Leary D. S. and Sampson P. D. (1989). Cardiovascular
control by arterial and cardiopulmonary baroreceptors in awake dogs with
atrioventricular block. The American journal of physiology 257, H2048-2058.
Rissanen J. (1982). Estimation of structure by minimum description length. Circ Syst Sig
Process 1, 395-406.
Rubinshtein R., Kuvin J. T., Soffler M., Lennon R. J., Lavi S., Nelson R. E., Pumper G.
M., Lerman L. O. and Lerman A. (2010). Assessment of endothelial function by
non-invasive peripheral arterial tonometry predicts late cardiovascular adverse
events. Eur Heart J 31, 1142-1148.
Ruttkay-Nedecký I. (1962). Vasoconstrictor reflex in acral skin regions after taking a
deep breath. Acta Neurovegetativa 24, 413-420.
Sagawa K. (1983). Baroreflex Control of Systemic Arterial Pressure and Vascular Bed.
In Comprehensive Physiology, pp. 453–496.
Sangkatumvong S. (2011). Modeling of Cardiovascular Autonomic Control in Sickle Cell
Disease. Ph.D. Dissertation, University of Southern California.
Sangkatumvong S., Khoo M. C., Kato R., Detterich J. A., Bush A., Keens T. G.,
Meiselman H. J., Wood J. C. and Coates T. D. (2011). Peripheral vasoconstriction
and abnormal parasympathetic response to sighs and transient hypoxia in sickle
cell disease. Am J Respir Crit Care Med 184, 474-481.
Saul J. P., Berger R. D., Albrecht P., Stein S. P., Chen M. H. and Cohen R. J. (1991).
Transfer function analysis of the circulation: unique insights into cardiovascular
regulation. The American journal of physiology 261, H1231-1245.
Scharf S. M., Brown R., Saunders N. and Green L. H. (1980). Hemodynamic effects of
positive-pressure inflation. J Appl Physiol 49, 124-131.
Schnall R. P., Shlitner A., Sheffy J., Kedar R. and Lavie P. (1999). Periodic, profound
peripheral vasoconstriction--a new marker of obstructive sleep apnea. Sleep 22,
939-946.
192

Seals D. R., Suwarno N. O. and Dempsey J. A. (1990). Influence of lung volume on
sympathetic nerve discharge in normal humans. Circulation research 67, 130-
141.
Seals D. R., Suwarno N. O., Joyner M. J., Iber C., Copeland J. G. and Dempsey J. A.
(1993). Respiratory modulation of muscle sympathetic nerve activity in intact and
lung denervated humans. Circulation research 72, 440-454.
Sleight P., La Rovere M. T., Mortara A., Pinna G., Maestri R., Leuzzi S., Bianchini B.,
Tavazzi L. and Bernardi L. (1995). Physiology and pathophysiology of heart rate
and blood pressure variability in humans: is power spectral analysis largely an
index of baroreflex gain? Clinical science 88, 103-109.
Sokolski M., Rydlewska A., Krakowiak B., Biegus J., Zymlinski R., Banasiak W.,
Jankowska E. A. and Ponikowski P. (2011). Comparison of invasive and non-
invasive measurements of haemodynamic parameters in patients with advanced
heart failure. J Cardiovasc Med (Hagerstown) 12, 773-778.
Stern M. D. (1975). In vivo evaluation of microcirculation by coherent light scattering.
Nature 254, 56-58.
Stern M. D., Lappe D. L., Bowen P. D., Chimosky J. E., Holloway G. A., Jr., Keiser H.
R. and Bowman R. L. (1977). Continuous measurement of tissue blood flow by
laser-Doppler spectroscopy. The American journal of physiology 232, H441-448.
Takahashi T., Nishizawa Y., Emoto M., Kawagishi T., Matsumoto N., Ishimura E., Inaba
M., Okuno Y., Shimada H. and Morii H. (1998). Sympathetic function test of
vasoconstrictor changes in foot arteries in diabetic patients. Diabetes care 21,
1495-1501.
Task Force. (1996). Heart rate variability. Standards of measurement, physiological
interpretation, and clinical use. Task Force of the European Society of Cardiology
and the North American Society of Pacing and Electrophysiology. Eur Heart J
17, 354-381.
Theiler J., Eubank S., Longtin A., Galdrikian B. and Farmer J. D. (1992). Testing for
nonlinearity in time series: the method of surrogate data. Physica D 58, 77-94.
Toska K. and Eriksen M. (1993). Respiration-synchronous fluctuations in stroke volume,
heart rate and arterial pressure in humans. J Physiol 472, 501-512.
Toska K., Eriksen M. and Walloe L. (1996). Short-term control of cardiovascular
function: estimation of control parameters in healthy humans. The American
journal of physiology 270, H651-660.
193

Truijen J., van Lieshout J. J., Wesselink W. A. and Westerhof B. E. (2012). Noninvasive
continuous hemodynamic monitoring. Journal of clinical monitoring and
computing 26, 267-278.
Ursino M. (1998). Interaction between carotid baroregulation and the pulsating heart: a
mathematical model. The American journal of physiology 275, H1733-1747.
Vallbo A. B., Hagbarth K. E., Torebjork H. E. and Wallin B. G. (1979). Somatosensory,
proprioceptive, and sympathetic activity in human peripheral nerves. Physiol Rev
59, 919-957.
van de Vooren H., Gademan M. G., Swenne C. A., TenVoorde B. J., Schalij M. J. and
Van der Wall E. E. (2007). Baroreflex sensitivity, blood pressure buffering, and
resonance: what are the links? Computer simulation of healthy subjects and heart
failure patients. J Appl Physiol 102, 1348-1356.
Wallin B. G., Esler M., Dorward P., Eisenhofer G., Ferrier C., Westerman R. and
Jennings G. (1992). Simultaneous measurements of cardiac noradrenaline
spillover and sympathetic outflow to skeletal muscle in humans. J Physiol 453,
45-58.
Wesseling K. H., Jansen J. R., Settels J. J. and Schreuder J. J. (1993). Computation of
aortic flow from pressure in humans using a nonlinear, three-element model. J
Appl Physiol (1985) 74, 2566-2573.
Yang C. C. and Kuo T. B. (1999). Assessment of cardiac sympathetic regulation by
respiratory-related arterial pressure variability in the rat. J Physiol 515 ( Pt 3),
887-896.
Zhong Y., Wang H., Ju K. H., Jan K. M. and Chon K. H. (2004). Nonlinear analysis of
the separate contributions of autonomic nervous systems to heart rate variability
using principal dynamic modes. IEEE transactions on bio-medical engineering
51, 255-262.

Abstract (if available)

Abstract The regulation of peripheral vascular resistance (RPV) is believed to be largely sympathetically‐mediated. Thus assessment of RPV control would allow us to infer valuable information regarding sympathetic nervous activity. Variability in RPV is generally attributed to the baroreflex control of total peripheral resistance (TPR). Although it is known that respiration affects sympathetic outflow and deep breaths, akin to sighs, can lead to peripheral vasoconstriction, the respiratory modulation of TPR has been little studied. In the present study, we utilized noninvasive surrogate measures of RPV to examine the two mechanisms that influence its variability: the baroreflex control of peripheral vascular resistance and the respiratory‐peripheral vascular resistance coupling. The first surrogate measure was obtained from peripheral arterial tonometry (PAT). PAT measured the changes in volume at the finger tip, reflecting the vasoconstriction response as the reduction in its signal amplitude. The other surrogate measure was obtained from laser Doppler flowmetry, which monitors microvascular perfusion. The results of this study suggest that RPV fluctuations were directly modulated by respiration rather than through indirect effect of respiratory modulation of arterial blood pressure (ABP). The simulation model developed based on previous literatures pertinent to short‐term blood pressure regulation could not reproduce the sigh‐vasoconstriction response as observed in the experimental data. The minimal modeling approach was employed to estimate this respiratory coupling effect, which would be incorporated into the simulation model. By means of both modeling approaches, we demonstrated that only after the direct respiratory modulation mechanism was added to the simulation model that a similar vasoconstriction response following a sigh could be reproduced. ❧ The linear and nonlinear dynamics as well as the interaction effect involved in the modulation of RPV through changes in ABP and respiration were investigated in obese pediatric subjects exposed to orthostatic stress and subjects with sickle cell disease before and after blood transfusion treatment. In the obese pediatric subject group, we found that the linear gains of both the TPR baroreflex as well as the respiratory coupling mechanisms diminished as a result of orthostatic stress. The reduction in these gains suggests that sympathetic modulation of TPR decreased in spite of a rise in sympathetic tone. Orthostatic stress was found to lead also to a reduction in the strength of the nonlinear behavior in obese pediatric subjects. Subjects with more severe degrees of metabolic syndrome and obstructive sleep apnea syndrome showed larger reduction in nonlinear TPR baroreflex gain. Transfusion therapy in the sickle cell disease subjects led to an increase in nonlinear TPR baroreflex gain as well as the interaction between ABP and respiration. ❧ In conclusion, through a combination of the structured and the minimal modeling approaches, we have developed an extended model of blood pressure variability that incorporates the respiratory modulation effect on RPV. Taking this respiratory modulation effect into account is important for achieving accurate TPR baroreflex estimation. Finally, the system gains derived from the estimated kernels may constitute potentially useful biomarkers of sympathetic nervous system function.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Role of the autonomic nervous system in vasoconstriction responses to mental stress in sickle cell disease: a bioengineering perspective

PDF

Modeling of cardiovascular autonomic control in sickle cell disease

PDF

Time-varying closed-loop modeling of circulatory control in sleep-disordered breathing

PDF

Cardiorespiratory variability in normal and sleep disordered breathing states: a comprehensive model

PDF

Photoplethysmogram-based biomarker for assessing risk of vaso-occlusive crisis in sickle cell disease: machine learning approaches

PDF

Parametric and nonparametric modeling of cerebral hemodynamics

PDF

Principal dynamic mode analysis of cerebral hemodynamics for assisting diagnosis of cerebrovascular and neurodegenerative diseases

PDF

Parametric and non‐parametric modeling of autonomous physiologic systems: applications and multi‐scale modeling of sepsis

PDF

Systems pharmacology for acute myeloid leukemia and feedback control pharmacodynamic models

PDF

Modeling anti-tumoral effects of drug-induced activation of the cell-extrinsic apoptotic pathway

PDF

Autonomic and metabolic effects of obstructive sleep apnea in childhood obesity

PDF

Integrative modeling of autonomic and metabolic interactions in sleep-disordered breathing

PDF

Pharmacokinetic/pharmacodynamic modeling for genetically polymorphic populations

PDF

Modeling of neonatal hemodynamics during PDA closure

PDF

Estimation of cognitive brain activity in sickle cell disease using functional near-infrared spectroscopy and dynamic systems modeling

PDF

Alpha globin and vascular reactivity in humans: loss of alpha globin genes in human subjects is associated with improved nitric oxide-mediated vascular perfusion

PDF

Data-driven modeling of the hippocampus & the design of neurostimulation patterns to abate seizures

PDF

Nonlinear modeling of causal interrelationships in neuronal ensembles: an application to the rat hippocampus

PDF

Investigating the role of muscle physiology and spinal circuitry in sensorimotor control

PDF

Experimental and computational models for seizure prediction

Asset Metadata

Creator Chalacheva, Patjanaporn (author)

Core Title Modeling autonomic peripheral vascular control

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Biomedical Engineering

Publication Date 05/21/2014

Defense Date 03/26/2014

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag minimal model,OAI-PMH Harvest,peripheral vascular resistance,sympathetic nervous activity

Format application/pdf (imt)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Khoo, Michael C.K. (committee chair), Coates, Thomas D. (committee member), Kato, Roberta (committee member), Marmarelis, Vasilis Z. (committee member), Wood, John C. (committee member)

Creator Email chalache@usc.edu,sang.chalacheva@gmail.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c3-415916

Unique identifier UC11296353

Identifier etd-Chalacheva-2521.pdf (filename),usctheses-c3-415916 (legacy record id)

Legacy Identifier etd-Chalacheva-2521.pdf

Dmrecord 415916

Document Type Dissertation

Format application/pdf (imt)

Rights Chalacheva, Patjanaporn

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA