NONLINEAR IDENTIFICATION OF THE TOTAL BAROREFLEX ARC By  Mohsen Moslehpour   A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of  Electrical Engineering - Doctor of Philosophy 2016     ABSTRACT NONLINEAR IDENTIFICATION OF THE TOTAL BAROREFLEX ARC By Mohsen Moslehpour The baroreflex is one of the most important regulatory mechanisms of blood pressure in the body, and the total baroreflex arc is defined to be the open-loop system relating carotid sinus pressure (CSP) to arterial pressure (AP).  This system is known to exhibit nonlinear behaviors.  However, few studies have quantitatively characterized its nonlinear dynamics.  The aim of this thesis was to develop a nonlinear model of the sympathetically-mediated total arc without assuming any model form in both healthy and hypertensive rats.   Normal and spontaneously hypertensive rats were studied under anesthesia.  The vagal and aortic depressor nerves were sectioned, the carotid sinus regions were isolated and attached to a servo-controlled piston pump.  CSP was perturbed using a Gaussian white noise signal.  A second-order Volterra model was developed by applying nonparametric identification to the measurements.  The second-order kernel was mainly diagonal.  Hence, a reduced second-order model was similarly developed comprising a linear dynamic system in parallel with a squaring system in cascade -43% better than conventional linear dynamic in response to new Gaussian white noise   CSP.  The model also predicted nonlinear behaviors including thresholding and mean responses to CSP changes about the mean.  The linear and nonlinear terms of the validated model between normotensive and hypertensive rats were compared.  While, the linear gain was similar between these two groups, the nonlinear gains for the hypertensive rats were significantly larger.  Hence, nonlinear dynamic functioning of the sympathetically-mediated total arc may enhance baroreflex buffering of AP increases more in spontaneously hypertensive rats than normotensive rats. The importance of higher-order nonlinear dynamics was also assessed via development and evaluation of a third-order nonlinear model of the total arc.  Third-order Volterra and Uryson models were developed by employing several nonparametric and parametric identification methods.  The R2 values between the measured AP and AP predicted by both the best third-order Volterra and the third-order Uryson model in response to new Gaussian white noise CSP were not statistically different from the corresponding values for the previously established second-order Uryson model neither in normotensive nor in hypertensive rats.  Further, none of the third-order models were able to predict important nonlinear behaviors better than the second-order Uryson model.  Additional experiments suggested that the unexplained AP variance was partly due to higher brain center activity.  In conclusion, the second-order Uryson model sufficed to represent the sympathetically-mediated total arc under the employed experimental conditions and the nonlinear part of this model showed significant changes in hypertensive rats compared to normotensive rats.  iv  Dedicated to my loving wife, Sanaz,  and my wonderful parents, Robabeh and Ali Akbar  v  ACKNOWLEDGEMENTS I would like to thank my advisor Dr. Ramakrishna Mukkamala for his continuous guidance and support since I joined Michigan State University at May 2010.  I will be forever grateful for the opportunities he provided me in my career.  I am also thankful for the trust he deposited in my work, for the motivation demonstrated along this arduous course, and for teaching the way of being an independent researcher.     I would also like to express appreciation to my dissertation committee members, Dr. Lalita Udpa, Dr. Selin Aviyente, and Dr. Gregory Fink who were more than generous with their expertise and precious time. The realization of this work was only possible due to the several people's collaboration, to which I desire to express my gratefulness.  I would like to thank Dr. Toru Kawada, Dr. Masaru Sugimachi, and Dr. Kenji Sunagawa for providing the data as well as their efforts and their help in discussion and interpreting the results presented in this thesis. I would also like to thank all my colleagues and friends for the all the discussions and great times that I had during my PhD studies at Michigan State University, Dr. Guanqun Zhang, Mingwu Gao, Jiankun Liu, Anand Chandrasekhar, and Keerthana Natarajan.  Finally, I would like to express my special thanks to my lovely wife, Sanaz Behbahani, for her unconditional love and support who has stood by me through the long nights and hard days during this journey.  I am also thankful to my parents, Robabeh Rasoolzadeh and Ali Akbar Moslehpour, my brother and my sister, Mehdi vi  and Mahnaz, for their everlasting supports.  I would never have been able to accomplish what I have without them.  vii  TABLE OF CONTENTS LIST OF TABLES ............................................................................................................. ix LIST OF FIGURES ............................................................................................................ x CHAPTER 1. INTRODUCTION ................................................................................... 1 Background ..................................................................................................................... 1 Arterial Baroreflex .......................................................................................................... 3 Quantifying Reflexes ....................................................................................................... 5 Scope and Organization .................................................................................................. 6 CHAPTER 2. DATA COLLECTION AND PRE-PROCESSING ................................. 8 Data Collection ................................................................................................................ 8 Data Pre-Processing ...................................................................................................... 10 CHAPTER 3. NONLINEAR IDENTIFICATIOND OF THE TOTAL BAROREFLEX ARC USING SECOND ORDER VOLTERRA MODEL ................................................ 11 Nonlinear Model and Identification Technique ............................................................ 13 Model Evaluation ...................................................................................................... 17 Statistical Comparison ............................................................................................... 17 Results ........................................................................................................................... 18 Gaussian White Noise Data ....................................................................................... 18 Total Arc .................................................................................................................... 19 Neural and Peripheral Arcs ....................................................................................... 25 Discussion ..................................................................................................................... 28 Gaussian White Noise Approach for Nonlinear System Identification..................... 29 Nonparametric Identification Method ....................................................................... 29 Total Arc Model ........................................................................................................ 31 Neural Arc and Peripheral Arc Models ..................................................................... 34 Study Limitations ...................................................................................................... 36 CHAPTER 4. NONLINEAR IDENTIFICATION OF THE TOTAL BAROREFLEX ARC: CHRONIC HYPERTENSION MODEL ................................................................ 38 Nonlinear Model and Estimation Method ..................................................................... 39 Model Testing ............................................................................................................ 40 Model Comparison to WKY ..................................................................................... 41 Results ........................................................................................................................... 41 Gaussian White Noise Data ....................................................................................... 41 Total Arc Model ........................................................................................................ 44 Neural Arc and Peripheral Arc Models ..................................................................... 48 Discussion ..................................................................................................................... 51 Nonlinear Identification Method ............................................................................... 51 Total Arc Model in SHR ........................................................................................... 52 Neural Arc and Peripheral Arc Models in SHR ........................................................ 54 viii  Potential Physiologic Mechanisms ............................................................................ 55 Study Limitations ...................................................................................................... 56 CHAPTER 5. NONLINEAR IDENTIFICATIOND OF THE TOTAL BAROREFLEX ARC: HIGHER-ORDER NONLINEARITY ................................................................... 57 Nonlinear Model and Estimation Method ..................................................................... 57 Model Evaluation ...................................................................................................... 63 Results ........................................................................................................................... 63 Discussion ..................................................................................................................... 67 CHAPTER 6. CONCLUSIONS.................................................................................... 73 APPENDIX ................................................................................................................ 76 BIBLIOGRAPHY ............................................................................................................. 78  ix  LIST OF TABLES Table 1: Group average (meanSE) of the mean and standard deviation of the pre-processed variables during Gaussian white noise CSP stimulation in the training data. .. 18 Table 2: Group average of the R2 values between AP predicted by three models of the total arc and measured AP during Gaussian white noise CSP stimulation. ...................... 23 Table 3: Group average of the R2 values between SNA/AP predicted by models of the neural and peripheral arcs and measured SNA/AP during Gaussian white noise CSP stimulation......................................................................................................................... 27 Table 4: Group average R2 values between arterial pressure (AP) predicted by models of the total baroreflex arc and measured AP. ........................................................................ 46 Table 5: Group average parameters of the kernels of the validated Uryson model of the total arc.............................................................................................................................. 47 Table 6: Group average R2 values between efferent sympathetic nerve activity (SNA)/AP predicted by models of the neural and peripheral arcs and measured SNA/AP. .............. 50 Table 7:  Group average R2 values between arterial pressure (AP) predicted by models of the total baroreflex arc and measured AP. ........................................................................ 67 Table 8: Kernel values of Uryson model of the total baroreflex arc in WKY, SHR120, and SHR160 ............................................................................................................................... 77  x  LIST OF FIGURES Figure 1: Autonomic and hormonal control of cardiovascular function (Adapted from [5])............................................................................................................................................. 2 Figure 2: The total baroreflex arc (open-loop system relating carotid sinus pressure (CSP) to arterial pressure (AP)) has been shown to exhibit nonlinear behaviours in previous studies including mean responses to input changes about the mean (adapted from [23]).  CSP here was controlled using a binary white noise signal of the same mean but increasing amplitude.  Mean AP and sympathetic nerve activity (SNA) both decreased with the increasing CSP amplitude. .................................................................................. 12 Figure 3: Total baroreflex arc can play a role on genesis of Hypertension [adapted from [3]]..................................................................................................................................... 12 Figure 4: Gaussian white noise training data from one subject.  The left plots show pre-processed CSP, AP, and calibrated SNA versus time, while the right plots illustrate the power spectrum and histogram of un-processed CSP.  CSP here was controlled using a Gaussian white noise signal of mean of 120 mmHg and standard deviation of 20 mmHg............................................................................................................................................ 19 Figure 5: Group average of first- and second-order kernel estimates of a Volterra model of the total arc.  The second-order kernel was approximately diagonal, and the diagonal differed in shape from the first-order kernel.  The solid and dashed lines here and in subsequent figures respectively represent mean and mean±SE over the ten subjects for study. ................................................................................................................................. 20 Figure 6: Reduced, second-order Uryson model of the total arc derived by examination of the Volterra kernel estimates in Figure 5.  In this model,  is the input,  is the output, and , and  are the zeroth-, first-, and second-order kernels of the system.  The model is a parallel connection of a linear dynamic system characterized by the first-order kernel and a Hammerstein system (squarer followed by a linear dynamic system characterized by the second-order kernel). ........................................................... 21 Figure 7: Group average of first- and second-order kernel estimates of the reduced Uryson model of the total arc.  Both kernels show integral or low-pass characteristics, but the second-order kernel is more sluggish than the first-order kernel. ............................... 22 Figure 8: Group average of squared coherence functions for the Uryson model and a standard linear model of the total arc in the testing data.  The improved AP prediction offered by the Uryson model was in the low frequency regime.  Hence, the system nonlinearity was at low frequencies. ................................................................................. 23 Figure 9: Group average of the static behavior of the total arc predicted by the Uryson model in response to the staircase CSP (wherein each CSP step or level was flat) and the xi  measured static behavior.  The Uryson model was able to predict thresholding but not saturation. .......................................................................................................................... 24 Figure 10: Predicted AP by the group average Uryson model of the total arc in response to the same binary white noise CSP signal from a previous study shown in Figure 3.  Like the measured AP in Figure 3, the model predicted reductions in mean AP with increasing CSP amplitude.  But, unlike this measured AP, the model predicted increases in AP variance as the CSP amplitude increased.......................................................................... 25 Figure 11: Group average of first- and second-order kernel estimates of a reduced, second-order Uryson model of the neural arc.  Both kernels show derivative or high-pass characteristics.  The kernels are much faster than those of the total arc (see Figure 7). .. 26 Figure 12: Group average of the kernel estimate of a linear model of the peripheral arc.  The kernel shows integral or low-pass characteristics and is similar in speed to those of the total arc (see Figure 7). ............................................................................................... 27 Figure 13:  Gaussian white noise training data from one subject.  SHR120 and SHR160 are spontaneously hypertensive rats during Gaussian white noise carotid sinus pressure (CSP) stimulation with mean of 120 and 160 mmHg, respectively. ........................................... 43 Figure 14:  Group average first-order (linear) and second-order kernel estimates (meanSE) of complete Volterra models (see Equation) of the total baroreflex arc in SHR120 and SHR160.  The inputs for the first-order and second-order kernels are CSP and CSP2, respectively, while the output for both kernels is AP.  Hence, in discrete-time, the units are mmHg/mmHg (unitless) for the first-order kernel and mmHg/mmHg2 (mmHg-1) for the second-order kernel.  Front view precisely means the view point with azimuth of 135° and elevation of 0° with respect to the axis origin. .................................................. 45 Figure 15:  Group average first- and second-order kernel estimates of reduced Uryson models of the total arc in SHR120 and SHR160. ................................................................. 46 Figure 16:  Group average first- and second-order kernel estimates of Uryson models of the neural arc in SHR120 and SHR160. ............................................................................... 49 Figure 17:  Group average first-order kernel estimates of linear models of the peripheral arc in SHR120 and SHR160. ................................................................................................ 49 Figure 18: Group average two-dimensional (2D) slices of the estimated (four-dimensional) third-order kernels ) in Eq. (5).  WKY, Wistar Kyoto rats during Gaussian white noise carotid sinus pressure (CSP) stimulation with mean of 120 mmHg; and SHR120 and SHR160, spontaneously hypertensive rats during the same CSP stimulation with mean of 120 and 160 mmHg, respectively.  The diagonal slice () was largest in magnitude. ......................................................................................................... 64 Figure 19: Group average energy of each 2D slice of the estimated second-order Volterra kernels from [58], [66] and third-order kernels, normalized by the diagonal energy, in descending order of value for WKY, SHR120, and SHR160. The second-order kernel was xii  always diagonal, while the third-order kernel was virtually diagonal for SHR120 and approximately diagonal for WKY and SHR160. ................................................................ 65 Figure 20: Group average first-, second-, and third-order kernel estimates (meanSE) of reduced Uryson models of the total arc for WKY, SHR120 and SHR160 via the frequency-domain (black) and Laguerre expansion (gray) methods for Gaussian inputs.  The three kernels were not proportional to each other, thereby indicating that the model could not be further reduced. ............................................................................................................ 65 Figure 21: Representative time series for arterial pressure (AP) and sympathetic nerve activity (SNA; measured from splanchnic nerve and then normalized as described elsewhere [58]) in response to fixed and Gaussian white noise CSP stimulations for one WKY (left panel).  A.u. is arbitrary units.  Group average power spectra (mean±SE) for CSP, AP, and SNA (right panel).  The gray and black lines indicate the fixed and Gaussian white noise CSP inputs, respectively. ............................................................... 72 1  CHAPTER 1. INTRODUCTION Background Cardiovascular system consists of the heart, blood, and the vessels and it is responsible for maintaining blood flow to the body tissues.  Heart and blood vessels are controlled in particular to provide the required cardiac output (CO) and arterial blood pressure (ABP).  The extrinsic or global control of cardiovascular system is almost entirely through the autonomic nervous system (ANS).  In fact, the main function of cardiovascular regulation is to maintain ABP within a narrow range under a wide range of situations.   Different sensory mechanisms are working together in a feedback loop to regulate arterial blood pressure through ANS.  These sensory monitoring mechanisms entail blood pressure sensors (arterial baroreceptors), blood volume sensors (cardiopulmonary baroreceptors or volume receptors), blood chemistry (chemoreceptors), and plasma osmolarity (osmoreceptors).  These receptors can be categorized in general two categories of mechanical (barosensory) and chemical (chemosensory) [1], [2].  Effects of these sensory mechanisms are called baroreflex and chemoreflex respectively.  Based on the inputs from these sensors, ANS regulates the blood pressure mostly via sympathetic nervous system (SNS) and less through parasympathetic nervous system (PNS).  SNS affects regulation by innervating the blood vessels and the heart.  It passes through two main routes toward the circulation, through specific sympathetic nerves that innervate mainly the vasculature of the 2  internal viscera and the heart, and also through peripheral portions of the spinal nerves where it distributes to the vasculature of the peripheral areas [1].  It affects heart activity by both innervating the sinoatrial node (SA node) to increase the heart rate and also innervating cardiac muscles to elevate the VC.  In addition to heart, SNS innervates small arteries and arterioles by increasing their resistance to blood flow and thereby increasing blood pressure.  It also innervates large vessels (e.g. veins) by decreasing their volume, thus translocating more blood to heart.  PNS plays a minor role in regulation of the circulation by mainly adjusting the heart function through innervation of the SA node.  PNS fibers are carried to the heart via the vagus nerves.  Principally, stimulation of PNS causes marked decrease in HR.  This has been illustrated in Figure 1 (adapted form [5]) along with other hormonal loops which are not in the scope of this thesis.  Figure 1: Autonomic and hormonal control of cardiovascular function (Adapted from [5])  3  Arterial Baroreflex The focus of this thesis is on the most important blood pressure regulator, i.e. the arterial baroreflex mechanism (i.e. regulation through arterial baroreceptors).  The arterial baroreflex system is primarily responsible for maintaining blood pressure in the short-term (seconds to minutes) and also appear to contribute to longer-term blood pressure regulation [3], [4], [14].  It is well known that the arterial baroreflex buffers changes in ABP via stretch receptors belonging to spray type ending nerves lying in the carotid sinus and aortic arch (i.e. carotid sinus baroreflex and aortic arch baroreflex) and maintains the blood pressure near its normal operating level by providing negative feedback to central nervous system.  Fluctuations in blood pressure cause changes in firing patterns of the arterial baroreceptors.  This signal is conveyed to the medulla oblongata within the brainstem (cardiac center in the autonomic nervous system) via afferent nerve fibers.  This signal is deciphered and compared to a set-point for arterial pressure and the proper command is send though efferent fibers.   When ABP increases (decreases), the arterial baroreceptor stretches more (less).  This stretch increases (decreases) the signal going to nervous system through afferent nerves.  Baroreceptors in carotid sinus are transmitted to brain through are transmitted to brain via afferent vagus nerves.  ANS responds to this reflex signal by decreasing (increasing) -sympathetic nerve fibers and reverse effect of efferent -sympathetic nerve fibers also -sympathetic nerve fibers 4  decrease (increase) total peripheral resistance (TPR) by affecting arterioles and increase (decrease) system venous unstressed volume (SVUV) by affecting veins.  These changes, decrease (increase) ABP as a result (Figure 1).  For instance, in postural change and when standing from supine position, due to blood pressure drop in head and upper body, the arterial baroreceptors initiate a reflex so that a reaction of ANS in entire body mediates this blood pressure drop in head and prevents fainting.   This system is negative feedback system.  Hence, it works in closed loop system and this makes it difficult to identify the system ([4]).  Assuming that carotid sinus baroreflex has the same effects as aortic arch baroreflex, one way to handle this problem is to isolate the carotid sinus region from the systemic circulation which causes an open loop preparation ([5][7]) and the pressure in that region is changed independently to identify the open loop system in anesthetized animals.  It should be noted that opening the arterial baroreflex loop can abolish Mayer waves and make the identification more accurate [8].  In this open loop preparation, the baroreflex system is divided into a controller and effector sub-systems, neural arc and peripheral arc respectively [19].  Neural arc represents the relationship between carotid sinus pressure (CSP) and efferent sympathetic nervous activity (SNA).  Peripheral arc is the system from SNA to arterial blood pressure (ABP).  In the open loop regime the cascade of these two systems is called total baroreflex arc  defined to be the open-loop system relating carotid sinus pressure (CSP) to arterial pressure (AP).   Linear dynamic analysis has been widely applied to identify these systems [19], [15], [20], [21].  This type of analysis results in linear transfer function of the system and not only it can elucidates the physiological roles of these arcs but also it 5  helps designing artificial vasomoter center for replacing natural neural arc [22].  Our collaborators others have identified the linear dynamics of the three baroreflex arcs in the form of transfer functions (i.e., gain and phase as a function of frequency) [15][18].  These linear models can capture the dynamic behavior of the systems to a significant extent.  They previously showed that the linear dynamics of the total arc are preserved in spontaneously hypertensive rats (SHR) despite resetting of mean AP [15].  However, the nonlinear dynamics of this system, which have been less investigated, could possibly respond differently to the chronic hypertension model.  Some recent studies showed that arterial baroreflex is not merely a linear system and there are nonlinear behavior with respect to that [20], [23].  Besides, the hypothesis that arterial baroreceptors can also play a role in genesis of hypertension reinforce nonlinear modeling of total arc system [3], [24], [25].   Quantifying Reflexes There are numerous diseases which are affecting neural cardiovascular regulatory reflexes and their functioning.  Therefore, quantifying these system is important in order to understand their functions in health and disease.  Besides, bionic devices can be designed to function as a replacement of neural regulatory mechanism in case of failure in the natural system.   The conventional approach for quantifying ANS regulation of blood pressure is to use an external stimulus (e.g. controlled carotid sinus pressure) to perturb the system and measure the regulatory response.  There are two general types of stimulus, selective and nonselective [27].  Selective stimulus (e.g. controlled carotid sinus pressure) helps quantifying the open loop system and change the normal operation, 6  therefore its interpretation should be done cautiously.  However, it helps characterizing one mechanism alone.  On the other hand, nonselective stimulus (e.g. supine and standing position) is done in a normal closed loop system.  Therefore, all mechanisms are functional and affect each other and since the physiological signals are usually highly correlated, these stimuli may not the reveal the mechanism under study perfectly.  Furthermore, the conventional methods could only quantify the steady-state gain values of the feedback mechanisms while their dynamic characteristics such as overall time constants and delays are undiscovered.  Here, we used the first approach and designed a specific complex experiment along with utilizing signal processing and system identification techniques to identify the system perfectly.  Briefly, we stimulated carotid sinus pressure independently while eliminating other reflexes and measured the regulated arterial blood pressure in an open loop preparation.  Then, the nonlinear model was developed to identify arterial baroreflex dynamics. Scope and Organization  The main goal of this thesis is to use system identification techniques reinforced with signal processing in order to quantitatively characterize the arterial baroreceptor reflexes and investigate the nonlinear dynamic model of the total baroreflex arc under open loop conditions. There are five chapters here. The current chapter (first chapter) gave an introduction and background for the total baroreflex arc.   In chapter two, we developed the second order Volterra and Uryson models of the total baroreflex arc in normotensive Wistar Kyoto rats using non-parametric identification method.  The goal of chapter three is to illustrate the importance of the 7  nonlinearity in chronic hypertension.  Hence, the second order Volterra and Uryson models were estimated for spontaneously hypertensive rats using the same identification method as in normotensive rats and the estimated model was compared with the estimated model for normotensive rats developed in chapter two.  In chapter four, we assessed the importance of higher-order nonlinear dynamics via development and evaluation of a third-order nonlinear model of the total arc using the same experimental data.  Third-order Volterra and Uryson models were developed by employing several nonparametric and parametric identification methods.  Finally, in chapter five perspective and significance of the work is explained.   8  CHAPTER 2. DATA COLLECTION AND PRE-PROCESSING Data Collection In this thesis we used the data that has been already collected by our Japanese collaborators.  Animals were studied according to a protocol that was approved by the Animal Subjects Committee at the National Cerebral and Cardiovascular Center of Japan.  The procedures are described in detail elsewhere [15].  Briefly, under general anesthesia (urethane and -chloralose mixture) and mechanical ventilation, the bilateral vagal and aortic depressor nerves were sectioned to eliminate confounding reflexes from the aortic arch and cardiopulmonary region.  (Hence, the model of the total arc developed herein precisely represents the sympathetically-mediated carotid sinus baroreflex.)  The carotid sinus regions were isolated from the systemic circulation to open the loop between CSP and AP/SNA.  A servo-controlled piston pump was interfaced to the carotid sinus regions filled with warmed Ringer solution via catheters to control CSP.  A femoral artery catheter was placed to measure AP.  A pair of electrodes was positioned on a postganglionic branch of the splanchnic sympathetic nerve to measure SNA.  The pre-amplified SNA was band-pass filtered with cutoff frequencies of 150 and 1,000 Hz.  It was then full-wave rectified and low-pass filtered with cutoff frequency of 30 Hz.  CSP was controlled for about 15 min using a Gaussian white noise signal with mean of 120 mmHg and standard deviation of 20 mmHg (with values more than three standard deviations from the mean being skipped).  So, the signal ranged from about 90 to 150 mmHg for 90% of the 9  stimulation.  A different realization of this signal was employed for each subject.  The switching interval of the noise was 0.5 sec to yield relatively flat input spectral power up to 1 Hz (see Figure 4).  To investigate static system behavior, CSP was also controlled using a staircase signal that started at 60 mmHg and then increased, step by step, in increments of 20 mmHg every 1 min up to 180 mmHg.  So, for example, CSP was held flat at 100 mmHg in the third step of this signal.  Thirteen normotensive Wistar-Kyoto rats (weight 397.8±18.5 grams) were studied according to this protocol.  All measurements were recorded at a sampling rate of at least 200 Hz (In half of the subjects, data were sampled at the rate of 1kHz and in the other half, they were sampled at the rate of 200Hz).  We refer to this stimulation as WKY.   To study the hypertension, another set of experiments were performed in eight SHR (22.2±4.5 weeks in age) under the same protocol.  CSP was stimulated with two different means but the same otherwise as described above.  First, to establish models for SHR at the same CSP level as the previous models for WKY, the Gaussian white noise CSP stimulation was at mean of 120 mmHg and standard deviation of 20 mmHg for about 15 min.  Again, the switching interval of the noise was 500 ms to produce relatively flat CSP spectral power up to 1 Hz (see Figure 13).  Second, to establish models for SHR at the normal CSP level of SHR, the Gaussian white noise CSP stimulation was at mean of 160 mmHg but the same otherwise.  Hereafter, we refer to the former stimulation as SHR120 and the latter stimulation as SHR160.  All signals were continuously recorded at a sampling rate of 200 Hz while CSP was controlled using two different Gaussian white noise stimulations.    10  Data Pre-Processing The measurements during Gaussian white noise stimulation were first low-pass filtered using a high-order filter and then down-sampled to 2 Hz.  For each subject, a 6-min segment of stationary data after linear de-trending was selected for model development or training, while a separate 3-min segment of stationary data after linear de-trending was selected for model testing.  Data from three of the WKY subjects, one SHR120 subject, and three SHR160 subjects were highly non-stationary and were thus excluded from further analysis.  In some of the subjects, a peak in the AP power spectrum around 0.7 to 0.8 Hz was visible.  This peak was likely caused by spontaneous respiratory effort rather than the CSP stimulation, so AP of all subjects was low-pass filtered again using a high-order filter but with a cutoff frequency of 0.7 Hz.  This filtering had no significant impact on the kernel estimates of those subjects that did not reveal such a peak (results not shown).  Finally, since the magnitude of the SNA measurement heavily depended on the electrode contact, SNA was calibrated per subject so that the models of the neural and peripheral arcs could be meaningfully averaged over the subjects.  In particular, SNA was calibrated per subject such that the average gain of the linear kernel of the neural arc was unity for frequencies < 0.03 Hz in the training data [15].  For the SHR subjects, this calibration was done based on the SHR120 training data. The measurements during the staircase stimulation were averaged over the last 10 sec of each step.  The average values of a system input and output were then plotted against each other. 11  CHAPTER 3. NONLINEAR IDENTIFICATIOND OF THE TOTAL BAROREFLEX ARC USING SECOND ORDER VOLTERRA MODEL  The total baroreflex arc  defined to be the open-loop system relating carotid sinus pressure (CSP) to arterial pressure (AP)  is a well-known contributor to cardiovascular regulation.  When stimulated in a controlled manner, this system exhibits thresholding and saturation (i.e., maximal and minimal AP responses) [20], [30], [31], mean responses (i.e., direct current or DC responses) to input changes about the mean (i.e., alternating current or AC changes) [23], [30], [32], [33], as shown in Figure 2, and other nonlinear behaviors [34][36].  Interestingly, the system onic baroreceptor unloading model of hypertension [3] as shown in Figure 3.  This model does not significantly alter mean CSP but does cause reductions in CS pulse pressure (i.e., a selective AC change), which, in turn, leads to a sustained, baroreflex-mediated increase in mean AP (i.e., a DC response).  Hence, total arc nonlinearity could possibly play a role in the genesis of hypertension.  Yet, most system identification studies of the total arc have been based on a linear dynamic model [19], [15], [20], [33], [37], [16][18].  Further, the few studies that have represented the total arc with a nonlinear dynamic model assumed a particular form for the model [20], [35].  12   Figure 2: The total baroreflex arc (open-loop system relating carotid sinus pressure (CSP) to arterial pressure (AP)) has been shown to exhibit nonlinear behaviours in previous studies including mean responses to input changes about the mean (adapted from [23]).  CSP here was controlled using a binary white noise signal of the same mean but increasing amplitude.  Mean AP and sympathetic nerve activity (SNA) both decreased with the increasing CSP amplitude. This latter behavior could play a role in the genesis of hypertension, as indicated by chronic unilateral baroreceptor unloading study [3] shown in Figure 3.  In this figure, it is shown that while CSP is kept constant while baroreceptors are unloaded, mean AP and CSP pulsatility increases in coarse of days.     Figure 3: Total baroreflex arc can play a role on genesis of Hypertension [adapted from [3]] 13  Our aim was to establish a second-order, nonlinear dynamic model of the total arc without making a priori assumptions about the model form.  To achieve this aim, we employed the powerful Gaussian white noise approach for nonlinear system identification [38].  In particular, we applied Gaussian white noise CSP stimulation, while measuring AP and sympathetic nerve activity (SNA), in an open-loop rat preparation followed by nonparametric identification to estimate first- and second-order kernels of a Volterra model of the total arc from the measurements.  We also likewise identified two sub-systems of the total arc, namely the neural arc, which relates CSP to SNA, and the peripheral arc, which relates SNA to AP.  Since this approach requires long data records to yield accurate kernel estimates [39], our strategy for the obtained short data records was as follows.  First, we examined the Volterra kernel estimates to define a reduced second-order, nonlinear dynamic model.  Then, we applied nonparametric identification to estimate the kernels of this reduced, yet potentially more predictive, model.  Finally, we assessed its output predictions.  The nonlinear model of the total arc that we report here significantly improved upon AP predictions over a standard linear model and helped reveal the structure of the total arc. Nonlinear Model and Identification Technique In general, a time-invariant system with fading memory can be written in the form of a Volterra series to within arbitrary precision [40].  For such systems that are also causal and discrete-time, the Volterra series is given as follows:            (1) 14  where  is discrete-time,  is the output, is the input,  is the lth order system kernel, L is the order of nonlinearity, and M is the system memory (which can be different for each kernel but is the same for all kernels here for convenience).  In this model, the output is expanded in terms of the input samples and the interactions amongst input samples of different lags.  These input terms affect the output through the kernels of the system.     In this study, the total arc and its sub-systems were assumed to be represented by a second-order Volterra series as follows:      (2) where  and  are the system input and output (i.e., CSP and AP for the total arc, CSP and SNA for the neural arc, and SNA and AP for the peripheral arc) with  precisely denoting the input after removing its mean value.  The zeroth-order kernel , which is the mean value of , along with the mean value of the input define the system operating point.  The first-order or linear kernel  indicates how the present and past input samples affect the present output sample.  The second-order kernel  indicates how the interaction or cross-talk between two input samples that are  and  samples in the past affect the present output sample.   The kernels of the total, neural, and peripheral arcs were estimated from the Gaussian white noise training data using a nonparametric, frequency-domain method [41].  This method was more effective than other nonparametric methods (see Discussion).  The memory M was set to 25 sec, which is twice the length of the linear 15  kernel of the total arc reported in previous studies [15].  This value was able to capture the memory of all systems (see Results).  The second-order kernel estimates were then visually examined to ultimately arrive at reduced, yet potentially more predictive, nonlinear models (see Results).  Note that nonparametric identification was employed, because it does not impose a particular form for the kernels.  However, the trade-off is that long data records are needed to accurately estimate higher-order kernels.  Since the training data here were relatively short, the Volterra series had to be limited to second-order.  However, this limitation may not be too serious, as many physiologic systems can be well represented with a second-order Volterra series [38].  The kernels of the Volterra model in Eq. (2), , , and , were estimated from the measured zero-mean input  and measured output  using a frequency-domain method.  This method is described in detail elsewhere [41].  Briefly, the kernels were estimated in succession.  First, the zeroth-order kernel was estimated as the mean of the output as follows:  where  is the expectation operator.  Next, the first-order kernel was estimated by first subtracting the contribution of the zeroth-order kernel from the output and then computing the cross-spectrum divided by the input spectrum (i.e., Wiener filter) as follows:   16  where  is the auto- or cross-correlation function between the indicated signals,  is the convolution operator, and  is the Fourier Transform operator.  Finally, the second-order kernel was estimated by first subtracting the contribution of the zeroth- and first-order kernels from the output and then computing a two-dimensional generalization of the Wiener filter as follows:     where  is again the correlation function amongst the indicated signals,  is the two-dimensional convolution operator, and  and  are one- and two-dimensional Fourier Transform operators, respectively.  Note that  and  above were computed via the standard sample mean and unbiased correlation function estimates. The kernels of the Uryson model of Figure 6, , , and , were estimated analogously.  First, the zeroth- and first-order kernels were estimated, as described above.  Then, the contribution of these kernels was subtracted from the output, also as described above.  Finally, the second-order kernel was estimated by first squaring the input and then computing the Wiener filter as follows:      17  While these steps actually yield the kernels of a Wiener model, for a second-order nonlinear system, the first- and second-order kernels of Volterra and Wiener models are the same [38]. Model Evaluation The merit of the resulting nonlinear models with the first- and second-order kernel estimates and linear models with only the first-order kernel estimates was evaluated as follows.  First, the inputs from the Gaussian white noise training and testing data were applied to the models, and R2 values between the predicted and measured outputs and squared coherence functions (the power spectrum of the predicted output divided by the power spectrum of the measured output) were computed.  Then, the inputs from the staircase data were applied to the models, and the predicted and measured outputs were compared qualitatively.  Finally, binary white noise CSP with mean of 95 mmHg but amplitudes of 5, 10, 20, and 40 mmHg and switching interval of 0.5 sec were applied to the models of the total and neural arcs, and the predicted outputs were qualitatively compared to the corresponding measured data from a previous study shown in Figure 3 [23].    Statistical Comparison The R2 values from linear and nonlinear models were compared using paired t-tests.  Before applying these tests, the R2 values were log transformed for more normally distributed data [42].  A p < 0.0125 was considered statistically significant based on a Bonferroni correction for up to four pairwise comparisons.   18  Results Gaussian White Noise Data Figure 4 shows the pre-processed CSP, AP, and calibrated SNA from the Gaussian white noise training data of one subject.  Table 1 shows the group average (meanSE) of the mean and standard deviation of these measurements and the pulse rate (PR).  The mean of CSP, AP, and SNA were 120.30.2 mmHg, 97.14.4 mmHg, and 80.711.8 au, respectively.  These levels define the operating points of the models of the total arc and its sub-systems developed herein.  The standard deviation of CSP, AP, and SNA were 16.50.3 mmHg, 6.6.0.2 mmHg, and 20.70.6 au, respectively.  These values indicate the range of validity of the models about their operating points.  The mean of PR was 39711 bpm, which corresponds to 6-7 Hz.  Hence, the systems were mainly stimulated at sub-pulsatile frequencies (i.e., frequencies beneath the PR).   Table 1: Group average (meanSE) of the mean and standard deviation of the pre-processed variables during Gaussian white noise CSP stimulation in the training data.  CSP AP SNA PR Mean 120.3±0.2 97.1±4.4 80.7±11.8 397.3±11.1 Standard Deviation 16.5±0.3 6.6±0.2 20.7±0.6 2.8±0.3 CSP is carotid sinus pressure (mmHg); AP, arterial pressure (mmHg); SNA, calibrated sympathetic nerve activity (au), and PR, pulse rate (bpm).    19   Figure 4: Gaussian white noise training data from one subject.  The left plots show pre-processed CSP, AP, and calibrated SNA versus time, while the right plots illustrate the power spectrum and histogram of un-processed CSP.  CSP here was controlled using a Gaussian white noise signal of mean of 120 mmHg and standard deviation of 20 mmHg.    Total Arc Figure 5 shows the group average (meanSE) of the first- and second-order kernels of a Volterra model of the total arc estimated from the Gaussian white noise training data.   20   Figure 5: Group average of first- and second-order kernel estimates of a Volterra model of the total arc.  The second-order kernel was approximately diagonal, and the diagonal differed in shape from the first-order kernel.  The solid and dashed lines here and in subsequent figures respectively represent mean and mean±SE over the ten subjects for study. There are two points to note.  First, the second-order kernel revealed small off-diagonal values.  Indeed, none of the off-diagonal static gains (i.e., sums of ,  were significantly different from zero based on one-sample t-tests, except for the static gain of  (p < 0.01).  However, the static gain of this fourth off-diagonal was less than 25% of that of the main diagonal.  Hence, the second-order kernel was approximately diagonal, thereby indicating that cross-talk between pairs of input samples of different lags (i.e.,  in Eq. (2)) hardly contributed to the output.  Second, the diagonal of the second-order kernel appeared different in shape from the first-order kernel.  Both of these findings were largely consistent for the individual subject kernel estimates 21  (results not shown).  The findings suggested to set  in Eq. (2) to arrive at the following reduced, yet potentially more predictive, second-order nonlinear model:            (3)  Figure 6: Reduced, second-order Uryson model of the total arc derived by examination of the Volterra kernel estimates in Figure 5.  In this model,  is the input,  is the output, and , and  are the zeroth-, first-, and second-order kernels of the system.  The model is a parallel connection of a linear dynamic system characterized by the first-order kernel and a Hammerstein system (squarer followed by a linear dynamic system characterized by the second-order kernel). Figure 6 shows a block diagram of the reduced model, which may be categorized as a Uryson model [43].  This model is a linear dynamic system in parallel with a squarer in cascade with another linear dynamic system.  The kernel of the former system is the linear kernel, whereas the kernel of the latter system is the second-order kernel.  The kernels of this reduced model were re-estimated using a nonparametric, frequency-domain method.  Figure 7 shows the resulting group average kernel estimates.  Both kernels showed low-pass or integral characteristics and similar dynamics.  That is, an impulsive increase in CSP at time zero would cause AP to initially decrease and then return to baseline.  The static gain of the first-order kernel (i.e., change in steady-state output divided by change in steady-state input) was -0.70 (unitless).  The static gain of the second-order kernel, unlike the first-order kernel gain, depends on the size and sign of the input change due to the squaring operation.  The static gain of this kernel was -0.22 (unitless) for an average step CSP 22  increase of 16 mmHg or +0.22 for a CSP decrease of 16 mmHg.  The dominant time constants (computed robustly via the kernel sum divided by the peak kernel amplitude [44]) of the first- and second-order kernels were about 4.1 and 6.2 sec, respectively.   Figure 7: Group average of first- and second-order kernel estimates of the reduced Uryson model of the total arc.  Both kernels show integral or low-pass characteristics, but the second-order kernel is more sluggish than the first-order kernel. Table 2 shows the group average of the R2 values between the AP predicted by the individual subject models when stimulated by the Gaussian white noise CSP in the training and testing data and the measured AP.  The training data results actually reflect model fitting rather than model prediction capabilities.  The Volterra model achieved the best model fit, as indicated by the higher R2 values, simply because it was a superset of the other models.  For the testing data, which do indicate model prediction abilities, the linear model achieved a fairly high R2 value of 0.640.03.  While the Volterra model did not improve upon this value, the Uryson model, which 23  is simpler and may thus include more accurate kernel estimates, attained an R2 value of 0.710.03.  So, the Uryson model improved AP prediction over the linear model by 12% (p < 0.01).   Table 2: Group average of the R2 values between AP predicted by three models of the total arc and measured AP during Gaussian white noise CSP stimulation.  Linear Second-order Volterra Second-order Uryson Training Data 0.73±0.03 0.85±0.01* 0.79±0.03* Testing Data 0.64±0.03 0.64±0.04 0.71±0.03* *denotes p < 0.01 for paired t-test between indicated linear and nonlinear models after log transformation of the R2 values. Figure 8 shows the group average of the squared coherence functions for the linear and Uryson models in the testing data.  As can be seen, the improved AP prediction afforded by the Uryson model was in the low frequency regime.  Hence, the system nonlinearity was at low frequencies.      Figure 8: Group average of squared coherence functions for the Uryson model and a standard linear model of the total arc in the testing data.  The improved AP prediction offered by the Uryson model was in the low frequency regime.  Hence, the system nonlinearity was at low frequencies.  24  Figure 9 shows the group average of the static behavior of the total arc predicted by the individual subject Uryson models in response to the staircase CSP (wherein each CSP step or level was flat) and the measured static behavior.  The model predicted thresholding (qualitatively) but not saturation.    Figure 9: Group average of the static behavior of the total arc predicted by the Uryson model in response to the staircase CSP (wherein each CSP step or level was flat) and the measured static behavior.  The Uryson model was able to predict thresholding but not saturation. Figure 10 shows the AP predicted by the group average Uryson model when stimulated by the binary white noise CSP of increasing amplitude.  Like the measured AP from a previous study [23] shown in Figure 3, the model predicted significant mean AP reductions with increasing amplitude.  But, unlike the measured AP, the model also predicted increases in AP variance as the amplitude increased.  Note that the linear model cannot predict any of these nonlinear behaviors. 25   Figure 10: Predicted AP by the group average Uryson model of the total arc in response to the same binary white noise CSP signal from a previous study shown in Figure 3.  Like the measured AP in Figure 3, the model predicted reductions in mean AP with increasing CSP amplitude.  But, unlike this measured AP, the model predicted increases in AP variance as the CSP amplitude increased.  Neural and Peripheral Arcs  The first- and second-order kernels of a Volterra model of the neural arc estimated from the Gaussian white noise training data also suggested a reduced Uryson model (results not shown).  Figure 11 shows the group average of the first- and second-order kernels of the Uryson model of the neural arc estimated from these data.  Both kernels showed high-pass or derivative characteristics and similar dynamics.  The static gain of the first-order kernel was -0.57 (au/mmHg), while the static gain for the second-order kernel was -0.12 (au/mmHg) for an average step CSP increase of 16 mmHg or +0.12 for a CSP decrease of 16 mmHg.  The dominant time constants of the first and second-order kernels were about 0.3 and 0.6 sec, 26  respectively.  While these small time constants may not have been accurately estimated due to the 2 Hz sampling rate, it is clear that the neural arc was much faster than the total arc.    Figure 11: Group average of first- and second-order kernel estimates of a reduced, second-order Uryson model of the neural arc.  Both kernels show derivative or high-pass characteristics.  The kernels are much faster than those of the total arc (see Figure 7).   Table 3 shows the group average of the R2 values between the SNA predicted by the individual subject models when stimulated by the Gaussian white noise CSP in the training and testing data and the measured SNA.  Again, as expected and indicated by the training data results, the Volterra model achieved the best model fit.  However, the R2 values were only modestly higher than those of the linear and Uryson models here.  For the testing data, the linear model attained a high R2 value of 0.770.02.  The nonlinear models did not significantly improve upon this value.  Hence, the neural arc was approximately linear.  Indeed, even though the neural arc exhibits the 27  nonlinear behaviors of thresholding and saturation [20] and DC responses to AC changes, as indicated via SNA in Figure 3, the nonlinear models of the neural arc showed predictions of these behaviors that were not that different from the linear model (results not shown). Table 3: Group average of the R2 values between SNA/AP predicted by models of the neural and peripheral arcs and measured SNA/AP during Gaussian white noise CSP stimulation.  Neural Arc Peripheral Arc  Linear Second-order Volterra Second-order Uryson Linear Training Data 0.80±0.01 0.84±0.01* 0.82±0.01 0.87±0.02 Testing Data 0.77±0.02 0.80±0.01 0.79±0.02 0.81±0.04 *denotes p < 0.01 for paired t-test between indicated linear and nonlinear models after log transformation of the R2 values.   Figure 12: Group average of the kernel estimate of a linear model of the peripheral arc.  The kernel shows integral or low-pass characteristics and is similar in speed to those of the total arc (see Figure 7). Figure 12 shows the group average of the first-order kernel of a linear model of the peripheral arc estimated from the Gaussian white noise training data.  A reliable second-order kernel could not be estimated, since SNA, which is the input to 28  the peripheral arc, was not Gaussian white noise.  The linear kernel showed low-pass or integral characteristics and the expected, positive open-loop dynamics.  That is, an impulsive increase in SNA at time zero would cause AP to initially increase and then return to baseline.  The static gain was 1.6 mmHg/au.  Its dominant time constant was 5.8 sec, so, as expected, the peripheral arc was more sluggish than the neural arc. Table 3 shows the group average of the R2 values between the AP predicted by the individual subject linear model when stimulated by the SNA in the Gaussian white noise training and testing data and the measured AP.  For both the training and testing data, the R2 value was high.  In particular, for the more meaningful testing data, the R2 value was 0.810.04.  Hence, the peripheral arc was likely approximately linear anyhow.   Discussion We developed a second-order, nonlinear dynamic model of the sympathetically-mediated total baroreflex arc by employing Gaussian white noise stimulation and nonparametric identification.  We validated the model by showing that it could predict (i) AP appreciably better than a standard linear model when stimulated by a new Gaussian white noise realization and (ii) the important nonlinear behaviors of thresholding and DC responses to AC changes.  The validated model illustrates that the form of the second-order nonlinearity of the total arc is close to diagonal.  This result, which represents the major finding herein, means that the square of zero-mean CSP (i.e., ) causes changes in AP (i.e., ), but the --mean CSP at different lags (e.g., 29   or ) have little impact on AP.  To shed light on the sources of nonlinearity, we likewise developed and tested models of the neural and peripheral arcs.  But, the models of these two sub-systems of the total arc showed approximately linear behaviors.     Gaussian White Noise Approach for Nonlinear System Identification     Application of white noise inputs allows systems to be accurately identified over the entire frequency range.  Such inputs are useful for both linear and nonlinear system identification.  However, white noise inputs specifically generated from a Gaussian distribution are needed for nonlinear identification.  In particular, since nonlinear systems do not obey the amplitude scaling property (e.g., if  results from , then  will result from ), a broad range of amplitude excitation is required for reliable identification of nonlinear systems.  Gaussian inputs provide such a range, whereas binary inputs, for example, only provide two amplitude levels.  Furthermore, Gaussian inputs facilitate nonlinear system identification in other ways including orthogonalization of the functionals so that each kernel can be estimated independently.  We refer the reader to [38] for more information on the powerful Gaussian white noise approach. Nonparametric Identification Method We estimated the first- and second-order kernels of the models using short periods (6 min) of Gaussian white noise stimulation.  Hence, the inputs were not exactly Gaussian and white.  We specifically applied a frequency-domain method to estimate the kernels without assuming any form [41].  Other nonparametric 30  -correlation method [45].  This method assumes that the input is strictly Gaussian white noise and thus requires long data records [39]also well known, efficiently solves the normal equations [46], [47].  This method only assumes that the input is broadband.  However, by not assuming a Gaussian input, it must compute higher-order correlations (to form the normal equations), which also requires long data records [39].  Also note that this method can only produce an estimate, if the data length is at least equal to the number of kernel samples for estimation.  The frequency-domain method for estimating the Volterra kernels assumes that the input is Gaussian but broadband and may thus be a good compromise between these methods.  We actually applied all three methods.  While the kernel estimates of the methods were similar on average, the frequency-domain method yielded the smoothest estimates (results not shown).   The frequency-domain method provides the optimal estimates of the linear and nonlinear kernels in the least squares sense without assuming a model form [48].  The resulting kernels of this nonparametric identification method are thus the best unbiased estimates.  In other words, the linear part of the nonlinear model and the best linear model are one in the same.  Parametric identification methods could provide better estimates of both the linear and nonlinear kernels by trading off bias for precision.  However, these methods assume a particular model form.  Our purpose here was not to assume a model form but rather to discover the form.  Nevertheless, we also applied standard autoregressive exogenous input (ARX) identification to estimate the linear kernels [28].  However, this parametric identification method did 31  not yield more predictive linear kernels than the frequency-domain method (results not shown).  Regardless of this finding, comparison between the linear and nonlinear models estimated by the frequency-domain method may be considered fair, since neither model assumes a particular form.          Total Arc Model We arrived at the second-order, nonlinear dynamic model of the total arc shown in Figure 6 and Figure 7 in two steps.  First, we applied the frequency-domain method to estimate the kernels of a Volterra model.  The resulting second-order kernel shown in Figure 5 was approximately diagonal.  Hence, only the past values of the square of the input, rather than the product of input samples of different lags, contributed to the output.  Further, this diagonal differed in shape from the first-order kernel (see Figure 5).  So, the Volterra model may be reduced to an Uryson model (see Figure 6).  Second, we again applied the frequency-domain method with the aim of more accurately estimating the kernels of the reduced model (see Figure 7).  The resulting model is a linear system in parallel with a squarer in cascade with another linear system.  The kernel of the former system () is the linear kernel, while the kernel of the latter system  is the second-order kernel, which was more sluggish.  We tried to further simplify the total arc model by assuming that .  However, the resulting Hammerstein model did not improve AP predictions over a linear model (results not shown).  We also tried to derive a more accurate total arc model by estimating the second-order kernel with a non-zero diagonal and a non-zero fourth off-diagonal (whose gain was statistically different from zero in the 32  Volterra model).  However, this model did not improve AP predictions over the Uryson model (results not shown). We assessed the validity of the Uryson model of the total arc.  In particular, we applied three different CSP inputs to the model, and compared the predicted AP to the measured AP.   First, we applied a Gaussian white noise input that was not utilized to develop the model.  As shown in Table 2, the linear model with only the first-order kernel was able to predict 64% of the measured AP variance, while the Uryson model with both kernels predicted 71% of the variance.  Although the linear model was quite explanatory, the nonlinear model significantly improved the AP prediction by 12% (p < 0.01).  The squared coherence function of Figure 8 indicated that this AP improvement was at low frequencies.  Hence, the system nonlinearity was in the low frequency regime.  While the Uryson model did improve the prediction, 29% of the measured AP variance remained unexplained.  These variations were not white (results not shown) and could be due to higher-order nonlinearity, non-stationarity, SNA from higher brain centers, and fast-acting hormonal loops.  Note that measurement noise may not have been a factor, as AP was invasively measured and then low-pass filtered and down-sampled all the way to 2 Hz.  Also note that low frequency AP may be similar regardless of the site of measurement [49] .   Second, we applied a staircase input to predict static system behavior.  As shown in Figure 9, the model was able to predict thresholding.  The mechanism for this prediction is as follows.  When CSP increases relative to its mean value, the second-order kernel enhances the magnitude of the AP drop.  However, when CSP 33  decreases relative to its mean value, the second-order kernel blunts the AP increase.  On the other hand, the model did not predict saturation. One reason is that the Gaussian white noise stimulation used to develop it (mean and standard deviation of 120 and 16.5 mmHg) hardly excited the saturation regime (CSP > 160 mmHg).  Quantitative differences between predicted and measured thresholding may have been due to differences in the operating points of the data used to develop the model (mean CSP and AP of 120 and 97 mmHg as shown in Table 1) and test the model (mean CSP and AP of 120 and 120 mmHg as shown in Figure 9).  Third, we applied binary white noise of increasing amplitude.  As illustrated in Figure 10, the model predicted reductions in mean AP, but increases in AP variance, with increasing amplitude.  While the corresponding measured AP from a previous study (see Figure 3) [23] likewise indicated mean AP reductions, it revealed little change in AP variance.  Higher-order nonlinearity may be needed to blunt the AP variance.  Another possible reason for the difference between the prediction and measurement may be variations in the range of CSP input used to develop and test the model.   As implied above, the Uryson model developed herein can only be expected to be valid over the range of data utilized in its development.  This range is defined by a CSP of 120±16.5 mmHg (mean±SD) and mostly < 1 Hz and an AP of 976.6 mmHg.  The model should be considered only over this range. A few other nonlinear dynamic models of the total arc have been previously conceived.  After finding that neither a Hammerstein model (a static nonlinearity 34  followed by a linear dynamic system) nor a Wiener model (a linear dynamic system followed by a static nonlinearity) could explain data, a sandwich model  a Wiener model (to represent the neural arc) in cascade with a linear model (to represent the peripheral arc)  was proposed to represent the total arc [20].  A more complicated nonlinear model, which may also be viewed as a sandwich model, was developed earlier [35].  The main difference between these previous efforts and the present study is that we did not assume a certain model form.  Rather, we let the data dictate the form via Gaussian white noise stimulation and nonparametric identification.  Indeed, we found that the total arc could be represented with an Uryson model, which contradicts a sandwich, Hammerstein, and Wiener model.  That is, a sandwich model with a static nonlinearity that is odd about the operating point (similar to a thresholding and saturation curve) would show an identically zero second-order kernel; a Hammerstein model would reveal identical first- and second-order kernels; and a Wiener model would show a second-order kernel with non-zero off-diagonal values [50].   Neural Arc and Peripheral Arc Models We arrived at the second-order Uryson model of the neural arc shown in Figure 11 using a similar two step approach and likewise assessed the model.  However, this model (and a Volterra model) displayed approximately linear behavior.  In particular, when stimulated by a new realization of Gaussian white noise, it could not predict SNA better than a linear model, as shown in Table 3.  Indeed, the linear model could already explain much (77%) of the SNA variance.  The unexplained 35  variations could be due to the aforesaid factors as well as measurement noise.  In addition, while the neural arc shows thresholding and saturation [20] and DC responses to AC changes (see SNA response in Figure 3), the Uryson model of the neural arc could not well predict these behaviors.  We could not develop a reliable nonlinear model of the peripheral arc, because the SNA input was not Gaussian white noise.  We thus settled upon the linear model shown in Figure 12.  When stimulated by new Gaussian white noise, this model was also able to predict much (81%) of the AP variance, as shown in Table 3. In sum, while the total arc exhibited appreciable nonlinear behaviors, its two sub-systems displayed approximately linear behavior.  Hence, identification of the neural and peripheral arcs did not shed light on the sources of total arc nonlinearity.  One possible explanation for this seemingly contradictory finding is that splanchnic SNA, which was used to construct the sub-system models, did not represent whole body SNA, which actually determined total arc behavior.  However, we performed pilot experiments and found similar cardiac and splanchnic SNA responses to controlled CSP stimulation (results not shown).  So, these experiments did not support this explanation.  Another explanation may be that the neural arc is actually nonlinear but its nonlinearity was not well identified.  In particular, as indicated in Figure 8, the nonlinearity was in the low frequency regime.  However, SNA power was predominantly in the high frequency band due to the derivative characteristics of the neural arc.  Hence, the little SNA power at low frequencies may have been dominated by SNA from higher brain centers rather than from the baroreflex.  As described in [38]zing effect in the identification 36  process such that neural arc nonlinearity was masked over the low frequency regime.  Note that afferent C-fibers could be responsible for neural arc nonlinearity.  These fibers are slow acting nerves [51], [52] that are highly nonlinear with respect to stretch [52].  Such behaviors are congruent with the findings here that the nonlinear behavior of the total arc was in the low frequency regime.   As indicated above, a Wiener model of the neural arc was previously proposed [20], [23].  We also tried to represent this system with Wiener and Hammerstein models.  But, neither model appreciably improved SNA predictions over a linear model (at most 4% when stimulated by new Gaussian white noise).  The reason for the difference between this and past studies could possibly be variations in the amplitude of the CSP stimulation employed.  Study Limitations Our study has several limitations.  First, the use of anesthesia and open-loop conditions surely impacted the models.  However, closed-loop identification has its own challenges [53].  Second, CSP excitation was limited to an amplitude range that hardly reached the saturation regime, a short time period, and frequencies mainly within 1 Hz.  Hence, saturation and DC responses to long-term, pulsatile changes could not be modeled.  However, note that increasing the amplitude, time period, and switching rate of CSP stimulation could damage the baroreceptors and cause major non-stationarity. Third, the vagal arm of the total arc was abolished.  However, inclusion of this arm would also bring in the confounding effects of the cardiopulmonary baroreflex.  Finally, the model was restricted to second-order 37  nonlinearity due to the short data records.  However, parametric identification methods, which assume a particular model form, may be required to estimate higher-order kernels from short data records [50].   38  CHAPTER 4. NONLINEAR IDENTIFICATION OF THE TOTAL BAROREFLEX ARC: CHRONIC HYPERTENSION MODEL In the previous chapter, we employed the Gaussian white noise approach for nonlinear system identification to develop a second-order, nonlinear dynamic model of the total arc in normotensive Wistar Kyoto rats (WKY).  The model predicted AP 12% better than a linear dynamic model in response to new Gaussian white noise and important nonlinear behaviors including baroreflex thresholding and mean responses to input changes about the mean.  The validated model revealed that the structure of the total arc is a linear dynamic system in parallel with a cascade combination of a squaring system and a different linear dynamic system, as shown in Figure 6.  This [56].   In this chapter, we aimed to likewise establish second-order, nonlinear dynamic models of the total arc as well as its sub-systems in spontaneously hypertensive rats (SHR) and to compare these models to our previously published models for WKY.  Our results indicate that the second-order nonlinear dynamics of the total arc in SHR, which showed the same structure as in WKY, are augmented significantly more than the linear dynamics.  Hence, nonlinear dynamic functioning of the total arc may enhance baroreflex buffering of AP increases more in SHR than WKY. 39  Nonlinear Model and Estimation Method Nonlinear models of the total arc and its sub-systems were developed as likewise outlined in the previous chapter.  As described therein, each system was assumed to be represented by a second-order Volterra series as Eq. (2) that we repeated here:       Again,  is discrete-time;  and  are the measured input and output (i.e., CSP and AP for the total arc, CSP and SNA for the neural arc, and SNA and AP for the peripheral arc) with precisely denoting the input after removing its mean value; and , ,  are the system kernels, with memory M, for estimation.  The zeroth-order kernel  is simply the mean value of .  As before, the first-order or linear kernel , which is the time-domain version of the conventional transfer function, indicates how the present and past input samples (e.g.,  and ) affect the current output sample .  The second-order nonlinear kernel  indicates how the interaction between, or product of, two input samples that are  and  samples in the past (e.g.,  or ) impact .  While this model neglects higher order nonlinearity, many physiologic systems can be well represented with a second-order Volterra series [38].   The kernels of the total, neural, and peripheral arcs were estimated from the pre-processed SHR120 and SHR160 training data using a nonparametric, frequency-domain method as explained in the previous chapter ([41]).  The memory M was set to 25 sec, which is twice the length of the linear kernel of the total arc reported in our 40  previous study [15].  This value captured the memory of all systems (see Results).  The second-order kernel estimates were then visually examined to ultimately arrive at reduced nonlinear models with potentially more accurate kernel estimates (see Results).  Note that the second-order Uryson model of Figure 6 is one example of a reduced nonlinear model.  In this simpler model, the product of the present and past input samples of the same lag (e.g., ) affect  but not the product of past input samples of different lags (e.g., ).  Hence, while the second-order Volterra kernel is a function of two variables (i.e.,  as indicated in the above equation), the second-order Uryson kernel is only a function of one variable (i.e.,  as indicated in Figure 6).  Further, in a second-order Uryson model, the second-order kernel (i.e.,  differs in shape from the first-order kernel (), as implied in Figure 6.  This means that a second-order Volterra model may be reduced to a second-order Uryson model, if (a) the off-diagonal values of the second-order Volterra kernel are zero (i.e.,  = 0 for   ) and (b) the diagonal values (i.e.,  for  = ) are not simply proportional to the first-order kernel (i.e.,   a, where a is an arbitrary constant).      Model Testing The resulting nonlinear models with the first- and second-order kernel estimates and linear models with only the first-order kernel estimates were evaluated.  First, the inputs from the SHR120 and SHR160 training and testing data were applied to the models. Then, R2 values between the predicted and measured outputs were computed.  Finally, the R2 values from linear and nonlinear models were compared 41  after log transformation via paired t-ultiple comparisons [57].  Model Comparison to WKY The resulting validated models for SHR120 and SHR160 were compared to models for WKY (last chapter as well as [58]).  The models for WKY were developed and validated by applying similar methodology to ten age-matched WKY except that the Gaussian white noise stimulation was employed only at a mean of 120 mmHg, which is the normal CSP level of WKY.  In particular, first- and second-order kernel estimates for SHR120, SHR160, and WKY were characterized in terms of three parameters: area, to indicate the steady-state gain; absolute peak amplitude, to indicate the maximal gain; and dominant time constant (via a robust rectangular method [44]), to indicate the speed in reaching steady-state.  The area and absolute peak amplitude of the second-order kernel estimates were scaled by the standard deviation of the input so that they could be meaningfully related to the corresponding parameters of the first-order kernel estimates.  The kernel parameters for SHR120, SHR160, and WKY were then compared after log transformation using unpaired t-tests  Results Gaussian White Noise Data Figure 13 shows the pre-processed CSP, AP, and calibrated SNA from the SHR120 and SHR160 training data of one subject.  The group average (meanSE) mean and standard deviation of the pre-processed AP in the training data were 176.3±15.5 42  and 9.4±1.1 mmHg for SHR120 and 143.7±15.1 and 10.6±1.7 mmHg for SHR160, respectively.  The corresponding values for WKY from the previous chapter were 97.1±4.4 and 6.6±0.2 mmHg [58].  The mean of AP for SHR120 and SHR160 was significantly higher than that of WKY, which suggests baroreceptor resetting in SHR.  The standard deviation of AP for SHR120 and SHR160 tended to be higher than that for WKY, which hints at enhanced total arc dynamics in SHR.  The group average mean and standard deviation of the pre-processed SNA in the training data were 123.6±17.5 and 25.0±2.0 arbitrary units (au) for SHR120, 90.6±16.7 and 28.9±2.8 au for SHR160, and, again from the previous chapter, 80.7±11.8 and 20.7±0.6 au for WKY, respectively.  (Note that SNA cannot be compared between different subjects due to the SNA calibration step.)  These mean values with corresponding CSP levels and the standard deviation values define the operating point and range of applicability, respectively, of the models of the baroreflex arcs reported herein. 43   Figure 13:  Gaussian white noise training data from one subject.  SHR120 and SHR160 are spontaneously hypertensive rats during Gaussian white noise carotid sinus pressure (CSP) stimulation with mean of 120 and 160 mmHg, respectively.  44  Total Arc Model Figure 14 shows the group average first- and second-order kernels of Volterra models of the total arc estimated from the SHR120 and SHR160 training data.  Note that the inputs for the first-order and second-order kernels are CSP and CSP2, respectively, while the output for both kernels is AP.  Hence, in discrete-time, the units are mmHg/mmHg (unitless) for the first-order kernel and mmHg/mmHg2 (mmHg-1) for the second-order kernel.  The kernels are qualitatively similar to those of WKY in two ways (last chapter as well as [58]).  First, the second-order kernels revealed small off-diagonal values.  In fact, the off-diagonal energies (i.e., sums of squares of , squares of ), and all of the off-diagonal energies were statistically smaller than the diagonal energy via t-tests (p < 0.0003).  Hence, the second-order kernels were approximately diagonal.  Second, the diagonals of the second-order kernels were different in shape from the first-order kernels (see Figure 15 for clear view).  As described above, these two attributes of the second-order kernel mean that the Volterra model may be reduced to the Uryson model of Figure 6 for both SHR120 and SHR160.  The kernels of the reduced models were re-estimated using a non-parametric frequency-domain method [58] with M again set to 25 sec.  (This procedure yielded somewhat different and likely more accurate second-order Uryson kernel estimates than the diagonal of the second-order Volterra kernel estimates.)  Figure 15 shows the resulting group average Uryson kernel estimates, while the Appendix provides the numerical values of the kernel estimates including for WKY.  The four kernels were 45  similar in shape to each other as well as to the corresponding kernels for WKY (see [58]).  These kernels indicated low-pass characteristics and similar dynamics with each other (i.e., an impulse increase in CSP at time zero would cause AP to decrease and then return to baseline without oscillation).  Figure 14:  Group average first-order (linear) and second-order kernel estimates (meanSE) of complete Volterra models (see Equation) of the total baroreflex arc in SHR120 and SHR160.  The inputs for the first-order and second-order kernels are CSP and CSP2, respectively, while the output for both kernels is AP.  Hence, in discrete-time, the units are mmHg/mmHg (unitless) for the first-order kernel and mmHg/mmHg2 (mmHg-1) for the second-order kernel.  Front view precisely means the view point with azimuth of 135° and elevation of 0° with respect to the axis origin.   46   Figure 15:  Group average first- and second-order kernel estimates of reduced Uryson models of the total arc in SHR120 and SHR160.  Table 4: Group average R2 values between arterial pressure (AP) predicted by models of the total baroreflex arc and measured AP.  Linear Second-order Volterra Second-order Uryson SHR120 Training Data 0.53±0.05 0.81±0.02* 0.69±0.04* Testing Data 0.45±0.04 0.52±0.07 0.64±0.04* SHR160 Training Data 0.63±0.05 0.82±0.03* 0.71±0.05* Testing Data 0.59±0.06 0.63±0.06 0.71±0.05* * denotes statistical significance for paired t-test comparison wcorrection for three comparisons.   Table 4 shows the group average R2 values between the AP predicted by the individual subject models when stimulated by the Gaussian white noise CSP in the SHR120 and SHR160 training and testing data and the measured AP.  The training data results actually indicate model fitting ability, whereas the testing data results truly indicate model prediction capacity.  In the testing data, the linear models achieved an R2 value of only 0.450.04 for SHR120 but 0.590.06 for SHR160.  The Uryson models significantly improved upon these values by 43% for SHR120 and 21% for SHR160.  The linear and Uryson model predictive capacities for WKY were more similar to 47  those for SHR160 than SHR120 (see [58]).  Note that the predictive capacity of the Volterra models was worse than that of the Uryson models in the testing data due to overfitting in the training data.         Table 5: Group average parameters of the kernels of the validated Uryson model of the total arc.   WKY SHR120 SHR160 Linear    Kernel Area/Gain (unitless) -0.70±0.11 -0.61±0.15 -0.76±0.09 Absolute Peak Amplitude (unitless) 0.085±0.011 0.127±0.022 0.165±0.029 Time Constant (s) 4.1±0.5 2.4±0.4 2.4±0.2 Second-order Uryson Kernel Area/Gain (unitless) -0.22±0.03 -0.37±0.05 -0.38±0.04 Absolute Peak Amplitude (unitless) 0.020±0.004 0.053±0.007 0.037±0.006 Time Constant (s) 6.2±0.8 3.5±0.4 5.5±0.6 p-values Linear Kernel  WKY vs SHR120 WKY vs SHR160 SHR120 vs SHR160 Area/Gain (unitless) 0.49 0.41 0.23 Absolute Peak Amplitude (unitless) 0.12 0.021 0.31 Time Constant (s) 0.10 0.007* 0.57 Second-order Uryson Kernel Area/Gain (unitless) 0.019* 0.006* 0.65 Absolute Peak Amplitude (unitless) <0.001* 0.013* 0.091 Time Constant (s) 0.006* 0.75 0.014* WKY, Wistar Kyoto rats during Gaussian white noise carotid sinus pressure (CSP) stimulation with mean of 120 mmHg; and SHR120 and SHR160, spontaneously hypertensive rats during the same CSP stimulation with mean of 120 and 160 mmHg, respectively.  The WKY values are from last chapter, and the p-values were obtained via unpaired t-  Table 5 shows group average parameters of the first- and second-order kernels of the validated Uryson models of the total arc for SHR120 and SHR160.  This table also includes the corresponding values for WKY from last chapter.  The absolute peak amplitude and time constant of the linear kernel for SHR160 were nearly twice as large and almost half as small as those for WKY, respectively.  These differences were either significant or close to significant.  However, the area (or gain) of the linear kernel for SHR160 and all three parameters of this kernel for SHR120 were not 48  significantly different from those for WKY.  The gains and absolute peak amplitudes of the second-order kernels indicated a 20 to 60% magnitude of effect relative to their linear kernel counterparts, whereas the time constants of the second-order kernels generally indicated a slower effect than the linear kernels.  The parameters of the second-order kernels for SHR120, SHR160, and WKY were more significantly different than those of the linear kernels.  In particular, the gains and absolute peak amplitudes of the second-order kernels for SHR120 and SHR160 ranged from about 170 to 270% larger than those for WKY.  In addition, the time constant of the second-order kernel for SHR120 was about 60% smaller than those for WKY and SHR160.  All of these differences were significant.  In sum, the second-order kernel of the total arc was augmented significantly more in SHR relative to WKY than the linear kernel.  Neural Arc and Peripheral Arc Models Like our finding in WKY (last chapter as well as [58]), the first- and second-order kernels of Volterra models of the neural arc estimated from the SHR120 and SHR160 training data likewise suggested reduced Uryson models (results not shown).  Figure 16 shows the group average kernels of the Uryson models of the neural arc estimated from these data.  The four kernels appeared similar in shape to each other as well as to the corresponding kernels for WKY (see last chapter as well as [58]).  These kernels indicated high-pass characteristics and similar dynamics with each other (i.e., an impulse increase in CSP at time zero would cause SNA to decrease and then return to baseline with oscillation) and were much faster than their total arc counterparts (see Figure 15).  49  Reliable second-order kernels of Volterra models of the peripheral arc could not be estimated, because the input of this system (SNA) was not Gaussian white noise.  Figure 17 shows the group average kernels of linear models of the peripheral arc estimated from the SHR120 and SHR160 training data.  The two kernels were similar overall to each other and to the corresponding kernel for WKY (see last chapter as well as [58]).  These kernels indicated low-pass characteristics and expected open-loop dynamics (i.e., an impulse increase in SNA at time zero would cause AP to increase and then return to baseline without oscillation) and were similar in speed to their total arc counterparts.  Figure 16:  Group average first- and second-order kernel estimates of Uryson models of the neural arc in SHR120 and SHR160.  Figure 17:  Group average first-order kernel estimates of linear models of the peripheral arc in SHR120 and SHR160. 50  Table 6 shows the group average R2 values between the SNA predicted by the individual subject models of the neural arc when stimulated by the Gaussian white noise CSP in the SHR120 and SHR160 training and testing data and the measured SNA.  In the testing data, the linear models achieved high R2 values of 0.650.03 for SHR120 and 0.820.01 for SHR160.  The Uryson and Volterra models significantly improved upon the value for SHR120 by 8 to 9% but did not appreciably improve the value for SHR160 (1%).  So, the neural arc for SHR160 was approximately linear.  Table 6 also shows the group average R2 values between the AP predicted by the individual subject linear model when stimulated by the SNA in the SHR120 and SHR160 training and testing data and the measured AP.  For both the training and testing data, the R2 values were high.  In particular, in the testing data, the R2 values were 0.830.03 for SHR120 and 0.780.10 for SHR160.  Hence, the peripheral arc was approximately linear anyhow.  The predictive capacities of the models of the neural arc for SHR160 and the models of the peripheral arc for both SHR120 and SHR160 were similar to those for WKY (see last chapter as well as [58]).       Table 6: Group average R2 values between efferent sympathetic nerve activity (SNA)/AP predicted by models of the neural and peripheral arcs and measured SNA/AP.  Neural Arc Peripheral Arc  Linear Second-order Volterra Second-order Uryson Linear SHR120 Training Data 0.65±0.03 0.76±0.02* 0.73±0.02* 0.82±0.04 Testing Data 0.65±0.03 0.71±0.03* 0.70±0.03* 0.83±0.03 SHR160 Training Data 0.82±0.01 0.85±0.01* 0.84±0.01* 0.83±0.07 Testing Data 0.82±0.01 0.83±0.01 0.83±0.01* 0.78±0.10 * denotes statistical significance for paired t-correction for three comparisons.   51   As described in our previous study, extensive report on the linear dynamics of the baroreflex arcs in SHR [15], the parameters of the kernels of the validated linear models of the neural and peripheral arcs for SHR120 and SHR160 were mostly not significantly different from those for WKY to the extent that they could be compared (i.e., the gains and absolute peak amplitudes of these kernels cannot be compared between different subjects due to the SNA calibration step).  Finally, the linear and second-order kernels of the validated Uryson model of the neural arc for SHR120 were characterized by gains of -0.550.09 and -0.270.06 au/mmHg, absolute peak amplitudes of 0.930.02 and 0.270.04 au/mmHg, and time constants of 0.300.05 and 0.480.06 s, respectively. Discussion The major finding of this study is that nonlinear dynamic functioning of the sympathetically-mediated carotid sinus baroreflex is enhanced significantly more in SHR than its linear dynamic functioning.  We arrived at this result by developing and validating dynamic models of the total baroreflex arc and its sub-systems in SHR via Gaussian white noise CSP stimulation and nonlinear system identification and then comparing these models to our previously established models for age-matched WKY [58]. Nonlinear Identification Method This study falls within the large body of literature on system identification analysis of cardiovascular variability interactions (see, e.g., [59][61]).  Amongst the 52  various identification methods that have been employed, we chose a frequency-domain method [41] to estimate the kernels of the nonlinear models.  This nonparametric method assumes that the input is Gaussian and broadband but makes no assumptions on the form of the kernels.  The method yields the best unbiased estimates of the linear and nonlinear kernels in the least squares sense.  Hence, the linear term of the nonlinear model and the optimal linear model are one in the same.  Parametric identification methods, which have been widely employed in this area [59], [61], [62], could provide better estimates of the kernels by assuming a particular kernel form so as to trade off bias for precision.  This possibility could especially hold in the identification of the peripheral arc whose input was not as broadband as the other investigated systems.  However, our goal was to determine the form of the kernel.  Even so, we did apply conventional autoregressive exogenous input identification to estimate the linear kernels, and this method did not yield more predictive linear kernels than the frequency-domain method (results not shown).        Total Arc Model in SHR We applied the frequency-domain method to develop second-order nonlinear dynamic models of the total arc for both SHR120 (SHR with Gaussian white noise stimulation at the normal CSP level for WKY) and SHR160 (SHR with the same simulation but at the prevailing CSP level for SHR) (see Figure 13).  These models were qualitatively similar to the corresponding model for WKY.  In particular, they generally indicated that the total arc may be represented as a linear dynamic system in parallel with a cascade combination of a squarer and a slower, linear dynamic system 53  (see Figure 6, Figure 14, and Figure 15).  Hence, total arc nonlinearity, which was captured by the squaring and slower, linear dynamic systems, was in the low frequency regime.    linear models by 43% for SHR120 and 21% for SHR160 (see Table 4).  The predictive capacity of the Uryson model relative to a linear model for SHR120 was superior to those for WKY and SHR160.  The reason for the relatively stronger nonlinearity in SHR120 may pertain to the operating point.  That is, the CSP levels for WKY and SHR160 are in the linear regimes of their respective static sigmoidal CSP-AP relationships, while the CSP level for SHR120 is near the nonlinear thresholding regime of its relationship [63].  In this sense, comparisons between WKY and SHR160, as opposed to SHR120, may actually be more appropriate.  The validated models of the total arc for SHR120 and SHR160 were, however, quantitatively different from the corresponding model for WKY (see Table 5).  More specifically, the linear kernel (time-domain version of the transfer function of the faster, dynamic system) for SHR160 showed significantly augmented transient dynamics in terms of magnitude and speed than the corresponding kernel for WKY.  The linear kernel for SHR120 showed some tendency for similarly enhanced transient dynamics.  However, the linear kernels for SHR120, SHR160, and WKY revealed steady-state gains, which are more meaningful, of very similar values.  By contrast, the nonlinear kernels (time-domain version of the transfer function of the slower, dynamic system) for SHR120 and SHR160 showed enhanced steady-state gains by about 170% and significantly augmented transient dynamics in terms of magnitude 54  (and speed for SHR120) relative to the corresponding kernel for WKY.  Note that while the linear kernels change AP in the opposite direction of CSP, the nonlinear kernels always reduce AP due to the squaring of CSP (see Figure 6).  Hence, in normal, closed-loop conditions, the nonlinear steady-state gain of the total arc may augment buffering of AP increases, while blunting buffering of AP decreases, to a greater extent in SHR than WKY.  Further, the nonlinear steady-state gain of the total arc may decrease mean AP in response to increases in the AP variance to greater extent in SHR than WKY. Neural Arc and Peripheral Arc Models in SHR We also developed second-order Uryson models of the neural arc for SHR120 and SHR160 (see Figure 16).  However, linear models of the neural arc showed good to excellent SNA prediction, so the Uryson model for SHR160 did not improve upon the SNA prediction (see Table 6).  Likewise, a linear model of the neural arc sufficed for WKY [58].  While the neural arc model for SHR120 was able to improve SNA prediction perhaps due to relatively stronger nonlinearity at the different operating point (see Table 6), the 8% improvement achieved was small compared to the 43% AP prediction improvement attained by the Uryson model of the total arc for SHR120.  Similar to WKY, we could only develop linear models of the peripheral arc for SHR120 and SHR160 (see Figure 17), as the SNA input to this system was not Gaussian white noise.  However, these standard models showed excellent AP prediction anyhow.  As we described in a previous, extensive report [15], the kernels of the validated linear models of the neural and peripheral arcs for SHR120 and SHR160 were 55  mostly not significantly different from those for WKY.  However, comparisons of these kernels were substantially limited due to the SNA calibration step.           In sum, the total arc models showed nonlinear behaviors, while the models of the neural and peripheral arc sub-systems, especially for WKY and SHR160, showed approximately linear behaviors.  As we discussed earlier [58], one explanation for this puzzling result is that the neural arc was nonlinear, but its nonlinearity was not well identified due to a linearizing effect caused by confounding SNA from higher brain centers.  Note that the improved predictive capacity of the Uryson model of the neural arc for SHR120, wherein nonlinearity may have been relatively stronger, supports the contention that the neural arc was nonlinear. Potential Physiologic Mechanisms This study does not reveal the mechanisms underlying the more significant enhancement of nonlinear dynamic functioning of the total arc in SHR.  However, previous studies shed a bit of insight.  Firstly, the carotid artery stiffens in SHR [64].  Since the carotid sinus baroreflex precisely responds to stretch, such stiffening alone would suggest blunted total arc functioning in SHR.  Hence, enhanced functioning downstream in the total arc must have occurred in terms of linear dynamics and, to a greater extent, nonlinear dynamics.  Secondly, baroreflex control of heart rate is blunted in SHR after vagal block [15], [65].  So, downstream dynamic functioning pertaining to the control of vascular properties and/or cardiac contractility must have specifically been enhanced.  Finally, we note that it may be possible that the linear component (faster dynamic system) and nonlinear component (squaring and slower 56  dynamic systems) of the Uryson model correspond to myelinated and unmyelinated fibers pathways of the total arc [51], [52] and that these pathways are differentially impacted in SHR relative to WKY. Study Limitations As we outlined previously, this study has experimental and mathematical limitations.  More specifically, our experimental procedures included the use of anesthesia, opening the baroreflex loop, and elimination of the vagal component of the baroreflex, while our mathematical procedures neglected higher-order nonlinearity.  These procedures are limitations for sure.  At the same time, they may have permitted as accurate identification of the sympathetically-mediated total arc as possible without assuming any form for the nonlinearity.   57  CHAPTER 5. NONLINEAR IDENTIFICATIOND OF THE TOTAL BAROREFLEX ARC: HIGHER-ORDER NONLINEARITY ---- Nonlinear Model and Estimation Method --           4 -- 58  -----           --, , , and -The first-order or linear kernel , which is the time-domain version of the conventional transfer function, indicates how the present and past input samples (e.g.,  and ) affect the current output sample .  The second-order kernel , which is Uryson in structure (i.e., diagonal) based on the previous work, indicates how the present and past squared input samples (e.g.,  and ) .  Note that because of the Uryson structure, products of pairs of present and past input samples of different lags (e.g., ) do not impact .  Finally, the third-order kernel 59   is scaled by and  is scaled by .--               (6) -, , , and -       60     The kernels , , , and  in Eq. (6) were estimated from  after removing its mean value and  using the frequency-domain method for Gaussian inputs.  This method, which minimizes the mean-squared output prediction error and yielded relatively smooth kernel estimates, is described in detail elsewhere [68].  Briefly, the kernels were estimated in succession according to the following four steps:   1.  where  is the expectation operator.  2.   where  is the auto- or cross-correlation function between the indicated signals, and  is the standard Fourier Transform operator. 3.    4.   where  is the three-dimensional Fourier Transform operator.  Note that  and  above were computed via the standard sample mean and unbiased correlation 61  function estimates.  We also tried the cross-correlation method for Gaussian white noise  to estimate kernels.  However since the results were not smooth, we just described the frequency domain method.   -----    The kernels , , , and  were again estimated from  after removing its mean value and  using the frequency-domain method for Gaussian inputs.  This method yielded relatively smooth kernel estimates.  The kernels were estimated in succession according to the first three steps in the last paragraph and then the following fourth step: 4.     62  The four kernels were also estimated from the same data using the Laguerre expansion method for Gaussian inputs [50].  This more parsimonious method produced similar kernel estimates.  More specifically, the first-, second-, and third-order kernels were assumed to be represented by a linear combination of  Laguerre basis functions as follows:   Here,  is the jth coefficient of the ith-order kernel for estimation, and  is the jth Laguerre basis function.  This function is given as follows:  where  () is a smoothing parameter that determines the rate of exponential decay of the Laguerre function.  The coefficients of each kernel were estimated in succession by minimizing the mean-squared output prediction error.  That is, for fixed  and  values, the  coefficients for each kernel were estimated via solution of the following  linear system of normal equations:    for  where  is the auto- or cross-correlation function between the indicated signals, and  is the output after subtracting the contribution of all previously estimated kernel.  To determine and , the same  value was assumed for all kernels, and the  value was assumed to be between 1 and 10.  For a fixed  value, the value of  was determined for each kernel by finding the value for which the mean-squared output prediction error no longer significantly decreased.  This calculation was repeated for 63  each  value between 0 and 1, and the value that minimized the mean-squared output prediction error was selected.  See the Discussion section for commentary on the various methods.  In all cases, the kernel  Model Evaluation ---  Results ----Figure 18-Figure 19-64  ------              Figure 18: Group average two-dimensional (2D) slices of the estimated (four-dimensional) third-order kernels () in Eq. (5).  WKY, Wistar Kyoto rats during Gaussian white noise carotid sinus pressure (CSP) stimulation with mean of 120 mmHg; and SHR120 and SHR160, spontaneously hypertensive rats during the same CSP stimulation with mean of 120 and 160 mmHg, respectively.  The diagonal slice () was largest in magnitude.  65   Figure 19: Group average energy of each 2D slice of the estimated second-order Volterra kernels from [58], [66] and third-order kernels, normalized by the diagonal energy, in descending order of value for WKY, SHR120, and SHR160. The second-order kernel was always diagonal, while the third-order kernel was virtually diagonal for SHR120 and approximately diagonal for WKY and SHR160.     Figure 20: Group average first-, second-, and third-order kernel estimates (meanSE) of reduced Uryson models of the total arc for WKY, SHR120 and SHR160 via the frequency-domain (black) and Laguerre expansion (gray) methods for Gaussian inputs.  The three kernels were not proportional to each other, thereby indicating that the model could not be further reduced.  --Figure 2066  ---------- Table 7--Table 2Table 4--------67  - Table 7:  Group average R2 values between arterial pressure (AP) predicted by models of the total baroreflex arc and measured AP. Values represent meanSE.  Third-order Volterra model refers to Eq. (5).  WKY, Wistar Kyoto rats during Gaussian white noise carotid sinus pressure (CSP) stimulation with mean of 120 mmHg; and SHR120 and SHR160, spontaneously hypertensive rats during the same CSP stimulation with mean of 120 and 160 mmHg, -three comparisons) with corresponding linear model and second-order Uryson model, respectively.  Discussion ---- Linear Second-order Uryson Third-order Volterra  Third-order Uryson (Frequency-domain) Third-order Uryson  (Laguerre expansion) WKY Training Data 0.73±0.03 0.79±0.03* 0.90±0.01* 0.80±0.03* 0.80±0.03* Testing Data 0.64±0.03 0.71±0.03* 0.69±0.03* 0.71±0.03* 0.71±0.03* SHR120 Training Data 0.53±0.05 0.69±0.04* 0.88±0.01* 0.72±0.05* 0.71±0.04* Testing Data 0.45±0.03 0.64±0.03* 0.64±0.04* 0.65±0.04* 0.66±0.04* SHR160 Training Data 0.63±0.05 0.71±0.05* 0.88±0.03* 0.72±0.05* 0.72±0.05* Testing Data 0.59±0.04 0.71±0.03* 0.70±0.03* 0.73±0.03* 0.72±0.03* 68  --    We used a third-order Volterra model but with a second-order Uryson kernel (see Eq. (5)).  We estimated the first- and second-order kernels using the frequency-domain method for Gaussian inputs, which we previously found to yield the smoothest kernels based on visual assessment [58], [66].  Since estimation of higher-order kernels from short data periods (6 min) is challenging, we applied several identification methods to estimate the third-order kernels.  The nonparametric identification methods, which do not assume a particular kernel form, included the frequency-domain method for Gaussian inputs and cross-correlation method [50], which assumes the input is Gaussian and white.  We again found that the frequency-domain method for Gaussian inputs was most effective, as the CSP input was not strictly white due to the finite data periods.  The resulting third-order kernel estimates likewise revealed a simpler Uryson structure for WKY, SHR120 (SHR with Gaussian white noise stimulation at the normal CSP level for WKY), and SHR160 (SHR with the same stimulation at the prevailing CSP level for SHR) but not as strongly for WKY and SHR160 (see Figure 18 and Figure 19).  It is possible that the diagonal nature of the third-order kernels for WKY and SHR160 was partially masked by noise arising from the identification of higher-order kernels from short data periods in the presence of stronger linearity (see Table 7).  Further, since the second-order kernel, which is easier to estimate, is surely Uryson (see Figure 19), the third-order kernel may indeed likewise show such structure.  At the same time, we acknowledge the possibility of nontrivial error in the third-order kernel estimates, which could have 69  confounded the interpretation of the results.  Parametric identification methods could possibly improve kernel estimation from short data periods by reducing the number of estimated parameters via assumption of a particular kernel form.  We also applied such methods including popular Laguerre expansion methods [50].  However, these methods artificially smoothed the kernels so as to obfuscate a diagonal appearance.  We therefore moved to a reduced, third-order Uryson model (see Eq. (7)).  We again applied various identification methods to re-estimate the kernels of this model with potentially greater accuracy.  The methods included the frequency-domain method for Gaussian inputs, cross-correlation method, and Laguerre expansion methods for Gaussian and arbitrary inputs.  We found that the frequency-domain method and Laguerre expansion method for Gaussian inputs yielded the smoothest estimates based on visual inspection, as the input was not strictly white and the requisite computation of higher-order correlations was more robust by leveraging the Gaussian nature of the inputs.  The resulting kernel estimates were not proportional to each other (see Figure 20), thereby indicating that the third-order model could not be further reduced to a Hammerstein model [43].  The Laguerre expansion method for Gaussian inputs represented the third-order Uryson model with only 19 parameters on average (6 parameters for each of the three kernels plus 1 smoothing parameter).  For comparison, the frequency-domain method for Gaussian inputs represented the second-order Uryson model with 100 parameters (M = 25 sec times 2 Hz sampling rate for each of the two kernels; see Eq. (7)).  Hence, the third-order Uryson model was likely well estimated.  Neither the third-order Uryson model nor the more complete third-order 70  nonlinear model of Eq. (5) predicted AP in response to new Gaussian white noise CSP better than the previously established second-order Uryson model (see Table 7).  Further, both of the third-order models could not offer added value in predicting well-known nonlinear behaviors over the second-order Uryson model (results not shown).  Since there is a possibility that the third-order kernel was actually not Uryson, we formed a set of third-order kernels using the highest energy slices of the third-order Volterra kernel estimate (see Figure 19).  These kernels are less simplified than the Uryson kernel but more simplified than the Volterra kernel.  However, the third-order models with these kernel estimates also did not improve the AP prediction in any way (results not shown).  Hence, even though the third-order kernels were not merely noise (see Figure 20), they may be interpreted as small in magnitude in the sense of contributing relatively little to AP prediction.  Note that a nonlinear model including only the first- and third-order kernels would not have performed better due to orthogonality arising from use of Gaussian inputs.  Also note that while fourth-order and even higher-order kernels could possibly contribute to AP prediction, since the first-order kernel, which is odd, and the second-order kernel, which is even, both contributed significantly, we doubt that kernels of order exceeding three would be important here. --71  --- shed We conclude that a second-order Uryson model of the total arc surely suffices compared to a third-order Uryson model and may likely suffice compared to higher-order nonlinear models in general.  This conclusion is only valid over the range of CSP and AP data utilized to develop the models.  Since this range did not include the baroreflex saturation regime, the nonlinear models could only predict baroreflex thresholding but not saturation [58].  Higher order nonlinear models would be needed to represent the total arc over the entire system operating range, with odd-order nonlinearity apparently required to account for thresholding and saturation.  In sum, this third part of a series of studies on nonlinear identification of the total arc justified 72  the earlier assumption of second-order nonlinearity under the employed experimental conditions.    Figure 21: Representative time series for arterial pressure (AP) and sympathetic nerve activity (SNA; measured from splanchnic nerve and then normalized as described elsewhere [58]) in response to fixed and Gaussian white noise CSP stimulations for one WKY (left panel).  A.u. is arbitrary units.  Group average power spectra (mean±SE) for CSP, AP, and SNA (right panel).  The gray and black lines indicate the fixed and Gaussian white noise CSP inputs, respectively.    73  CHAPTER 6. CONCLUSIONS The baroreflex was long believed to regulate AP only on the time scales of seconds to minutes [54].  However, chronic baroreceptor stimulation and other studies have now indicated that this system could contribute to long-term AP regulation [3], [4], [25], [55], [14].  Hence, the baroreflex could play a causative or protective role in hypertension and heart failure.  As indicated by Thrasunloading model of hypertension [3], nonlinearity of the baroreflex, in particular, could induce sustained increases in AP.  However, baroreflex nonlinearity is not well understood.   In the first part of this thesis, we developed a second-order, nonlinear dynamic model of the sympathetically-mediated total baroreflex arc by employing Gaussian white noise stimulation and nonparametric identification .  We validated the model by showing that it is able to predict AP appreciably better than a conventional linear model and some important nonlinear behaviors including thresholding and DC responses to AC changes.  A key advantage of nonlinear identification over linear identification is that it can indicate the structure of the system under study.  For example, consider a system that is composed of a static nonlinearity in cascade with a linear dynamic system.  Nonlinear identification of the overall system can not only yield a more accurate model than linear identification but also reveal whether the static nonlinearity precedes or follows the linear dynamic system.  The validated nonlinear model likewise provides information about the structure of the internal components of the total arc.  In particular, the model 74  illustrates that the structure is not a previously proposed cascade connection of systems (e.g., Hammerstein, Wiener, or sandwich model) but rather a parallel connection of a linear system and Hammerstein system.  In addition to providing information about the structure of the total arc, the model also indicates that the second-order nonlinear dynamics are slower than the linear dynamics and not insignificant in magnitude of effect compared to these dynamics.  While system identification cannot reveal the physical basis of the nonlinearity, it is interesting to speculate that the parallel connection of a fast linear system and a slower nonlinear system corresponds to the structure of anatomical components.  In particular, the linear system may correspond to the fast, afferent A-fiber pathway, whereas the nonlinear system could correspond to slower, afferent C-fiber pathway.  However, the model is only valid over the CSP range covered by the Gaussian white noise.  This range essentially did not include the saturation regime or pulsatile frequencies.   In the second part, we showed the pitfall in ignoring baroreflex nonlinearity by developing the second-order, nonlinear dynamic model of the sympathetically-mediated total baroreflex arc utilizing Gaussian white noise stimulation and nonparametric identification in spontaneously hypertensive rats.  Results herein indicate that the nonlinear gain of the sympathetically-mediated total baroreflex arc was enhanced in SHR relative to WKY, while its linear gain was preserved.  Hence, the nonlinear dynamic functioning of this system may enhance steady-state baroreflex buffering of AP increases more in SHR than WKY, perhaps to compensate for malfunctioning of other regulatory systems.  If the common linearity assumption were invoked here, the story would have been different.  This study is not the first to 75  chronic baroreceptor unloading model of hypertension [3], [4], mean CSP did not change but carotid sinus pulse pressure decreased, which, in turn, led to a sustained, baroreflex-mediated increase in mean AP.  This nonlinear behavior of the carotid sinus baroreflex, by contrast, played a causative rather than protective role in the hypertension model.   In the last part, we buttressed the major previous assumption that the baroreflex nonlinearity did not extend beyond second-order under the employed experimental conditions.  We --- Future studies are needed to:  (1) elucidate baroreflex nonlinearity over a wider system operating range than that attained by the experimental conditions herein (i.e. cover saturation, thresholding, and pulsatile frequencies); (2) investigate baroreflex nonlinearity in closed-loop conditions (via parametric identification [29] and the established second-order Uryson model); and (3) discover the mechanism of baroreflex nonlinearity (which may be related to the slow, afferent C-fiber pathway [51], [52]) including the unexplained AP variations.  Such future nonlinear modeling endeavors may enhance our understanding of the baroreflex in health and disease.  Specifically, further investigations of baroreflex nonlinearity in hypertension may improve our understanding of the role of the baroreflex in this prevalent disease process.   76  APPENDIX The numerical values specifying the group average kernel estimates of the Uryson model (i.e.  and ) of the total arc for SHR120 and SHR160 as well as for WKY in [58] are provided in Table 8.  Please note that only 17 sec of values were needed to capture these kernels.                  77  Table 8: Kernel values of Uryson model of the total baroreflex arc in WKY, SHR120, and SHR160 Time (s)    First-order (unitless) Second- order (mmHg-1) First-order (unitless) Second- order (mmHg-1) First-order (unitless) Second- order (mmHg-1) 0 0.0053441 0.0000888 0.0106978 0.0004215 0.0141826 0.0003052 0.5 -0.0113258 -0.0000093 -0.0317528 -0.0005482 -0.0305524 0.0001441 1 -0.0606674 -0.0004854 -0.1130225 -0.0024845 -0.1424666 -0.0005875 1.5 -0.0846313 -0.0009561 -0.1186430 -0.0032550 -0.1650641 -0.0015852 2 -0.0690056 -0.0010924 -0.0644120 -0.0027865 -0.1032691 -0.0021265 2.5 -0.0594982 -0.0010686 -0.0475934 -0.0024280 -0.0786241 -0.0021722 3 -0.0595114 -0.0010640 -0.0530959 -0.0022385 -0.0783448 -0.0021070 3.5 -0.0531564 -0.0010605 -0.0433527 -0.0019858 -0.0578963 -0.0020268 4 -0.0487779 -0.0010096 -0.0361096 -0.0016735 -0.0485531 -0.0018949 4.5 -0.0458009 -0.0009050 -0.0308038 -0.0012624 -0.0435933 -0.0017367 5 -0.0385619 -0.0007954 -0.0182875 -0.0010134 -0.0269673 -0.0015018 5.5 -0.0335319 -0.0007345 -0.0118544 -0.0009023 -0.0166916 -0.0013480 6 -0.0287025 -0.0006717 -0.0086836 -0.0007416 -0.0123483 -0.0012249 6.5 -0.0227861 -0.0005876 -0.0047388 -0.0006102 -0.0073592 -0.0009416 7 -0.0198903 -0.0005125 -0.0026564 -0.0004946 -0.0053071 -0.0007900 7.5 -0.0164859 -0.0004362 0.0015160 -0.0003245 -0.0015483 -0.0007129 8 -0.0117555 -0.0003785 0.0052325 -0.0002953 0.0016648 -0.0004655 8.5 -0.0089331 -0.0003444 0.0029448 -0.0002883 0.0016336 -0.0003715 9 -0.0063109 -0.0002790 0.0018112 -0.0001505 0.0026040 -0.0004274 9.5 -0.0033516 -0.0002317 0.0018556 -0.0000881 0.0023051 -0.0003040 10 -0.0010963 -0.0002079 -0.0005557 -0.0000920 0.0017833 -0.0001827 10.5 0.0016547 -0.0001710 -0.0010800 -0.0000583 0.0040070 -0.0002130 11 0.0020230 -0.0001363 -0.0019769 -0.0000554 0.0051331 -0.0001643 11.5 0.0008485 -0.0001073 -0.0036355 0.0000160 0.0042096 -0.0000540 12 0.0018865 -0.0000624 -0.0031750 0.0000851 0.0046813 -0.0000425 12.5 0.0020187 -0.0000454 -0.0030819 0.0000508 0.0038249 -0.0000237 13 -0.0001120 -0.0000472 -0.0040712 0.0000144 0.0027484 -0.0000296 13.5 -0.0011690 -0.0000349 -0.0024378 0.0000264 0.0029665 -0.0000446 14 -0.0014246 -0.0000237 -0.0006728 0.0000344 0.0022949 -0.0000434 14.5 -0.0019218 -0.0000164 -0.0005096 0.0000342 0.0016824 -0.0001204 15 -0.0017010 -0.0000233 -0.0005132 0.0000203 0.0011306 -0.0001845 15.5 -0.0018849 -0.0000378 -0.0015177 0.0000214 0.0003206 -0.0001653 16 -0.0023700 -0.0000380 -0.0025900 0.0000599 -0.0000272 -0.0001353 16.5 -0.0021489 -0.0000316 -0.0028518 0.0000561 -0.0001561 -0.0000861 17 -0.0020097 -0.0000260 -0.0031153 0.0000217 -0.0000403 -0.0000537   78  BIBLIOGRAPHY79  BIBLIOGRAPHY [1] J. E. Hall, Guyton and Hall Textbook of Medical Physiology, 12th ed. Saunders, 2010. [2] R. Klabunde, Cardiovascular Physiology Concepts, vol. 3. Lippincott Williams & Wilkins, 2011. [3] Am J Physiol Regul Integr Comp Physiol, vol. 282, no. 4, pp. R104453, Apr. 2002. [4]  in Am J Physiol Regul Integr Comp Physiol, vol. 288, no. 4, pp. R86371, Apr. 2005. [5] B. G. Katzung, Basic & Clinical Pharmacology, 11th ed. Lange Medical Books/McGraw Hill, 2009. [6] Handbook of Physiology. The Cardiovascular System. Peripheral Circulation and Organ Blood Flow, A. P. Soc., Ed. Bethesda, MD, 1983, pp. 497555. [7] T. H. Desai, J. C. Collins, M. Snell, and R. Mosqueda-arterial and carAm J Physiol, vol. 272, no. 5 Pt 2, pp. H234352, May 1997. [8] control by arterial and cardiopulmonary baroreceptors in awake dogs with Am J Physiol, vol. 257, no. 6 Pt 2, pp. H204858, Dec. 1989. [9] R. M. Oren, H. P. Schobel, R. M. Weiss, W. Stanford, and D. W. Ferguson, normal humJ Appl Physiol, vol. 74, no. 6, pp. 267280, Jun. 1993. [10] J Physiol, vol. 227, no. 1, pp. 24360, Dec. 1972. [11] cardiopulmonary baroreceptor control of splanchnic and forearm vascular J Physiol, vol. 286, no. 1, pp. 173184, Jan. 1979. [12] -pulmonary vein baroreflex on Am J Physiol, vol. 233, no. 5, pp. H58791, Nov. 1977. [13] J Physiol, vol. 50, no. 2, pp. 6584, 1915. 80  [14] Hypertension, vol. 57, no. 5, pp. 8806, May 2011. [15] T. Kawada, S. Shimizu, A. Kamiya, Y. Sata, K. Uemura, and M. Sugimachi, Am J Physiol Regul Integr Comp Physiol, vol. 300, no. 1, pp. R15565, Jan. 2011. [16] T. Sato, T. Kawada, M. Inagaki, T. Shishido, M. Sugimachi, and K. Sunagawa, Am J Physiol Regul Integr Comp Physiol, vol. 285, no. 1, pp. R26270, Jul. 2003. [17] T. Kawada, T. Sato, M. Inagaki, T. Shishido, T. Tatewaki, Y. Yanagiya, C. Zheng, -loop identification of carotid sinus Jpn J Physiol, vol. 50, no. 3, pp. 37180, Jun. 2000. [18] T. Kubo, T. Imaizumi, Y. Harasawa, S. Ando, T. Tagawa, T. Endo, M. Shiramoto, Am J Physiol, vol. 270, no. 3 Pt 2, pp. H105462, Mar. 1996. [19] Y. Ikeda, T. Kawada, M. Sugimachi, O. Kawaguchi, T. Shishido, T. Sato, H. Am J Physiol, vol. 271, no. 3 Pt 2, pp. H88290, Sep. 1996. [20] T. Kawada, M. Li, A. Kamiya, S. Shimizu, K. Uemura, H. Yamamoto, and M. -loop dynamic and static characteristics of the carotid sinus baroreflex in rats with chronic heart failuJ Physiol Sci, vol. 60, no. 4, pp. 28398, Jul. 2010. [21] H. Miyano, Y. Nakayama, T. Shishido, M. Inagaki, T. Kawada, T. Sato, H. regulation of lAm J Physiol, vol. 275, no. 2 Pt 2, pp. H4008, Aug. 1998. [22] T. Sato, T. Kawada, T. Shishido, M. Sugimachi, J. Alexander, and K. Sunagawa, aroreflex failure: a bionic baroreflex Circulation, vol. 100, no. 3, pp. 299304, Jul. 1999. [23] T. Kawada, Y. Yanagiya, K. Uemura, T. Miyamoto, C. Zheng, M. Li, M. -size dependence of the baroreflex neural arc tAm J Physiol Hear. Circ Physiol, vol. 284, no. 1, pp. H40415, Jan. 2003. [24] C. Heymans and E. Neil, Reflexogenic Areas of the Cardiovascular System, vol. 51, no. 9. Boston, MA: Royal Society of Medicine Press, 1958. [25] T. E. Lohmeier, E. D. Irwin, M. A. Rossing, D. J. Serdar, and R. S. Kieval, Hypertension, vol. 43, no. 2, pp. 30611, Feb. 2004. 81  [26] R. Mukkamala, J.-K. Kim, Y. Li, J. Sala-Mercado, R. L. Hammond, T. J. Scislo, resistance baroreflex gain values: validation by chronic arterial baroreceptor Am J Physiol Hear. Circ Physiol, vol. 290, no. 5, pp. H18306, May 2006. [27] Handbook of Physiology, American Physiology Society, 1983. [28] L. Ljung, System Identification: Theory for the User, 2nd ed. Englewood Cliffs: Prentice Hall, 1999. [29] X. Xiao, T. J-signal Physiol Meas, vol. 26, no. 3, pp. R4171, Jun. 2005. [30] le Circ Res, vol. 61, no. 5, pp. 64858, Nov. 1987. [31] Circ Res, vol. 46, no. 1, pp. 110, Jan. 1980. [32] R. M. Schmidt, M. Kumada, Am J Physiol, vol. 223, no. 1, pp. 17, Jul. 1972. [33] Circ Res, vol. 12, pp. 15262, Mar. 1963. [34] Am J Physiol, vol. 257, no. 2 Pt 2, pp. H46572, Aug. 1989. [35] W. H. Levison, G. O. Barnett, and W. D. JacksoCirc Res, vol. 18, no. 6, pp. 67382, Jul. 1966. [36] Ann N Y Acad Sci, vol. 156, no. 2, pp. 77986, May 1969. [37] H. Miyano, Y. Nakayama, T. Shishido, M. Inagaki, T. Kawada, T. Sato, H. Am J Physiol Hear. Circ Physiol, vol. 275, no. 2, pp. H400408, Aug. 1998. [38] P. Z. Marmarelis and V. Z. Marmarelis, Analysis of physiological systems: the white-noise approach. New York: Plenum Press, 1978. [39] of nonlinear biological Ann Biomed Eng, vol. 24, no. 4, pp. 25068, 1996. [40] IEEE Trans. Circuits Syst, vol. 32, no. 82  11, pp. 11501161, Nov. 1985. [41] Technometrics, vol. 3, no. 4, pp. 563567, Nov. 1961. [42] BMJ, vol. 312, no. 7038, p. 1079, Apr. 1996. [43] cascade of Hammerstein models for the description of nonlinearities in vibrating J Sound Vib, vol. 330, no. 5, pp. 10181038, Feb. 2011. [44] B. P. Lathi, Linear Systems and Signals. New York, NY: Oxford University Press, USA, 2004. [45] -linear System by Cross-Int J Control, vol. 2, no. 3, pp. 237254, Sep. 1965. [46] Ann Biomed Eng, vol. 16, no. 2, pp. 20114, Jan. 1988. [47] Ann Biomed Eng, vol. 16, no. 1, pp. 12342, Jan. 1988. [48] S. O. Haykin, Adaptive Filter Theory, 5th ed. Prentice Hall, 2013. [49] R. Mukkamala, A. T. Reisner, H. M. Hojman, R. G. Mark, and R. J. Cohen, IEEE Trans Biomed Eng, vol. 53, no. 3, pp. 45967, Mar. 2006. [50] V. Z. Marmarelis, Nonlinear Dynamic Modeling of Physiological Systems, 1st ed. Wiley, 2004. [51] J. L. Seagard, F. A. Hopp, H. A. Drummond, and D. M. Van Wynsberghe, blood pressurCirc Res, vol. 72, no. 5, pp. 101122, May 1993. [52] baroreceptor C-Acta Physiol Scand, vol. 166, no. 3, pp. 16774, Jul. 1999. [53] T. Kawada, M. Sugimachi, T. Sato, H. Miyano, T. Shishido, H. Miyashita, R. -loop identification of carotid sinus baroreflex open-loop transfer characteristics in Am J Physiol, vol. 273, no. 2 Pt 2, pp. H102431, Aug. 1997. [54] A. C. Guyton, Textbook of Medical Physiology. Elsevier Saunders, 2006. [55] -term control of Curr Hypertens Rep, vol. 8, no. 3, pp. 24954, Jun. 2006. [56] R. K. Pearson, Discrete-time Dynamic Models, vol. 8. Oxford University Press, 83  USA, 1999. [57] Scand J Stat, vol. 6, no. 2, pp. 6670, 1979. [58] M. Moslehpour, T. Kawada, K. Sunagawa, M. Sugimachi, and R. Mukkamala, Am J Physiol Regul Integr Comp Physiol, vol. 309, no. 12, pp. R1479R1489, Dec. 2015. [59] rocessing for Phil Trans R Soc A, vol. 367, no. 1887, pp. 391409, Jan. 2009. [60] -Based Signal Processing for the Analysis of Short-Term CardioProc IEEE, vol. 94, no. 4, pp. 805818, Apr. 2006. [61] J. A. Sala-Mercado, M. Moslehpour, R. L. Hammond, M. Ichinose, X. Chen, S. baroreflex enhances ventricular contractility in awake dogs: a mathematical Am J Physiol Regul Integr Comp Physiol, vol. 307, no. 4, pp. R45564, Aug. 2014. [62] A. Porta, L. Faes, A. Marchi, V. Bari, B. De Maria, S. Guzzetti, R. Colombo, and gling cardiovascular control mechanisms during head-down tilt via joint transfer entropy and self-Front Physiol, vol. 6, p. 301, Jan. 2015. [63] Y. Sata, T. Kawada, S. Shimizu, A. Kamiya, T. Akiyama, and M. Sugimachi, role of neural arc in sympathetic baroreflex resetting of Circ J, vol. 79, no. 3, pp. 5929, Jan. 2015. [64] Hypertension, vol. 3, no. 4, pp. 48595, Jan. . [65] function and end-Am J Physiol, vol. 277, no. 3 Pt 2, pp. H12006, Sep. 1999. [66] M. Moslehpour, T. Kawada, K. Sunagawa, M. Sugimachi, and R. Mukkamala1, Am J Physiol Regul Integr Comp Physiol, vol. Jan 16, no. [Epub ahead of print], 2016. [67] J. F. Barrett, -linear Physical J Electron Contr, vol. 15, no. 6, pp. 567615, Dec. 1963. [68] between fluctuations with nonlineProc IEEE, vol. 68, no. 8, pp. 10261027, 1980. [69] M. Schetzen, The Volterra and Wiener theories of nonlinear systems. Wiley, 1980.