In this study we develop a modeling framework for predicting baroreceptor firing rate as a function of blood pressure. We test models within this framework both quantitatively and qualitatively using data from rats. The models describe three components: arterial wall deformation, stimulation of mechanoreceptors located in the BR nerve-endings, and modulation of the action potential frequency. The three sub-systems are modeled individually following well-established biological principles. The first submodel, predicting arterial wall deformation, uses blood pressure as an input and outputs circumferential strain. The mechanoreceptor stimulation model, uses circumferential strain as an input, predicting receptor deformation as an output. Finally, the neural model takes receptor deformation as an input predicting the BR firing rate as an output. Our results show that nonlinear dependence of firing rate on pressure can be accounted for by taking into account the nonlinear elastic properties of the artery wall. This was observed when testing the models using multiple experiments with a single set of parameters. We find that to model the response to a square pressure stimulus, giving rise to post-excitatory depression, it is necessary to include an integrate-and-fire model, which allows the firing rate to cease when the stimulus falls below a given threshold. We show that our modeling framework in combination with sensitivity analysis and parameter estimation can be used to test and compare models. Finally, we demonstrate that our preferred model can exhibit all known dynamics and that it is advantageous to combine qualitative and quantitative analysis methods.
The main role of the cardiovascular (CV) system is to provide adequate oxygenation of all tissues, a function which is achieved by maintaining homeostasis of blood flow and pressure. When a mammal is subjected to an orthostatic maneuver (e. g. , running, jumping, etc.), its blood volume is redistributed, moving the system state away from homeostasis [1]. To re-establish homeostasis a number of control mechanisms are activated regulating vascular resistance and compliance, and cardiac pumping efficiency and frequency. An important contributor to this control system is the baroreflex, which uses specialized neurons called baroreceptors (BR) for signaling [2]. The BR neurons originate in the arterial wall and terminate in the nucleus solitary tract (NTS), where sensory information is integrated. These neurons are continuously stimulated via activation/inhibition of mechanosensitive receptors responding to changes in arterial wall stretch imposed by pulsating blood pressure [3]. This stimulus modulates the formation of action potentials propagating along the BR nerves terminating in the NTS, where efferent signals are generated to regulate heart rate, cardiac contractility, as well as vascular resistance and compliance. It is known that the baroreflex system contributes to short-term blood pressure regulation, operating on a time-scale of seconds to minutes [4]. For example, upon head-up tilt, blood is pooled in the lower extremities, increasing blood pressure in the lower body, while decreasing it in the upper body, causing an imbalance, which persists until the baroreflex system is activated. Figure 1 shows a schematic representation of the baroreflex pathways. While the BR pathways are generally well established, analysis of the complete control system, including afferent and efferent signaling, is hindered by the difficulty of measuring the activity of each component without disrupting the feedback loop. For example, in vivo, only macroscopic quantities can be measured non-invasively including heart rate and blood pressure. From such measurements it is difficult to examine how the individual components of the system interact and consequently it is difficult to determine which sub-systems are compromised in subjects experiencing baroreflex failure [5] or decreased arterial baroreflex sensitivity [6]. These difficulties limit the development of targeted diagnosis procedures and treatment plans aiming to alleviate symptoms for patients. Mathematical modeling is an eminent tool for gaining more insight into this complex feedback loop, offering a stringent and systematic way to identify underlying mechanisms of the system. For example, the only way to estimate model parameters and thereby suggest essential biomarkers, which may not be directly measurable, is by using models in combination with direct measurements. Modeling also offers a way to understand complex systems, as it makes the inaccessible accessible, a concept denoted the “mathematical microscope” [7]. This paper focuses solely on the afferent part of the baroreflex system, while future studies will address efferent signaling and integration of the two parts within a system level model. Since the 1950s researchers have put forward numerous mathematical models [8]–[19], which tried to integrate known dynamics with hypothesized mechanisms in order to provide more understanding of the system as a whole. Many insights have been gained, however, most of these models were developed to describe BR response to a particular stimulus, rather than to a range of stimuli eliciting all known responses. Therefore they all lead to different hypotheses explaining the system mechanisms. Inspired by shortcomings of previous studies, we developed a modeling framework containing model components reflecting physiological pathways. This framework splits the afferent signaling into three parts describing vessel wall deformation, mechanoreceptor stimulation, and the frequency of action potential generation. For each component we propose multiple models, which we test both qualitatively and quantitatively. This new approach allows us to understand the contribution of each component model to the overall signal. For example, if the objective is to build a BR model that can reflect the response to a sinusoidal pressure stimulus observed experimentally, the modeling framework can be used to identify which combinations of components are sufficient to describe the experimental outcome, and which component models may be excluded from possible explanations of observed features. Moreover, we show how our framework may be used to inform hypotheses, by suggesting a particular component mechanism responsible for generating a given pressure-response feature of BR firing.
In this section we describe the main qualitative characteristics of BR firing rate as well as the data used for quantitative model tests. To model the dynamics, which produce the BR firing rate in response to given blood pressure stimuli, we include three components separating distinct physiological pathways, and for each component we develop a number of linear and nonlinear models. The three components (Figure 3) include: arterial wall deformation, mechanoreceptor stimulation, and action potential generation. As a driving force for the models we use arterial pressure, which determines arterial wall deformation quantified by the wall strain. The wall deformation stimulates the stretch sensitive mechanoreceptors found in the BR nerve endings within the arterial wall. Thus changes in blood pressure modulate the opening of these channels, and thereby the current flowing through them, which determine the rate at which action potentials are formed. The time between subsequent action potentials determines the firing rate, and thus our models relate the receptor strain to the frequency of action potentials, thereby allowing us to predict the BR firing rate. For each model component, described below, we review previous modeling methodologies and use these to inform the design of the new component models, collectively used to describe the firing rate of afferent BR neurons in response to an applied blood pressure stimulus.
Models will be tested quantitatively using three types of pressure stimuli: sinusoidal at a fixed frequency, a step-increase, and a step-increase followed by a step decrease (Figure 2A–C). We investigated six linear and nonlinear models summarized in Table 3. For the wall strain three models were investigated, the simplest assumes the wall strain has a spring-like response (denoted). The second model (denoted) accounts sigmoidally for increased stiffening with increased pressure, and finally we investigate a viscoelastic model (). The mechanoreceptor strain, is modeled using one, two, and three Voigt bodies, respectively, in series with the spring (). Finally, two models were used for determining the BR firing rate, a linear model () and an integrate-and-fire model (). As mentioned above, these models can all be described as a system of algebraic and differential equations. For all models the model input is pressure and the model output is BR firing rate, initial conditions were computed to ensure that model solutions start at steady state. The objective was to estimate model parameters minimizing the least squares error between the model and data. This is calculated from the point wise residual error between model and datawhere is the average firing rate of the specific data set considered and denotes the parameter vector. To estimate the parameters we minimize the sum of squares cost function (referred to as RMSE in Tables 4,5, and 6) Since data is only available for the BR firing rate and the pressure stimuli, for most models not all parameters are identifiable. We denote as identifiable parameters, those that are sensitive and not correlated, given the model output and the associated available data [65]. In this study, identifiability of parameters was determined using sensitivity based methods [66]. Subsequently, for models completely characterized by smooth functions, the Levenberg-Marquardt method was used to estimate model parameters, while for models not fulfilling this requirement (the integrate-and-fire models), parameters were estimated using the Nelder-Mead method. Both used optimization algorithms from Kelley [67]. Below we first describe the methodology used for sensitivity analysis and parameter identification and subsequently we discuss results obtained using nonlinear optimization, the latter is separated according to the input stimulus. Sensitivity analysis: For any smooth model of the form (23), the sensitivities [68]–[70] can be computed asFollowing Pope et al. [71], we use a finite difference approximation to compute where is the unit vector in the component direction and is a small number. The BR firing rate is obtained computationally, with an integration tolerance of imposed on solution of the differential equations, thus is bounded by. To satisfy this requirement we let. Sensitivities are ranked by averaging time-varying functions using the two-norm. For each model, this ranking was used to separate parameters into two groups: one group consisted of parameters to which the model output was sensitive, and the other group consisted of parameters to which the model output was insensitive. Estimating only sensitive parameters allows more reliable estimation of parameters [72]. Not all sensitive parameters are practically identifiable [65], [66]. To identify parameter correlations, we used the QR-SVD subset selection method [71], [73], [74]. We also used a method based on covariance analysis to identify pairs of correlated parameters [66]. For each pair of correlated parameters the least sensitive parameter was kept fixed at its nominal value while the other was included in the subset. Parameter correlations were computed fromwhere is the variance of the assumed noise in the data, is the covariance matrix, and is the correlation coefficient. Parameters for which are labeled as correlated. For the models studied in this work we let. Once a set of uncorrelated sensitive parameters were identified, we used either the Levenberg-Marquardt or the Nelder-Mead method to estimate the subset of practically identifiable model parameters [67]. The Levenberg-Marquardt method was used for models that can be described using smooth functions, while the Nelder-Mead method was used for models including the leaky integrate-and-fire component. Since this model contains a discontinuity the gradient based Levenberg-Marquardt method is not applicable. Sinusoidal stimulus: Now we present results obtained using sinusoidal forcing allowing us investigate asymmetry of the model response. Results (Figure 6) show BR firing rate as a function of time and BR firing rate as a function of stimulus. For both graphs model results are marked with red lines and data with black. The associated pressure stimulus is depicted in Figure 2A. For this stimulus we analyzed five models. We first describe results obtained with the three linear models, analyzing the impact of including one, two, or three Voigt bodies. Second we discuss results obtained with the nonlinear models analyzing the impact of including more advanced description of the wall strain. For this stimulus we did not analyze the integrate-and-fire model, since we did not anticipate any added effect of this model because of the input rage of the pressure stimulus. The three linear models include a component determining the wall strain, described using a linear elastic function of pressure, a component representing mechanoreceptor stimulation, described using one, two, and three Voigt bodies, and a component predicting the BR firing rate. The three models have, , and parameters, respectively, as well as two additional parameters and associated with the sinusoidal stimulus. In [20, p. 695] the authors indicated that phase measurements are less accurate than amplitude measurements due to the inaccuracies associated with assigning interspike intervals to bins. Thus, the parameters and were added to the parameter set. Sensitivity analysis together with subset selection allowed us to identify four uncorrelated parameters including, and, which were estimated for all three models. The nominal values for the model parameters (listed in Table 4) were computed as follows. The parameter (), where is Young' s modulus (mmHg), (mm) is the wall thickness, and is the zero pressure radius as described in the Methods section (see also Table 2). We use approximating a lower bound to values observed in a previous study [75]. In [76, Figure 1] Feng et al. provide detailed measurements of the external diameter and thickness for the rat aortic arch, measured in adult male Sprague-Dawley rats. They found that in the region with aortic BR endings the average values of and. Using these values we compute (). We note this parameter and were highly correlated indicating equivalent fits could be achieved through adjustment of either parameter. No direct experiments exist allowing estimation of nominal values for the elastic and viscous constants associated with mechanoreceptor strain. These parameters appear only in the dynamic part of the model and determine the adaptation time-scales. To ensure that the three models are distinct, it is essential that parameters representing time-scales are separated, otherwise the models would essentially reduce to one. This knowledge, along with values chosen in the study by Bugenhagen et al. [56] motivated our choice for nominal parameter values. To avoid the problem of structural nonidentifiability [65] we rescaled the parameters as follows and for. The full list of the model parameters together with their initial conditions, units and literature reference is provided in Table 2. As for the stimulus, the average pressure (127 mmHg) and the amplitude (5 mmHg) was provided in [20]. To compute the frequency and the shift of the pressure, we digitized the stimulus provided in [20], Figure 2A, and then fitted to a sinusoidal function, obtaining and. As noted in Figure 6, results of parameter estimation with each of the three models were indistinguishable, though estimated parameter values varied significantly, the latter is due to added complexity associated with adding more Voigt bodies. The fact that graphs were almost identical was also reflected by the least squares cost RMSE (and the coefficient of determination) for models, and we obtained 2. 522 (0. 949), 2. 507 (0. 950), and 2. 495 (0. 951), respectively, see Table 4. Next, we investigated the impact of including more complex wall models. Additionally, we incorporated a nonlinear response wall model, and a viscous wall model. To be more precise we compare the BR response of the following three models, and described using 7,8, 9 parameters plus the two parameters associated with the stimulus. We examined the ability of each of these models to fit the sinusoidal stimulus. Sensitivity analysis and subset selection allowed us to estimate 4–6 parameters. All models allowed us to estimate, , , and. In addition, for the nonlinear elastic model was added to the subset and for the viscoelastic model and were added to the subset. Given that the more complex nonlinear models allows estimation of more parameters, one should anticipate better results. But due to the limited dynamics embedded within the pressure stimulus, adding more complex wall models did not improve results as reflected by the least squares cost RMSE (and the coefficient of determination), which for, and gave 2. 507 (0. 950), 2. 517 (0. 950), and 2. 458 (0. 952), respectively; see Table 5. Step-increase stimulus: This section presents results with the same five models previously used for prediction of the BR response with the sinusoidal pressure stimulus. As with the sinusoidal stimulus we do not test the integrate-and-fire model, due to the nature of the input stimulus. Again, we first discuss results obtained with the three linear models, and followed by results obtained using the more complex nonlinear and viscoelastic wall models. Studies were done to capture the effect of overshoot and adaptation in response to four input stimuli varying in the magnitude of the pressure step. All stimuli start at the same baseline pressure, and the step-increase was imposed at the same time. As before the three models have, , and parameters, respectively, but functions describing the “smooth” step pressure increase (2) only involve one additional parameter, representing the onset of the step-increase. This parameter was not provided in [20]. Subset selection together with efforts to make model comparison possible resulted in. As reported in [20, Figure 5] the baseline pressure associated with the step-increase stimulus was set to 115 mmHg, and the step-increases (from the baseline) to 128,134,137,143 mmHg, respectively. Figure 7 (A–D) shows the ability of the three linear BR models to reflect observed overshoot and adaptation. Each panel shows the optimized firing rate. The least squares cost RMSE (and the coefficient of determination) of model for the optimized values of its parameters with respect to the four step-increases 128,134,137, and 143 mmHg were: 1. 860 (0. 899), 2. 677 (0. 919), 2. 420 (0. 969), and 1. 832 (0. 983). Marginal improvements were obtained with, which gave: 1. 800 (0. 905), 2. 702 (0. 917), 2. 390 (0. 970), and 1. 823 (0. 983), and finally, for the values were: 1. 764 (0. 909), 2. 700 (0. 918), 2. 390 (0. 970), and 1. 809 (0. 983), see Table 4. Similar to the sigmoidal stimulus, no improvements (results not shown) were obtained with the more advanced nonlinear and viscoelsatic wall models. Square stimulus: The square stimulus is characterized by a constant pressure input followed by a step-increase after which the pressure is decreased to its baseline value. This type of stimulus primarily tested the models' ability to reflect PED followed by recovery, although other features including adaptation and overshoot are also shown. Similar to previous studies we first investigated the simpler linear models including one, two and three Voigt bodies. For the square input stimulus, in Figure 8A, we plot BR firing rate data extracted from Saum et al. [32] (circles) and the corresponding optimized fit using (solid line), changing the number of Voigt bodies did not improve the model response. This model has 7 parameters and additional two and related with the input stimulus (3). Subset selection together with our effort to make model comparisons possible made us estimate the parameters. The least squares cost RMSE (and the coefficient) with optimized parameters was 7. 384 (0. 862 for), see Table 6. While the model, as anticipated, was able to produce overshoot and adaptation, this model was not able to capture PED accurately. We hypothesize that the inability to show PED is due to the simple linear firing rate model, which does not allow the BR firing rate to cease for sub-threshold stimuli. Thus, we first investigated the impact of exchanging the linear BR firing rate model with the integrate-and-fire model. Including the integrate-and-fire model clearly improved results (not shown) though with the linear wall model it was difficult to accurately fit the data both during adaptation and recovery. Subsequently, we analyzed the impact of exchanging the linear wall model with the nonlinear wall model, keeping the integrate-and-fire model. Results with this model () is shown in Figure 8B. This figure shows the recorded BR firing rate (circles) and the model fit (solid line) in response to the square pulse stimulus. Model parameters estimated include. Optimized parameter values and units are given in Table 6 together with the and RMSE errors. Finally, we investigated the impact of adding a viscoelastic wall model, which did not provide any additional improvements. Simultaneous fits: Figure 7 showed that linear models can exhibit overshoot, adaptation, and can fit the firing rate data for all four step-increases, though as reported in Table 4, each step-increase resulted in significantly different parameter estimates. However, data are extracted from experiments done within the same fiber, thus we expected only small variation in parameter values. We performed additional optimizations to investigate if the observed differences in the parameter estimates, were simply a result of performing optimizations for one stimulus at the time. To remedy this problem, we estimated one set of parameters for all four step-increases. Results of this simulation are shown in Figure 9A (computed with the model). This simulation confirms that the simple linear model cannot estimate one set of parameters that allows simultaneous fit of the response to all four pressure stimuli. Similar results were obtained with the other models. In particular, it should be noted that the overshoot is diminished for the smaller step-increases, and that the model was unable to capture the correct baseline firing rate. In contrast, when including a nonlinear elastic wall we were able to estimate one set of parameters that allowed us to simultaneous fit the response to all four pressure stimuli. This model accurately reproduced the baseline firing rate as well as the overshoot and adaptation observed in response to the step-increase (Figure 9B). We hypothesize that this difference is due to larger range of pressure within the applied stimuli, where the known nonlinear behavior of the arterial wall deformation plays an important role. It is known that arteries appear stiffer at higher pressures than at lower pressure. Thus the nonlinear wall model significantly improves the fit. In the previous section we showed the ability of our proposed linear and nonlinear BR models to fit the firing rate data measured from rats. It is well known (see section Methods) that the BR firing rate can exhibit a number of qualitative characteristics including saturation, threshold, adaptation, overshoot, PED and rectification. The quantitative data used to test the model in the previous section showed adaptation, overshoot, and PED, in response to a sinusoidal (with fixed amplitude) and step changes (increase/decrease) in blood pressure. However, these stimuli did not test saturation, threshold, or rectification. Although the models show adaptation, no clear conclusion could be drawn to determine how many Voigt elements (time-scales) were needed to reflect known BR firing rate observations. Now we show our preferred model with estimated parameters, including nonlinear deformation of the elastic wall, two Voigt bodies for computing nerve ending stimulation, and a leaky integrate-and-fire model for predicting firing rate, exhibits the features not yet studied experimentally. This was done using ramp and sinusoidal (with varied amplitude and frequency) pressure stimuli. Rectification: Figure 10A presents the model' s response to a sinusoidal wave pressure stimulus with various amplitude. This simulation is motivated by the observation of Brown et al. [20, Figure 6] that a 2. 5 increase in amplitude of the sinusoidal stimulus resulted in an increased amplitude of the firing rate, with a lower mean firing rate. Moreover, it was noted that for large amplitude stimulation the firing rate ceases during the trough of the pressure wave. These two observations are referred to as rectification. One could question if the simpler linear model is able to display this phenomena. The linear wall model would certainly be able to reproduce the increased amplitude for a single stimulus, but again, if multiple stimuli were tested, correct predictions require the nonlinear wall model. Moreover, the ability of the firing rate to cease requires the threshold built into the integrate-and-fire model. With the simple linear neuron model, the firing rate would become negative, which does not represent what happens physiologically. Threshold and saturation: Two other prominent firing characteristics are threshold and saturation. In [27, Figure 5] Seagard et al. noted that BRs with a higher threshold pressure were less sensitive, had lower discharge rates, and had higher values for saturation. Receptors with higher discharge rates were also more sensitive and were found to have afferent fibers with greater conduction velocities. In Figure 10 B we show that our model is able to reproduce qualitatively similar saturation features. Adaptation: Even though our quantitative models were able to capture adaptation, it was noted that results with one, two, or three Voigt bodies were similar, in other words, the models could not clearly distinguish if the adaptation process included one or three time-scales. Yet, several authors (e. g. , [11], [12], [30], [77]) have hypothesized that adaptation occurs with more then one time constant. It is also known that the muscle spindle can produce a response of this kind to a clipped-off ramp stretch [78]. Figure 10 C shows that the studied model admits the fast adaptation and the slow adaptation in agreement with experiments. We also plot an exponential fit and show that a similar adaptation is not possible by only one exponential function. This qualitative feature made us include two Voigt bodies in our preferred model, a conclusion that could not have been made strictly from quantitative simulations presented in the previous section. Asymmetry: In Figure 10 D we show that our preferred model clearly exhibits asymmetry when exposed to a ramp-up followed by a ramp-down pressure stimulus, which agrees with experiments (see e. g. , [23]).
The objective of this study was to develop a mathematical framework for constructing computationally efficient and accurate BR models, which in contrast to the existent models, are able to reflect all known qualitative BR firing features as well as fit quantitative data. Our overall aim was not to focus on a concrete experimental species but rather to formulate a family of BR models, which could potentially be included in a more comprehensive model of CV system. Quantitative computations were done comparing our models to experimental measurements by Brown et al. [20] and Saum et al. [32]; while qualitative studies were performed to show that our preferred generic model is able to exhibit all known firing rate responses. All models used blood pressure as an input and computed the BR firing rate as an output. Although our procedure was designed to be generically applicable to various species and multiple types of baroreceptors, we tested our models using only quantitative data from experiments preformed using aortic baroreceptors from rats. We believe that this is the first work that offers a systematic approach to building and evaluating BR models with the objective to provide the simplest possible family of generic models. Our modeling framework first analyzed the known physiology and common features of the firing rate observed in the BRs of various species. Second we generated submodels describing each stage of the physiological response: arterial wall deformation, stimulation of mechanosensitive channels found in the BR nerve endings, and generation of action potentials. Finally we modeled the BR system by combining the submodels in various configurations (summarized in Table 3). Each of these configurations was tested in order to determine the contributions of each component to the transduction of the BR signal. This process allowed identification of the importance of nonlinear effects of two critical sub-systems in the BR response, the arterial wall and the neuron itself. This framework advanced the state of BR modeling by first evaluating models comparatively with respect to the same data and features, second by generating a model which fits all known characteristics of BR firing qualitatively, and third by developing a model which is capable of fitting multiple data sets of BR firing rates quantitatively. A particular insight was revealed by consideration of BR models with various descriptions of the arterial wall. Applying our framework demonstrated the insufficiency of linear wall models' representations of the response of a single BR neuron to multiple step-pressure inputs (see Figure 9A). A nonlinear elastic wall model was required to implement a model capable of accurately fitting the BR response to multiple pressure levels with one set of parameter values (see Figure 9B). The choice of this model is further motivated by the well known fact that arteries exhibit nonlinear deformation with saturation at both high and low pressures [23], [27]. Additionally by applying our framework and considering the effects of including the viscoelastic wall model, we found that the additional complexity did not contribute to better definition of BR dynamics, despite previous studies having shown wall deformation does have viscous components [45], [79]. This is likely due to our modeling choice for nerve ending stimulation. This portion was modeled using two Voigt bodies in series to allow adaptation at multiple time-scales. Data is not available to separate the viscoelastic part of the wall-deformation with the viscoelastic deformation associated with stimulation of the mechanosensitive channels, thus indirectly our model exhibits both features. One explanation would consider the first Voigt body to be associated with wall deformation while the second is associated with nerve ending deformation. Moreover, it should be emphasized qualitative simulations were needed to show that the two Voigt bodies allow multiple time-scales, a feature we were not able to extract from simulations alone. These considerations, and our studies, affirm the importance of viscoelastic effects; however, in terms of simplicity it is advantageous to isolate the viscoelastic components within the model, and further we note linear viscoelastic effects are sufficient to capture the dynamics of BR firing when coupled with a nonlinear elastic total deformation of the arterial wall. To our knowledge, this study provides the first direct measure of the importance of incorporating various time-scales in BR models. It is believed that various time-scales in the adaptation process are due to the viscoelastic coupling of the nerve ending to the arterial wall. We chose to emphasize this in our modeling process by considering different numbers of Voigt bodies in series with a spring. In Table 4 we show the results of testing three models, , and differing only with respect to their nerve ending models, , and, respectively. Our findings indicate that no more then two timescales in the adaptation process are needed in order to achieve a very precise fit to the data. This conclusion is closely related to the fact that we tested our models using rat data with fairly limited pressure-stimulus response as only this type of experiments are currently available. To test this component more carefully, it is essential to analyze data recorded over longer time-scales. Another insight afforded by this investigation highlights the importance of nonlinearities in the neural response to mechanoreceptor strain. As hypothesized previously [30], our study affirms the nonlinearities of action potential generation, even for the leaky integrate-and-fire model are sufficient to produce the hysteretic phenomenon of PED. In contrast the simple linear model of firing in response to mechanoreceptor strain does not allow for the asymmetric responses seen in PED as well as in the response to sinusoidal stimulus with high amplitude. The nonlinear-elastic wall in combination with two Voigt bodies modeling mechanoreceptor stimulation responds in an equal but opposite manner to rising and falling pressure, thus the change in firing rate with the linear model is symmetric to step-increase and step-decrease, which is not reflective of the data. We affirm the hypothesis that the neuron itself is responsible for generating PED, as this feature was robustly represented by the leaky integrate-and-fire model regardless of the mathematical description for arterial wall strain. This would provide a good explanation for the observation of PED in multiple species, many of which have a high degree of variability in the viscoelasticity in their respective arterial walls. The results and insights generated through application of our proposed modeling framework are not limited to those presented in this study. In addition it provides a means to identify which features and what level of detail of the underlying physiological systems are of greatest significance in generating BR dynamics. This ability is useful in developing experiments which may be able to isolate physiology responsible for a given phenomenon, such as the responsibility of the neuron in generating PED. Further this approach provides evaluative power to make design decisions when developing a model for a specific data interpretation or simulation task. An example of this follows from our insights into the role of the arterial wall in BR signal transduction. Although the arterial wall may best be modeled using viscoelastic theory, our framework allows a modeling decision to be made in favor of simplicity if only the output dynamics are of interested. This investigation further suggests a methodology for integrating a model generated in this manner into a model of larger scope. Suppose a mathematical representation of an overall baroreflex system (see Figure 1) is desired to reflect only normal physiological conditions, then it may be sufficient to use only simplified description of the BR signal. For example a simple linear firing rate model may be adequate for simulations operating in the range above the firing rate threshold. However, to reflect heart rate at various abnormal physiological conditions a more complex model combining nonlinear deformation with the leaky integrate-and-fire model may be necessary. Additionally, application of our modeling approach to a larger CV model might reveal features of the BR subsystem with importance in maintaining homeostasis. We hypothesize that overshoot, adaptation and recovery, features of the BR firing in response the extremes of pressure waves, are critical for regulation of blood pressure during stressful situations, such as a head-up-tilt experiment.
Many people have experienced lightheadedness when standing up, yet the exact cause of this phenomenon remains unknown. For some people, lightheadedness occurs because of anomalies in the blood pressure control system, which keeps blood flow and pressure at homeostasis. One way to explore this system is via mathematical modeling, which can offer valuable insights into the complex dynamic processes. This study develops a framework for modeling activity of the baroreceptor neurons. The models consist of three components reflecting three physiological mechanisms relating blood pressure to the baroreceptor firing rate: modulation of arterial blood pressure causes dilation of the arterial wall, stimulating mechanoreceptors within the baroreceptor nerve endings, emanating from the aortic arch and carotid sinus, which in turn modulates the firing rate of the baroreceptor neurons. This signal is integrated in the brain stem, stimulating baroreflex efferents to counteract the pressure increase. In this study, we review the main static and dynamic features of the baroreceptor firing activity, and show, using a combination of modeling techniques and rat aortic baroreceptor data, how to build a computationally efficient, yet biologically correct model. These models are important components for describing efferent responses, such as: heart rate, contractility or stroke volume.
Modeling the Afferent Dynamics of the Baroreflex Control System