Salt-sensitive hypertension is known to be associated with dysfunction of the baroreflex control system in the Dahl salt-sensitive (SS) rat. However, neither the physiological mechanisms nor the genomic regions underlying the baroreflex dysfunction seen in this rat model are definitively known. Here, we have adopted a mathematical modeling approach to investigate the physiological and genetic origins of baroreflex dysfunction in the Dahl SS rat. We have developed a computational model of the overall baroreflex heart rate control system based on known physiological mechanisms to analyze telemetry-based blood pressure and heart rate data from two genetic strains of rat, the SS and consomic SS.13BN, on low- and high-salt diets. With this approach, physiological parameters are estimated, unmeasured physiological variables related to the baroreflex control system are predicted, and differences in these quantities between the two strains of rat on low- and high-salt diets are detected. Specific findings include: a significant selective impairment in sympathetic gain with high-salt diet in SS rats and a protection from this impairment in SS.13BN rats, elevated sympathetic and parasympathetic offsets with high-salt diet in both strains, and an elevated sympathetic tone with high-salt diet in SS but not SS.13BN rats. In conclusion, we have associated several important physiological parameters of the baroreflex control system with chromosome 13 and have begun to identify possible physiological mechanisms underlying baroreflex impairment and hypertension in the Dahl SS rat that may be further explored in future experimental and modeling-based investigation.
- computational model
- baroreceptor reflex
- autonomic nervous system
the dahl salt-sensitive (SS) rat is an extensively studied model of human salt-sensitive hypertension. Consomic substitution studies involving the introgression of whole chromosomes from a salt-resistant rat strain into the isogenic background of the SS strain, congenic substitution studies involving the introgression of smaller regions of chromosomes, and other genetic studies have led to the identification of many quantitative trait loci (QTLs) associated with this disease (12, 36, 44, http://www.pga.mcw.edu). However, phenotypes collected from many of these studies often involve relatively gross measurements such as blood pressure and heart rate, which provide little information on the underlying physiology. Many protective QTLs have been shown to interact epistatically (i.e., in a nonlinear or nonadditive manner) so that these measurements are unable to distinguish various combinations of introgressed QTLs (36). It is also likely that many QTLs have been left unidentified because of these interactions. Thus, these types of measurements become diminishingly informative with an increased degree of genetic nonlinearity.
It seems, then, that more detailed phenotypic measurements are required to understand the underlying etiology and to make sense of the genetics associated with this complex disease. Of course, this is not always possible; many measurements of interest are either inaccessible or simply not practical to obtain. In addition, many of these measurements are operating-point dependent and are influenced to a high degree by physiologic state. Methods of obtaining these measurements often require invasive techniques that introduce stressors (surgical, pharmacological, etc.) that may themselves alter physiological state and therefore the observed measurements. Thus, differences detected in such experimental measurements may not always indicate differences in underlying physiology but can rather indicate differences in confounding variables related to experimental conditions and/or methods.
Mechanistic mathematical models offer a powerful complement to laboratory measurements (5). By accounting for the dynamics of physiological control systems they can provide quantitative assessments of unmeasured variables and associated regulatory mechanisms under given experimental conditions. Here, we analyze high-level phenotype measurements using mathematical models to test hypotheses about the physiological mechanisms underlying differences in the baroreflex system between the SS and SS.13BN rat on high- and low-salt diets.
Baroreflex dysfunction is known to be associated with hypertension in human subjects (7) and has been associated with various rat models of hypertension including the Dahl SS rat and the spontaneously hypertensive (SHR) rat. Baroreflex dysfunction in the SHR model is characterized by both decreased pressure sensitivity (defined as the slope of the curve relating blood pressure to firing rate) and an increased pressure threshold (defined as the minimum blood pressure required for a response) of the afferent baroreceptor nerve fibers (8). It has been shown that this increased pressure threshold (i.e., baroreceptor resetting) is regulated by a mechanism (or mechanisms) independent of vessel distensibility (3). The gain in the overall baroreflex was shown to be blunted by a factor of about four in the SHR compared with the Wistar-Kyoto (WKY) strain (52), but studies of sympathetic baroreflex dysfunction in the SHR are conflicting (19, 52). Of greater interest to the present study are studies of baroreflex dysfunction in the SS rat. Similar to the SHR rat, afferent baroreceptor nerve fiber dysfunction in the SS rat is characterized by elevated pressure thresholds and reduced pressure sensitivity (4). A high-salt diet was shown to elevate pressure threshold but did not alter pressure sensitivity, in SS rats as well as in normal Dahl salt-resistant (SR) rats (2). Defects in both the sympathetic (33) and parasympathetic (61) baroreflex responses have been observed in SS rats on high-salt diets, but normal sympathetic baroreflex responses have been observed in rats on low-salt diets (20).
Kendziorski et al. (25) used a mathematical modeling approach to compare baroreflex function in the SS vs. Brown Norway (BN) rat and found a significant attenuation of baroreflex gain in the SS rat compared with the BN rat. Through linkage analysis, their study was able to map components of baroreflex function to chromosomes 4 and 10. In the present study, differences in baroreflex function between the SS rat strain and a consomic strain, the SS.13BN, are characterized. The BN chromosome 13 has been shown to be associated with a marked reduction in blood-pressure salt sensitivity when introgressed into the SS genetic background (12, 36), but the physiological origins associated with its protective effects are yet to be fully understood. Because of the known associations of baroreflex dysfunction with hypertension along with the marked reduction in blood pressure salt sensitivity associated with the BN chromosome 13, the SS.13BN strain provides a useful control in the present study. Using this design, the goals are twofold: 1) to better characterize baroreflex dysfunction in the SS rat and 2) to identify any association of baroreflex heart rate control parameters with the chromosome 13.
The baroreflex is a nonlinear negative feedback control system known to be responsible for the reflex control of heart rate and believed to be involved in the short-term control of arterial blood pressure (11). Because this reflex plays such a fundamental role in cardiovascular physiology, it has been the subject of intensive investigation over the last several decades which has led to a much better understanding of the mechanisms underlying the various components of the reflex. There remains, however, much to be learned about both individual components and the overall operation of this control system. Most models describing the overall reflex are empirical in nature and ignore many of the known nonlinearities associated with the reflex. To dissect the origins of baroreflex dysfunction in the SS rat we have developed a model based on known physiological mechanisms associated with certain model components and parameterized using data available from the literature. Other model components are parameterized to match data from individual SS and SS.13BN rats on high- and low-salt diets to elucidate functional differences in the reflex system between these four groups. A block diagram of the model is shown in Fig. 1.
Afferent baroreceptor component.
Afferent baroreceptor nerve fibers are strain sensors (21) found in both the aortic arch and the carotid sinus that respond to arterial blood pressure perturbations via changes in firing rate. Both myelinated and unmyelinated fibers are known to exist. Differences in both the static and dynamic characteristics of these two fiber types have been demonstrated (9, 57), but the significance of these differences is unknown. In addition, two general baroreceptor discharge patterns have been identified by Seagard et al. (50), which are summarized by the following characteristics: the type I pattern is characterized by a “discontinuous hyperbolic pattern” with relatively high threshold frequency and saturation firing rate whereas the type II pattern is characterized by a “continuous sigmoidal pattern” with a relatively low threshold frequency and saturation firing rate. Because the impact of fiber type differences in the overall reflex control of heart rate is not well characterized, we do not account for these differences in our model. Furthermore, blood pressure data collected for this study come from aorta only, and we therefore ignore contributions from the carotid sinus baroreceptors. This is not expected to present any major shortcomings given that baroreceptors of the aortic arch have recently been shown to play a dominant role in the baroreflex control of heart rate (42).
The nonlinear features of the baroreceptor response include threshold, saturation, asymmetry, and adaptation. Several basic approaches have been used to capture these important features [see review by Taher et al. (56)]. While the biophysical origins of these nonlinearities are not completely characterized, it is clear that viscoelastic processes of the baroreceptor nerve endings contribute (16). In this study, the modeling approach of Srinivasan et al. (54, 55) is adopted, incorporating refinements proposed by Alfrey (1). The Srinivasan model is based on known mechanisms underlying mechanoreception and divides the afferent baroreceptor into aortic wall, nerve ending, transducer, and encoder components. The aortic wall model used by Srinivasan et al. is based on a complicated pressure-volume relationship derived for cylindrical tubes with elastomeric walls (26). As an alternative to their model, we may represent the static aortic wall using a three-parameter empirical pressure-area relationship proposed by Langewouters et al. (27): (1) where A is the aorta cross-sectional area, Am is the maximal cross-sectional area, Po is the aortic pressure at which the compliance curve reaches a maximum, and P1 represents the steepness of rise, or compliance, of the curve and is called the half-width pressure. This function is used to fit an average pressure-radius curve from a group of SHR rats obtained from data of Andresen (3) and plotted in Fig. 2. Fits of Eq. 1 to the data in Fig. 2, plotted as a solid line in the figure, yield parameter estimates Po = 87.52 mmHg, P1 = 123.76 mmHg, and Am = 28.76 mm2. (Radius R and area A are related by
Following Srinivasan et al. (54, 55), the static expression of Eq. 1 can be used to generate a dynamic model: (2) where Bwall represents the viscosity associated with the aortic wall. Aortic wall strain, εwall, can be computed from aortic radius R according to (3) where Ro is the unstressed radius of the aorta.
Since, over the pressure ranges of interest, the pressure-radius curve shown in Fig. 2 is approximately linear, we chose to simplify the aortic wall component by modeling the pressure-area relationship as (4) where Cwall represents the slope of the pressure-radius curve. The dynamic wall model then becomes (5)
The fit to this simple linear model is shown as a dashed line in Fig. 2, for parameter values Ro = 1.6 mm and Cwall = 0.006 mm/mmHg.
In the Srinivasan model, the aortic wall strain and dynamic strain associated with friction at the nerve ending are uncoupled, and a separate parameter is varied as a function of static input pressure to achieve fits to a variety of data. A more physiological model of baroreceptor nerve-aortic wall coupling based on models of muscle spindles is proposed by Alfrey (1). This model includes a rapid adaptation phase modeled by a single Voigt body and a slower adaptation phase, also known as the rapid resetting phase, modeled by a single dashpot that was allowed to completely reset. However, Brown et al. (8) found that 90% of baroreceptor adaptation to pressure steps occurred within 2–3 s after the maximum firing frequency and that a steady-state level was reached in ∼10–20 s. This suggests that two time constants (and therefore two Voigt bodies) would better account for the rapid adaptation phase. In addition, Munch et al. (38) found that the extent of rapid resetting averaged only 33% and had a time constant of 180–250 s. Thus, the dashpot model of rapid resetting in the Alfrey model (which resets 100%) would be better represented with an additional Voigt body since complete resetting does not occur in the rapid resetting phase. Our proposed model of baroreceptor coupling dynamics is illustrated in Fig. 3 and is described by the following system of ordinary differential equations: (6) where the εs are strains, the Ks are elastic constants, and the Bs are viscous constants. The subscript ne denotes nerve ending, and subscripts 1, 2, and 3 denote various components of the surrounding connective tissue. The strain dynamics governed by Eq. 6 are driven by the wall strain, εwall, which changes in response to changes in pressure according to Eqs. 4 and 5. The strain sensed by the baroreceptor nerve ending is given by (7)
In the Srinivasan model, the baroreceptor threshold and “jump frequency” phenomena are accounted for in the transduction component and are modeled by a set of empirical functions. The encoder model consists of a simple integrate-and-fire mechanism. We have adapted a simpler form of the Srinivasan transducer and encoder components given by the equation below: (8) Here n represents afferent baroreceptor firing rate, S represents baroreceptor strain sensitivity, δth represents strain threshold, and ζ is a parameter introduced to account for the jump frequency at threshold. The parameter ζ is dimensionless with its value constrained between 0 and 1.
Because of the wealth of dynamic baroreceptor response data from the WKY rat available in the literature, we have parameterized a generic (rat strain-independent) model of afferent baroreceptor dynamics using pressure step-response data from this strain. That is, we treat parameters from this component of the model as constants when we analyze data from other rats. This seems reasonable given that Brown et al. (8, 9) were unable to identify any differences in the dynamic baroreceptor response in normotensive WKY rats vs. hypertensive SHR rats. We assume this holds true in the Dahl SS and related consomic strains as well. Model fits to step-response data from Brown et al. (8) and ramp-response data from Andresen (3) are shown in Fig. 4 with associated parameter estimates provided in the legend.
On the other hand, baroreceptor static response parameters, including pressure/strain sensitivity and threshold, have been shown to be quite different between rat strains (2–4). In the present study, we are interested in differences in baroreflex sensitivities and thresholds between the rat strains studied. However, the afferent baroreceptors are not the only components contributing to the sensitivity and threshold properties of the overall baroreflex, and the heart rate data used to parameterize the overall model cannot inform the model of the relative contributions of the various subcomponents to these system properties. We therefore parameterize the sensitivity parameter S of our afferent baroreceptor model using ramp-data available for the SS rat (4) and set the threshold δth to zero. (Threshold is instead accounted for in the peripheral nervous system component of our model.) The full set of parameter values used in the afferent baroreceptor model is listed in Table 1.
Central and peripheral nervous system components.
The least well-understood component of the baroreflex control of heart rate is the integration of baroreceptor afferent signals at the nucleus tractus solitarius (NTS) in the medulla of the brain stem. Important aspects of NTS integration have only recently been discovered. For example, the NTS is now known to be sensitive to many characteristics of the input signal (46, 47), and NTS network architecture is thought to play an important role in determining its response patterns (45). In addition, higher brain centers are thought to play a modulatory role in the integration of sensory signals at the NTS (51, 53, 60) and are thought to continuously modulate baroreflex gain to “optimize the response of the cardiovascular system to daily life challenges” (13, 14). Because of these many complexities, most existing models of the central component of the baroreflex are empirical in nature, many of which are derived through transfer function analysis. These models usually lump the central and peripheral nervous system components together. For instance, Petiot et al. (41) modeled the frequency response of sympathetic nervous activity to aortic depressor nerve stimulation in rat using derivative gain in combination with an overdamped second-order low-pass filter. Their analysis was limited to the linear range of the sympathetic response. Likewise, Kawada et al. (23, 24) successfully modeled the dynamic sympathetic nervous system regulation by the arterial baroreflex using derivative and second order low-pass filters followed by a sigmoidal nonlinearity in rabbit.
For simplicity, we model the central component as an all-pass filter given by (9) where αcns is central nervous system activity and Gcns is central gain. Here, Gcns is set to 1. The peripheral nervous system is modeled by two parallel sigmoidal nonlinearities (one for the sympathetic pathway and the other for the parasympathetic pathway) given by the following 4-parameter logistic functions: (10) where Ts and Tp represent sympathetic and parasympathetic tones, Ts,min and Tp,min represent minimum tones, Ts,max and Tp,max represent maximum tones, and Gs and Gp represent sympathetic and parasympathetic gains, respectively. Here, gain is defined as the steepness of the sympathetic and parasympathetic sigmoidal tone curves. Note that this definition is in contrast to the more standard definition of gain, which is the ratio of the magnitude of an output signal to an input signal (49). Unlike the ratios of the magnitudes of Ts and Tp to αcns, Gs and Gp are operating point independent. The parameters αs,o and αp,o represent baseline or offset sympathetic and parasympathetic activity (defined as horizontal shifts in the sympathetic and parasympathetic sigmoidal tone curves of Eq. 10), respectively, and account for the “central threshold” described by Eckberg (15). A similar approach has been used in previous studies (30, 32).
Since we set the baroreflex strain sensitivity S to a constant value and we model the central component as an all-pass filter, the sensitivities to change in pressure are lumped into the gain parameters, Gs and Gp. We allow Gs and Gp to vary independently so that selective differences in the sympathetic and parasympathetic limbs may be detected. Similarly, because baroreceptor threshold is ignored and assumed to be accounted for by higher level components in the model, αs,o and αp,o are interpreted as representing thresholds. Similar to the gains, we allow these offset parameters to be varied independently. Because the parasympathetic limb has been shown to operate at a higher pressure than the sympathetic limb (17, 52), we constrain αp,o to be greater than αs,o. The parameters involved in our central and peripheral nervous system components are summarized in Table 2.
The baroreflexes exert their chronotropic effects through the reciprocal action of the neurotransmitters norepinephrine and acetylcholine on the sinoatrial node of the heart. Norepinephrine is released from the cardiac sympathetic nerve and acts to accelerate heart rate, whereas acetylcholine is released from the vagus nerve and acts to decelerate heart rate. The secretion rates of norepinephrine and acetylcholine depend on sympathetic and parasympathetic stimulation frequencies, respectively, which are modulated by the baroreflexes via the mechanisms described above. Parasympathetic stimulation is known to exert beat-by-beat control of heart rate, whereas sympathetic control of heart rate is a more gradual process that is unable to exert beat-by-beat control (28, 34, 35).
Mokrane and Nadeau (35) studied the dynamics of the sympathetic heart rate response in dog using a cardiac sympathetic nerve electrical stimulation protocol and developed a model that successfully simulated the dynamic heart rate response to sympathetic stimulation. We assume that sympathetic heart rate dynamics in dog match those in rat and have adapted the Mokrane model in our analysis. Release of norepinephrine from the cardiac sympathetic nerve is given by (11) where cnor is proportional to norepinephrine concentration at the sinus node, Ts is sympathetic tone, qnor is the norepinephrine secretion rate constant, and τnor is the time constant of norepinephrine reuptake processes. The value of τnor = 9.1 s used in our model is taken from Mokrane and Nadeau (35). Levy et al. (28) demonstrated that the secretion rate of norepinephrine is proportional to the rate of reuptake; qnor was therefore assigned the value τnor−1. Sympathetic tone is measured in this model in arbitrary units (AU). Therefore, the variable cnor is measured in arbitrary units as well. (It is possible but not necessary to introduce an additional arbitrary conversion factor to express cnor in concentration units.)
The heart rate response to sympathetic stimulation is known to saturate in the physiological range (6). As in the Mokrane-Nadeau model, this saturation effect is accounted for in the static response of the heart to sympathetic stimulation (ΔHRs,s), and is modeled by the Hill equation with a Hill coefficient of 2: (12) Knor is the level of norepinephrine producing a half-maximal response, and ΔHRs,max is the maximum increase in heart rate due to sympathetic stimulation: (13) where HRmax is the maximum heart rate and HRo is the intrinsic heart rate. The dynamic response of the change in heart rate due to sympathetic stimulation is modeled by the first order process given below (14) where ΔHRs is the dynamic increase in heart rate due to sympathetic stimulation and τHR,nor is the time constant associated with sympathetic heart rate dynamics. The value τHR,nor = 2.1 s is taken from Mokrane and Nadeau (35). The value of Knor = 1.12 AU is estimated using dynamic sympathetic heart rate response data from Warner and Cox (58) as shown in Fig. 5. The values of HRo and HRmax are adjusted to match data from individual rats in the overall model.
In our model, the parasympathetic heart rate response assumes the same form as the sympathetic heart rate response model described above: (15) where cach is proportional to the concentration of acetylcholine at the sinus node, Tp is parasympathetic tone, qach is the acetylcholine secretion rate, and τach is the time constant associated with acetylcholine hydrolysis. The value of τach = 0.2 s used in our model is taken from Levy et al. (28). As in the sympathetic heart rate model, qach is assigned the value τach−1. Parasympathetic tone TP and the variable cach are measured in this model in AU.
There is a known relationship between heart rate variability (HRV) and parasympathetic activity that has been described as “a function in which there is an ascending limb where HRV increases as parasympathetic effect increases until it reaches a plateau level; HRV then decreases as parasympathetic effect increases” (18). This suggests that, like the sympathetic heart rate response, the parasympathetic heart rate response saturates in the physiological range. Using data from Warner and Cox (58) (Fig. 6A), we found that a function of identical form to Eq. 12 may be used to characterize the static heart rate response to parasympathetic stimulation: (16) Here, Kach is the level of acetylcholine producing a half-maximal response, and ΔHRp,max is the maximum decrease in heart rate due to parasympathetic stimulation, and is given by (17) where HRmin is the minimum heart rate and HRo is the intrinsic heart rate.
Mokrane et al. (34) used a transfer function to characterize the dynamic heart rate response to vagal stimulation. The vagal heart rate response was characterized by a combination of two different filter behaviors: a low-pass filter with a mean time constant of 2.5 s and an all-pass filter. These two filter behaviors are likely due to acetylcholine acting through both a fast G protein-coupled pathway as well as a slower second-messenger pathway (28). We account for these two filter behaviors using an algebraic and a first order differential equation: (18) where ΔHRp,fast and ΔHRp,slow represent the parasympathetic mediated decreases in heart rate due to fast and slow pathways, γ is the proportion of the dynamic response acting through the fast pathway, and τHR,ach is the time constant associated with the slow pathway (assigned the value 2.5 s). The overall parasympathetic heart rate response, then, is given by the sum of the individual fast and slow pathways as shown below: (19) The time-dependent response of HR (from Ref. 58) to the dynamic parasympathetic stimulation protocol shown in Fig. 6B allows us to estimate values for Kach and γ of 0.65 AU and 0.75, respectively. Values for the parameters HRo and HRmin are obtained on an individual basis from experiments in the present study.
It has long been known that combined sympathetic and parasympathetic stimulation does not produce an algebraically additive effect at the sinoatrial node. Combined autonomic effects are instead complicated by an interaction between the two pathways (29). Rosenblueth and Simeone (48) accounted for this interaction using a model that multiplied two factors, one for each autonomic pathway, concluding that the sympathetic and parasympathetic pathways exerted their effects independently of one another. Warner and Russell (59) proposed a different model, dismissed the findings of Rosenblueth and Simeone, and determined that the two pathways do indeed interact. Katona et al. (22) later recognized that the two models were essentially the same and that the different conclusions drawn from the models were a matter of interpretation. Here, we adapt the model of Warner and Russell to account for the combined sympathetic and parasympathetic effects on heart rate: (20) where HRs and HRp are given by (21) In our model, β is an adjustable parameter constrained between 0 and 1.
Parameters used in our heart rate effector component are summarized in Table 3. Values of fixed parameters are indicated in the table. Parameters identified as adjustable are estimated on an individual basis in the analysis below. Yet, before individual data are analyzed, the heart rate effector model may be validated based on independent data from Warner and Russell (59). Figure 7A shows a dynamic stimulation protocol combining vagal and sympathetic stimulation. The associated measured and model-predicted heart rate data are shown in Fig. 7B. These model predictions, which involve no additional parameter adjustments (besides the adjustable parameters β, HRo, HRmax, and HRmin), demonstrate that the model is able to accurately simulate the major features of the heart rate response to combined stimulation.
Generation of consomic population.
The consomic rat line was derived from inbred normotensive BN/SsNHsd/Mcw (BN) rats and salt-sensitive hypertensive Dahl SS/JrHsdMcwi (SS) in which the chromosome 13 of the BN was introgressed into the background of the SS as we have described previously (12).
Chronic phenotyping protocol.
Blood pressure was measured by radiotelemetry using a Dataquest ART 3.1 system (Data Sciences International; DSI, St. Paul, MN). A gel-filled catheter connected to a transmitter (model TA11PA-C40, DSI) was surgically implanted in the femoral artery with the transmitter anchored under the skin over the flank area of the rat. Rats used were adult, male and female of the SS or the consomic SS.13BN strain obtained from in-house colonies. These animals were maintained from weaning on a custom AIN-76 diet (Dyets, Bethlehem, PA) containing 0.4% NaCl with water provided ad libitum. Transmitters were implanted at 9 wk of age, and the animals were allowed to recover for 1 wk before pressures were collected for 3 days to ensure a steady state. The diet was then switched from 0.4% NaCl to 8% NaCl for 2 wk, and pressures again were collected for analysis. Pressure was collected at each salt level for 2 min at 100 Hz between 9 AM and 12 noon. All measurements were made in a sound attenuated room where rats were maintained in their home cage, placed on the receiver for the duration of the study. All protocols were approved by the Institutional Animal Care and Use Committee.
Because heart rate was not obtained in the above phenotyping protocol, experimental heart rate data used to parameterize the model was extracted from the arterial blood pressure time series data. This was done by first high-pass filtering the time series data and then detecting the time of zero-crossings of the pulse-pressure upstroke. The high-pass filter used was of the FIR class and from the Kaiser family of window functions. The fundamental frequency of the arterial pressure pulse corresponds to heart rate. By filtering frequency components lower than this fundamental frequency, both the arterial pressure offset as well as any slow oscillations are eliminated so that the pulse pressure is made to cross zero, allowing us to detect the time of zero-crossings of the pulse-pressure upstroke. We assume that the rat heart rate does not go <240 beats/min (4 Hz) under normal baseline physiological conditions. Therefore, the stop-band frequency of the Kaiser window was specified as 3 Hz and the pass band as 4 Hz. The stop-band attenuation and pass-band ripple were specified as 1 × 10−3 and 5.75 × 10−2, respectively. The filter order and coefficients required to meet these specifications were calculated using MATLAB's Signal Processing Toolbox. The filter was applied to the data using MATLAB's filtfilt function. This function performs zero-phase digital filtering by processing the input data in both the forward and reverse directions; thus, the effects of filter phase distortion are minimized. Zero-crossings were detected using conditional statements in an iteration loop.
Noise in the arterial pressure recordings occasionally resulted in relatively large spikes appearing in the heart rate data. We assumed that, under normal physiological conditions, instantaneous changes in heart rate cannot exceed 50 beats/min, and we therefore removed data from any spikes exceeding this upper limit.
The overall model of baroreflex heart rate regulation is comprised of a system of eight ordinary differential equations. The differential equations are solved using MATLAB's solver ode23, an implementation of an explicit Runge Kutta 3(2) pair. The model contains a total of 13 adjustable parameters that are estimated for individual rats using the pressure heart rate data derived as described above. [In addition, three initial conditions ε2(0), ε3(0), and cnor(0) were estimated for each time series to compare model simulations to the data.] The parameters were estimated by fitting the model to the heart rate data using a genetic algorithm and least squares objective criterion given by (22) where mse is the mean square error, Nd is the number of data points, HRdata is experimental heart rate at time ti, and HRmodel is model predicted heart rate at time ti. A single time course contained insufficient data to precisely identify all 13 model parameters, so we therefore used multiple datasets (n = 6 from the SS.13BN strain and n = 9 from the SS strain for both dietary conditions) from a number of individual rats from each of the four experimental groups studied to determine probability distributions and associated confidence bounds on parameter values to explain the data. Because initial conditions on variables with short time constants were arbitrarily selected, the objective criterion of Eq. 22 was computed only for t ≥ 10 s.
Parameter bounds were similarly set so that all possible parameter values fell within reasonably realistic physiological constraints as determined by available published data. For example, upper bounds of tone parameters were set so that the saturation effects of tone on heart rate (6, 18) could be reproduced without going to unrealistically high values. Similarly, upper and lower bounds on minimum, maximum, and intrinsic heart rate values were set using data from Bolter and Atkinson (6) and Nylander et al. (39). Upper bounds on strain initial conditions could be calculated using the pressure-radius relationship derived for rat aortas (see Fig. 2 and Table 1) and baroreceptor mechanics model parameters (see Fig. 3 and Table 1) for the highest expected arterial pressure estimated from the data. Where appropriate data was not available, upper and lower bounds were set to cover the broadest range of possible values that can be reasonably expected from initial optimization studies. All parameters were constrained to positive values. A list of initial parameter values, initial conditions, and upper and lower bounds, is given in Table 4. These initial values and bounds were used in the optimizations of each dataset from each experimental group of rats.
Because of the several nonlinear interactions between various components of our overall baroreflex model, typical gradient-based search methods could not be used. Such methods generally fail when applied to optimization problems involving nonsmooth nonconvex parameter landscapes. We instead relied on the stochastic search algorithms to fit our model to the experimental data. Unlike gradient-based methods, stochastic search algorithms such as simulated annealing and genetic algorithms are able to escape local minima and are ideally suited for these types of nonlinear global optimization problems. Genetic algorithms, in particular, are a powerful class of optimizers that can easily be parallelized. To fit the heart rate time course data, we have implemented a genetic algorithm based on PIKAIA (10) in the MATLAB environment. PIKAIA is an adaptive genetic algorithm employing decimal-encoding, crossover, elitism, fitness-based roulette-wheel selection, and an adaptive mutation rate. All optimizations in the present study were performed using our implementation of PIKAIA on a 128-node computational cluster.
Using aortic blood pressure as a model input, we fit model predicted heart rate to experimental heart rate with a genetic algorithm nonlinear optimization protocol and the least-squares objective criterion. The ability of our model to reproduce important features of the experimental heart rate time series data in the two strains of rat studied (results shown are from rats on low-salt diets) is illustrated by the examples shown in Fig. 8. Figure 8 demonstrates that the model is able to capture both the large relatively slow peaks (Fig. 8, C and D) as well as the relatively small but rapid oscillations (Fig. 8, E and F). We note, however, that the results shown in Fig. 8 represent cases where the model is able to reproduce these rapid oscillations particularly well and do not necessarily reflect model performance for all datasets. In some cases, the model was unable to reproduce these rapid transients, and the optimizer tended produce smooth fits through the data to minimize the sum of square error of the objective criterion. In most cases the model is able to fit the large slow peaks reasonably well. The fact that the model was able to capture both the slow and fast heart rate variations is striking given that most all parameters related to baroreflex dynamics (i.e., time constants) were fixed based on data from a different species.
The model was used to assess baroreflex function in four different experimental groups of rats including two genetic strains, the SS and SS.13BN, on both low- and high-salt diets. The model was fit to a total of 12 datasets from the SS.13BN rats under each of the two experimental conditions (6 datasets for each condition), and a total of 18 datasets from SS rats under each of the two experimental conditions (9 datasets for each condition). The model-predicted heart rate time-series data (associated with the optimized parameter sets) is plotted along with experimentally derived heart rate data for each of the 12 datasets of the SS.13BN rats in Fig. 9 and for each of the 18 datasets of the SS rats in Fig. 10. As illustrated in Figs. 9 and 10, the model was able to capture most of the important features of each dataset. In some cases, the model tended to underestimate the magnitude of certain peaks or missed certain peaks entirely; however, in nearly every case the model was able to capture the general pattern of the data. A notable exception in which the model fails to capture the pattern of the data entirely is shown in Fig. 10B, dataset 4. The reasons for the poor fit of this dataset are not entirely clear.
Because the model explicitly accounts for the various subcomponents comprising the overall baroreflex, intermediate variables and parameters of interest can be estimated from the information contained in the data. The parameters estimated from the individual time courses are listed in Table 5. These include estimates of important quantities such as sympathetic and parasympathetic gains and offsets. The precision of parameter estimates is limited both by the amount of information contained in the heart rate data, which varies from dataset to dataset, as well as by the amount of interindividual variability of physiological parameters in each group (population) of rats. The averages and standard deviations of parameter estimates for each experimental group of rats, along with the optimal set of parameter values used to fit each individual dataset, are summarized in Table 5. Experimental mean arterial blood pressure and heart rate phenotype data for each group of rats is also given.
Figure 11 shows model-predicted (unmeasured) afferent baroreceptor firing rate, autonomic tone, and neurotransmitter concentrations associated with dataset 2 of Fig. 9A. On the basis of model predictions for all individuals simulated, we may compare the levels of sympathetic and parasympathetic outflows in the baseline resting state between individuals and between experimental groups. Doing this, we assume that the telemetry recordings are taken from rats in a resting state. That is, it is assumed that no stressful stimuli are affecting the measurements. The predicted baseline autonomic outflows of each rat along with group means and standard deviations are summarized in Table 5.
To determine whether there are any statistically significant differences in mean parameter values between the SS and SS.13BN rats on either high- or low-salt diet, we performed a one-way ANOVA on parameter estimates from each group of rats. To determine specifically which pairs of means differ statistically, we followed the ANOVA with the Dunn-Sidak multiple-comparisons test. This statistical analysis was also applied to group means of phenotype data (mean arterial blood pressure and heart rate) as well as model predictions (sympathetic and parasympathetic tones). Two outliers were removed prior to performing any statistics on the data. Here, an outlier was defined as a data point that is >1.5 times the interquartile range above the 75th or below the 25th percentile of the sample (the outliers are indicated in Table 5 with an asterisk). All parameter values associated with a given outlier dataset were excluded from statistical analysis, even if other parameters associated with the dataset were not outliers. One of the outliers was associated with dataset 4 of the SS high-salt group; since the model was not able to fit this dataset, it was not surprising that its associated parameter values were outliers. The reason for the other outlier (dataset 5 of the SS.13BN high-salt group) excluded from the analysis was less clear. The model appeared to produce a reasonable fit to the data; however, several important parameter estimates (most notably, the sympathetic and parasympathetic gains) of this dataset were quite different from estimates of other datasets of the same group.
Pairs of groups showing statistically significant differences are indicated symbolically in Table 5. Important results from the analysis are summarized in Fig. 12. In terms of experimental phenotype data, statistically significant differences in mean arterial pressure were detected between the SS rat on a high-salt diet and every other group as well as between the SS.13BN high- and low-salt diet groups. No statistical differences were detected in mean heart rate between any of the four groups. Previously documented experimental findings in SR and SS rats suggest differences in the sympathetic and parasympathetic gain and offset parameters, Gs, Gp, αs,o, and αp,o, between these strains (2, 4, 20, 33, 61). Here, statistical differences in estimates of αs,o were detected between the SS high-salt group and SS low-salt group, and between the SS high-salt group and the SS.13BN low-salt group. There was also a trend toward a greater αs,o in the SS.13BN high-salt group compared with the two low-salt groups, but this was not found to be significant. The αp,o parameter demonstrated a similar trend with an increase with high-salt diet in both strains. There also appeared to be a difference in αp,o between the two strains with higher values in the SS strain; however, the only statistically significant difference detected in αp,o was between the SS high-salt and SS.13BN low-salt groups. These model-based predictions that a high-salt diet elevates sympathetic and parasympathetic offset in the SS rat and that there exists a defect in parasympathetic offset in SS compared with SS.13BN rats require further experimental investigation.
The gain parameters also revealed interesting physiological differences between the different groups of rats. Our results show a mild attenuation in Gs in SS.13BN rats on a high-salt diet and a much larger attenuation of Gs in SS rats on a high-salt diet. The SS rats on a low-salt diet were found to have the highest sympathetic gain (Gs) among all four groups, and statistically significant differences between the SS low-salt group and the two high-salt groups were detected. It is interesting that SS, and not the SS.13BN rat, which is known to be protected from developing hypertension, has the highest gain on a low-salt diet. No statistical differences were detected in the parasympathetic gain (Gp) between any of the four groups. This suggests that the central/peripheral nervous systems may play some role in the baroreflex dysfunction observed in the SS rat on a high-salt diet. Specific conclusions drawn from these data are that a high-salt diet selectively attenuates the sympathetic gain and that SS.13BN rats are protected from this attenuation. There may be an intrinsic difference in sympathetic gain that is associated with a lower gain in the SS.13BN rats, but this finding is not significant. One parameter we were surprised to detect a significant difference in was the minimum heart rate, which was found to be significantly lower in the SS.13BN high-salt group compared with both SS groups. It is unclear why this parameter is lower in this particular group and what the implications of this difference may be. Again, these model-based predictions require experimental verification.
While model predictions failed to reveal any statistically significant differences in sympathetic or parasympathetic outflows between any of the four groups, an interesting trend emerges. Figure 12G reveals that sympathetic tone (Ts) is highest in the SS rat on a high-salt diet. Furthermore, a high-salt diet was shown to elevate Ts in the SS strain and to reduce Ts in the SS.13BN strain. On the other hand, as shown in Fig. 12H, parasympathetic tone (Tp) was higher in the SS strain on both high- and low-salt diet, and a high-salt diet elevated Tp in both strains. These results are in line with the finding of a selective sympathetic defect in the SS rat. It would seem that the SS.13BN rat is able to compensate any potential increases in blood pressure that would otherwise be induced by a high-salt diet through sympathetic withdrawal. The SS rat, on the other hand, is unable to exert this sympathetic withdrawal compensation, possibly because of the observed defects in sympathetic gain with high-salt diet, which therefore manifests as an increased sympathetic tone with a high-salt diet. This increase in sympathetic tone may then contribute to the increased blood-pressure phenotype observed in this strain with a high-salt diet. On the other hand, because parasympathetic gains are identical between the two strains and unaffected by diet, the parasympathetic outflows are effectively increased with high-salt diet in both strains as would be expected.
Using a mechanistic modeling approach, we are able to explain a number of previously published experimental findings related to baroreflex dysfunction in the Dahl SS rat. Furthermore, we have uncovered a number of important differences in baroreflex physiology between the SS and consomic SS.13BN rats on high- and low-salt diets and have for the first time associated chromosome 13 of the Dahl rat with dysfunction of the baroreflex control system. Studies have demonstrated a sympathetic defect in SS rats on high-salt diets (33) but not on low-salt diets (20). This defect with high-salt diet can be explained by the attenuated sympathetic gain and therefore blunted sympathetic withdrawal in the SS rat on high-salt diet as demonstrated in this study. Furthermore, previous studies have demonstrated an increased afferent baroreceptor pressure threshold in both the SS and SR rat on a high-salt diet (2). It is difficult to make any conclusions related directly to the afferent baroreceptor component in this study because parameters of that component were fixed in our analysis. However, the sympathetic and parasympathetic offset parameters were adjustable and were informed by experimental data, and these parameters reflect composite differences in thresholds of both the afferent and central components of the reflex. We found increased sympathetic offset with high-salt diet in the SS rat which may be related to an increased afferent baroreceptor pressure threshold as described previously.
An important difference in SS and SS.13BN physiology uncovered through our analysis includes a markedly reduced attenuation of sympathetic gain with high-salt diet in the SS.13BN rat. In addition, model predictions of sympathetic and parasympathetic outflows revealed an interesting trend that we used to develop a new hypothesis related to baroreflex dysfunction in the SS rat. Model analysis predicts that high-salt diet reduces sympathetic tone in the SS.13BN rats, whereas it causes elevation of sympathetic tone in the SS rat. If real, this apparent defect in sympathetic withdrawal in the SS rat on high-salt diet may be due to a defect in sympathetic gain. Thus, a defect in sympathetic gain, leading to increased sympathetic outflow, may be playing a role in the blood pressure phenotype of the SS rat on a high-salt diet. This hypothesis requires further testing by collecting a larger sample size to increase the power of our analysis and also by collecting direct measurements of autonomic outflows and validating our model against these measurements.
This hypothesis implicates at least some neural contribution of chromosome 13 to the salt-sensitive hypertension phenotype of the Dahl SS rat. Though originally hypothesized to be primarily a defect in renal and humoral factors, neural mechanisms were later explored and shown to play an important role in the pathogenesis of the disease (31, 37). We believe that Dahl SS hypertension cannot be explained by a single mechanism or defect, and that renal, humoral, endocrine, neural, and hemodynamic factors are all involved in some way (31, 37, 43). However, the relative contributions of each of these factors and their genetic associations are not well known. A major advantage of the modeling approach in general is the ability to integrate the effects of multiple physiological controllers on whole system behavior so that the relative roles of the various controllers in salt-sensitive hypertension may be better elucidated. The baroreflex model developed here represents an important starting point in this effort.
Though we were able to confirm a number of previous findings, our results conflict with some previously published findings. For example, Whitescarver et al. (61) concluded that baroreflex dysfunction in SS rats is due to a parasympathetic defect. However, our results failed to identify any differences in parasympathetic physiology between any of the strains studied. In addition, studies have shown a defect in afferent baroreceptor sensitivities in SS rats, even on low-salt diet (4, 33). We were unable to definitively identify any dysfunction in baroreflex gain in the SS rat on low-salt diet. There are many possible explanations for these discrepancies. For example, all of these previous studies used the SR rat as a control, which has very different genetics and possibly very different physiology than the consomic SS.13BN control rat used in the present study. Thus, differences between SS and SR rat may not exist between the SS and SS.13BN rat. (Note that the SS.13BN retains some of the salt-sensitivity of the SS strain.) In addition, the measurements performed in these studies often relied on pharmacological manipulations or invasive surgical procedures involving anesthesia. Such perturbations may themselves be expected to alter the function of the physiological control systems of interest or, at the very least, may alter the physiological state of the animal. Thus, the findings may depend on the experimental conditions and may not necessarily hold in the animals' baseline resting state. Our experimental approach relied on a telemetry protocol to minimize the effects of stressors on our measurements. In addition, the parameters of our model are, theoretically, less dependent on operating points since they explicitly account for many of the nonlinearities of the control systems of interest. In this case, the model parameters would be expected to hold under a wider variety of experimental conditions.
Despite that fact, it is crucial to note that the precision of our parameter estimates is highly dependent on the sensitivity of our model to those parameters. This sensitivity is dependent on both the data used to estimate the parameters as well as the structure of the model itself. Concerning the structure of the model, we have accounted for a number of known nonlinear features of the baroreflex system that tend to increase the complexity of the parameter landscape. Nonlinearities such as saturation can lead to a number of local minima in the solution space. We have addressed these difficulties by using optimizers ideally suited to overcome these local minima, but the likely possibility exists that our parameter estimates do not correspond to the global minima. With regard to the data, it may be difficult to obtain sensitive estimates of all parameters using heart rate data alone. Heart rate datasets containing a number of large peaks and valleys are rich in information pertaining to the baroreflexes and would be expected to confer a higher sensitivity of the model to the parameters, whereas a dataset containing a relatively constant heart rate contains less information about the baroreflexes and could potentially be explained by several alternative combinations of parameter values (and therefore, would confer less sensitivity of the estimates to the data). Most previous modeling studies of the overall baroreflex addressed this problem by perturbing the system in some way such as a sit-stand maneuver (40) or the Valsalva maneuver (30). Pharmacological interventions such as phenylephrine and sodium nitroprusside have also been used as a means of perturbing the system (20, 33); such interventions allow the response of the baroreflex system to be assessed under extreme physiological conditions and offer a means of challenging the model at these extremes. In the present study, we identify parameters using only the naturally occurring spontaneous fluctuations in heart rate. The current model can be used to design experiments using dynamic perturbations to improve the precision of parameter estimates.
In addition to the physiological insights into Dahl salt-sensitive hypertension gained through our modeling analysis as described above, the model developed here can serve as a powerful tool for analyzing and integrating data from different components of the baroreflex and may be valuable in helping us better understand the general function of and interaction between various components of the overall reflex. As far as we are aware, this work represents the most mechanistically detailed available model of the overall baroreflex heart rate control system in the rat. Ideally, parameters of our model would be identified with detailed measurements of each of the components of the baroreflex for each of our experimental rat strains. Although these data are not available for the rat strains used in the present study, such detailed measurements exist from other strains and from other species. With the assumption that there are no differences between strains or species in certain model parameters, we are able to take advantage of these available data. Specifically, we have demonstrated that baroreceptor dynamics parameters from the WKY rat and heart rate dynamics parameters from mongrel dogs allow us to effectively simulate heart rate data from SS and SS.13BN rats. We expect that many of the parameter values used in the present model are conserved across species; it would therefore be interesting to see whether the model is able to explain differences in the baroreflex control system in different animal models of baroreflex dysfunction or in human studies.
A remarkable feature of these fits is their ability to capture both the slow and fast variations in heart rate as illustrated in Fig. 8. These two types of variation are believed to be mediated separately by the sympathetic and parasympathetic nervous systems and thought to be a consequence of the different time constants associated with the two nervous systems. It is believed that the shorter time constants of the parasympathetic pathway allow it to exert beat-by-beat control, while the longer time constants associated with the sympathetic pathway are associated with relatively smoother and slower variations in heart rate (28, 34, 35). By accounting for these differences in time constants, we are able to simulate these disparate effects of the two pathways on heart rate.
However, the model is unable to capture the fast dynamics of the parasympathetic nervous system in all of the datasets. We can speculate a variety of explanations for this, the simplest being that there is some variation in either baroreceptor or PNS/sinus node dynamics among different individuals and that additional parameters must be made adjustable to capture these dynamics. The other explanation is that the smoother fits through the data represent local minima in the solution space. This would also explain why the parameters estimates of the parasympathetic nervous system were less precise than those of the sympathetic nervous system and may be a reason why we were unable to detect differences in parameters of the parasympathetic nervous system between the different rat strains studied. Furthermore, our model neglects the delays associated with both the sympathetic and parasympathetic nervous systems because they were assumed to be negligibly small in the rat (28, 41). However, even small delays in the parasympathetic nervous system could greatly affect the phase of the fast heart rate variations. If these variations were out of phase, the smooth fits observed with many datasets could represent actual global minima. Finally, even in those datasets for which the model is able to capture the fast variations, the model was never able to perfectly capture all of the data. It is not clear whether peaks not captured by the model represent important physiological processes, physiological noise, or noise associated with the processing of the heart rate data from the blood pressure data. It will be interesting to see whether better fits of these fast dynamics can be achieved by accounting for system delays to correct for phase differences, by analyzing heart rate in the frequency domain, or by attaining actual ECG heart rate measurements to eliminate any noise introduced through our signal processing protocol.
The performance of the model in capturing the smoother slower peaks associated with the sympathetic nervous system was generally better. This may be reflected in the higher precision of the sympathetic gain estimates compared with the parasympathetic estimates. However, although the model was able to capture the general pattern of most of these peaks, the magnitude of the peaks was often underestimated. In addition, not all of the peaks could be simulated by the model. This suggests that the model may be missing some important dynamic of the baroreflex control of heart rate, or that some other physiological control mechanism(s) is at play. We have, for example, made several simplifications in the baroreflex control system such as neglecting contributions from unmyelinated fibers. Brown et al. (9) demonstrated differences in frequency responses between myelinated and unmyelinated fibers. In addition, Seagard et al. (50) identified two types of baroreceptor responses (type I and type II), but the model accounts for only type I response. Another simplification of our model was the assumption that the CNS acts as an all-pass filter. Many studies have demonstrated that this is not the case and suggest that filter behavior at the CNS may play an important role in baroreflex dynamics (23, 24, 41, 46, 47). The model also ignores contributions from the baroreceptors in the carotid sinus. Pressures in the carotid sinus can be quite different from those in the aortic arch, and accounting for the relative contributions of the baroreceptors of the carotid sinus may be required to achieve better fits. Finally, many other physiological controllers are known to be involved in the autonomic reflexes including the chemoreceptors, cardiopulmonary receptors, and others. It is therefore unrealistic to expect to simulate all features of the heart rate data using a model of only the baroreflex system. The model developed in this study may easily be adapted to include contributions from other fiber types and from carotid sinus baroreceptors. In addition, the model serves as a module that may be incorporated into a more comprehensive cardiovascular model. In this manner, the relative roles of different physiological systems in the control of heart rate may be dissected more completely, and we may begin to elucidate the relative roles and interactions between the various components of arterial blood pressure regulation so that a deeper understanding of the physiological origins of blood pressure salt-sensitivity becomes within closer reach.
This work was supported by National Heart, Lung, and Blood Institute Grant HL-82798.
No conflicts of interest, financial or otherwise, are declared by the authors.
We are grateful to Terry Kurth for technical assistance with the chronic phenotyping protocol.
We are grateful to The Milwaukee Institute for providing resources for high-performance computation.
- Copyright © 2010 the American Physiological Society