Hysteretic behavior of bladder afferent neurons in response to changes in bladder pressure

Background Mechanosensitive afferents innervating the bladder increase their firing rate as the bladder fills and pressure rises. However, the relationship between afferent firing rates and intravesical pressure is not a simple linear one. Firing rate responses to pressure can differ depending on prior activity, demonstrating hysteresis in the system. Though this hysteresis has been commented on in published literature, it has not been quantified. Results Sixty-six bladder afferents recorded from sacral dorsal root ganglia in five alpha-chloralose anesthetized felines were identified based on their characteristic responses to pressure (correlation coefficient ≥ 0.2) during saline infusion (2 ml/min). For saline infusion trials, we calculated a maximum hysteresis ratio between the firing rate difference at each pressure and the overall firing rate range (or Hmax) of 0.86 ± 0.09 (mean ± standard deviation) and mean hysteresis ratio (or Hmean) of 0.52 ± 0.13 (n = 46 afferents). For isovolumetric trials in two experiments (n = 33 afferents) Hmax was 0.72 ± 0.14 and Hmean was 0.40 ± 0.14. Conclusions A comprehensive state model that integrates these hysteresis parameters to determine the bladder state may improve upon existing neuroprostheses for bladder control.

threshold is reached and gradually increase their firing as the bladder fills and pressure rises [6,10,12,16,17]. Pressure thresholds for these neurons are typically between 5 and 20 cmH 2 0 [3,10,17,18], with neurons of different pressure thresholds being successively recruited as the bladder fills [16]. Some of these fibers decrease their firing or plateau at high pressures [14]. However, the relationship between bladder afferent activity and intravesical pressure is not a simple linear one. Firing rate responses to pressure can differ depending on prior activity, demonstrating hysteresis in the system [8,16].
Hysteresis is a nonlinear phenomenon in which the output of a system depends on both the current input and recent history. It commonly occurs in ferromagnetic and ferroelectric material but is present in other systems including many mechanosenory systems [19], for example, in stretch receptors in cray fish [20] and muscle spindles in cats [21]. Hysteresis is also a property of smooth muscle, which makes up the detrusor layer in the bladder wall [22].
Hysteresis in the bladder pressure-bladder afferent relationship has been highlighted in previous LUT literature [8,16,23] but has not been quantified or compared across different bladder states. The goal of this study was to quantify the hysteretic relationship between bladder afferent activity and bladder pressure during non-micturition bladder contractions as the bladder is being filled and also when it is at an isovolumetric state. Neural activity was recorded from sacral dorsal root ganglia at the S1 and S2 levels, which contain the cell bodies of pelvic neurons, in cats. Neurons that responded to bladder filling demonstrated a quantifiable hysteresis that was similar to examples found in the published literature. A better understanding of this hysteretic relationship could be utilized to implement a comprehensive state model for a closed-loop bladder neuroprosthesis, in which pressure is estimated from afferent activity and stimulation is delivered for bladder control accordingly. Such a system could be used clinically for patients suffering from spinal cord injury and other neurogenic bladder disorders.

Subjects
Five intact adult male cats (age: 0.9-1.4 years old, 4.2-5.2 kg, domestic short-haired, Liberty Research, Inc, Waverly, NY) were used in non-survival experiments in this study with one cat used per experiment. All procedures were approved by the University of Michigan Institutional Animal Care and Use Committee, in accordance with the National Institute of Health's guidelines for the care and use of laboratory animals.

Experimental setup and surgical procedure
Animals were initially anesthetized with a ketamine (6.6 mg/kg)-butorphanol (0.66 mg/kg)-dexmedetomidine (0.033 mg/kg) intramuscular (IM) dose, intubated, and then maintained on isoflurane anesthesia (2-4 %) during surgical procedures. Respiratory rate, heart rate, end-tidal CO 2 , O 2 perfusion, temperature, and intraarterial blood pressure were monitored continuously using a Surgivet vitals monitor (Smiths Medical, Dublin, OH). Intravenous (IV) lines were inserted into one or both cephalic veins for drug and fluid infusions. Intravenous fluids (1:1 ratio of lactated Ringers solution and 5 % dextrose) were infused at a rate of 5-10 ml/kg/h and increased up to 30 ml/kg/h during surgery as needed.
A catheter was inserted into the bladder for intravesical fluid infusion and pressure monitoring. In experiments 1, 3, and 4, an abdominal midline incision was performed to expose the bladder and a 6 Fr supra-pubic dual-lumen catheter (Laborie, Williston, VT) was inserted into the dome of the bladder and secured with a purse-string suture. A single-lumen 3.5 Fr catheter (Utah Medical Products, Midvale, UT) and dual-lumen 3.5 Fr catheter was inserted into the bladder via the urethra in experiments 2 and 5, respectively.
Following bladder line placement, a midline dorsal cut was made to expose vertebrae from L7 to S3. The spinal muscles were reflected from the vertebrae and a laminectomy was performed to access sacral DRG (S1-S2) in the cat [18]. Iridium oxide microelectrode arrays (5 × 10 and 4 × 10 ICS-96 MultiPort split planar arrays, Blackrock Microsystems, Salt Lake City, UT) were implanted into the DRG either bilaterally or unilaterally using a pneumatic inserter (Blackrock Microsystems). These types of multielectrode arrays are a standard approach for identifying and recording from dozens of neurons simultaneously in DRG [18,[23][24][25][26][27] and peripheral nerves [28][29][30]. For experiments 1 and 5 the 5 × 10 array was inserted in the left S1 DRG and the 4 × 10 into the left S2 DRG. For experiments 2, 3, and 4, 5 × 10 arrays were inserted bilaterally in S1 and 4 × 10 arrays were inserted bilaterally in S2. Microelectrode shank lengths were either 0.5 or 1.0 mm with 0.4 mm inter-shank spacing. The ground wires were connected to a stainless steel needle inserted in the skin lateral and caudal to the laminectomy incision site and reference wires were placed near the spinal cord. Animals were then transitioned to alpha-choloralose (70 mg/kg induction; 20 mg/kg maintenance; doses given IV) for subsequent testing. Alpha-chloralose anesthesia was augmented with buprenorphine (0.01 mg/kg; given every 8-12 h IV).
DRG neural data was acquired at a rate of 30 kHz and band-passed filtered (250 Hz-7.5 kHz) using a Grapevine neural interface processer and Trellis recording system (Ripple, Salt Lake City, UT). A global amplitude threshold, between −20 and −35 µV (depending on the noise amplitude), was set for all electrode channels. Any signal crossing of the threshold was captured as a spike snippet and stored for offline analysis. Bladder pressure was monitored using a pressure transducer (DPT-100, Utah Medical Products, Midvale, UT) and transducer amplifier (World Precision Instruments, Sarasota, FL). The bladder pressure signal was recorded with the Grapevine system at 1 kHz. During testing (see "Experimental procedures" section) saline was infused into the bladder using a syringe pump (New Era Pump Systems, Inc., Farmingdale, NY). Figure 1 shows the experimental set-up.

Experimental procedures Slow fill trials
The bladder was first emptied using the bladder catheter. Sacral DRG neural activity was recorded while saline was infused into the bladder at a near physiological rate of 2 ml/min [31]. In most trials, this was done until dripping from the external meatus or around the urethral catheter (for those experiments where the catheter was inserted via the urethra) was observed. This infused volume was defined as the "leak volume" for a given experiment. For trials in which infusion was stopped before leaking was observed, then the leak volume from a prior fill sequence was assumed. For experiments 1-4, room-temperature saline (22 °C) was used; whereas for experiment 5, bodytemperature saline (41 °C) was used. Two infusion trials per experiment (cat) in which there were only non-voiding bladder contractions were used in the analysis.

Isovolumetric trials
In experiments 3 and 4, isovolumetric trials were performed with the bladder volume within 20-50 ml, while assuming negligible urine generation. Neural activity and bladder pressure were recorded for non-voiding bladder contractions.
After completion of all testing, animals were euthanized with a 3 ml intravenous dose of sodium pentobarbital (390 mg/ml) while under deep anesthesia.

Data analysis
After data collection, spike snippets were sorted in Offline Sorter v3.3.5 (Plexon, Dallas, TX), using principal component analysis followed by manual review to identify unique spike clusters. In MATLAB (Mathworks, Natick, MA), instantaneous firing rates for each Arrays were implanted in S1 and S2 DRG. Saline was infused into bladder either via a supra-pubic line or an intraurethal line. Intravesical pressure was monitored with a pressure transducer and amplifier. Both neural data and pressure were recorded with a Grapevine data acquisition system. Image modified from Bruns et al. [57] unit were calculated at intervals of 0.5 s. Discrete spike events were converted into a smoothed time series of firing rates using a non-causal linear filter with triangular kernel of width 1 s [27]. Bladder pressure was filtered (4 Hz low pass). Units whose firing rates highly correlated with bladder pressure (correlation coefficient, ρ ≥ 0.2) over the course of a saline infusion trial were identified as bladder units. These units were confirmed visually to increase firing with increasing bladder pressure as shown in the literature [6,10,16,18]. Hysteresis in the bladder pressure-firing rate relationship was calculated during bladder contractions using a method derived from Kosmulski et al. [32] for electrochemical capacitors. Sets of three contractions (pressure change ≥ 10 cmH 2 O, stratified by 25 % intervals of the leak volume for slow fill trials) were used for hysteresis calculations (Fig. 2a). The start and end of each contraction was determined based on pressure inflection points. For each contraction, the pressure trace was divided into 2 cmH 2 0 bins and the mean firing rate corresponding to pressure within each bin was calculated. The mean firing rate and pressure were plotted against each other (Fig. 2c). Pressure ranged from a minimum value (P min ) to a maximum value (P max ) with a corresponding firing rate range from FR min to FR max . Figure 3 shows a stylized diagram of pressure plotted against firing rate demonstrating how the hysteresis values are calculated. For each binned pressure value (P 1 , P 2 , P 3 , and P 4 ), the difference in firing rate (ΔFR) is calculated and divided by the firing rate range (FR max − FR min ). This ratio is defined as ΔFR rel . The following two hysteresis indices were computed: Hmax and Hmean. Hmax is the maximum ΔFR rel and Hmean is the average ΔFR rel . The start and end points were excluded in the Hmean calculation to avoid overrepresentation of the narrow ends of the pressure-firing rate curves. Hmax and Hmean were then averaged across the 3 contractions. Hmax and Hmean are dimensionless values ranging from 0 to 1, where 0 represents no hysteresis and 1 represents maximum hysteresis.

Statistical analysis
A linear mixed models ANOVA with volume range as the fixed effect and a random intercept of animals (experiments) was carried out to determine if there was any statistical difference in Hmax or Hmean values between volume ranges for both slow fill and isovolumetric trials. Given the small number of animals and variation in the number of bladder units observed per experiment, we went with the linear mixed models with experiment as the random term in the model. P values (significance level α < 0.05) were obtained by likelihood ratio tests of the full model with the effect in question (volume) against the model without the effect in question [33]. To test for significance between experiments, an ANOVA linear model was performed. Post hoc analysis was done using the Tukey-Kramer multiple comparison test. Relationships between Hmax and Hmean, Hmax and ρ, and Hmean and ρ were assessed using regression analysis. MATLAB and Microsoft Excel (Redmond, WA) were used to perform statistical tests. Values are reported as mean ± standard deviation.

Bladder units
Seventy units from five experiments were identified as bladder units (ρ ≥ 0.2), with 57 units in S1 DRG and 13 units in S2 DRG. These units were quiescent at very low pressures and increased firing as bladder pressure increased. Four of those units did not correlate with bladder contractions and were excluded from further analysis. Forty-six units were present in the slow fill trials with average ρ = 0.57 ± 0.16 (correlation analysis done over the infusion period, n = 5 experiments). Correlation coefficients did not differ significantly across experi-  For slow fill trials, there was no difference in hysteresis observed at different volume ranges (Hmax: χ 2 (3) = 5.91, p = 0.12; Hmean: χ 2 (1) = 0.83, p = 0.36) (Fig. 5). Note, given the variability in the number of contractions produced for the different animals, some units only had hysteresis values for one given volume range, whereas others had hysteresis values for multiple volume ranges. Figure 6 shows the breakdown of Hmax and Hmean across different volumes for the isovolumetric trials. There was no difference between volumes (Hmax: χ 2 (1) = 0.70, p = 0.40; Hmean: χ 2 (1) = 3.04, p = 0.08).
The effect of different firing rate calculation approaches was analyzed. In addition to the non-causal linear filter with triangle kernel (described in "Methods"), two other  ii Firing rate for a bladder unit that correlates well with pressure during the same time period (symbols correspond to those used in i, with correlation coefficient (ρ) indicated). iii Firing rate for a bladder unit that is weakly correlated with bladder pressure. b Waveform and interspike interval (isi) histogram, with bin width = 5 ms, for neurons in a. c Bladder pressure plotted against firing rate of respective neurons from a for each contraction during the designated periods at 25-50 % leak volume (left column) and 50-75 % leak volume (right column). The Hmax and Hmean values are given above the plots. Arrows represent direction of cycle methods were used to calculate firing rates: instantaneous firing rate and boxcar. At every 0.5 s, the instantaneous firing rate was determined by calculating the reciprocal of the interval between the two spikes surrounding the time point [34]. For the boxcar method (rate histogram) the spike rate at every 0.5 s was calculated as the spike count within 1 s bins divided by the bin size [11,14]. The correlation coefficient, Hmax, and Hmean were compared for all units for each firing rate approach for all slow fill trials. Overall, little difference was seen in the hysteresis values obtained using the different firing rate methods. Figure 8 shows an example from one experiment comparing the hysteresis values and correlation coefficients for 3 units. Across all experiments, the triangle kernel method resulted in the highest ρ between firing rate and pressure for a given unit and was the method selected for our analysis.

Comparison to published literature
The hysteresis indices, Hmax and Hmean were calculated for bladder pressure-firing rate relationships from published literature (Table 1). In general, even though the data in these papers were collected under a variety of conditions and from different types of bladder neural signals and sources, similar Hmax and Hmean values were obtained, with Hmax = 0.64 ± 0.13 and Hmean = 0.39 ± 0.14.

Discussion
The goal of this study was to quantify the hysteresis observed in bladder pressure-bladder afferent relationships. A total of 66 bladder afferents across S1 and S2 DRG in 5 cats were identified that correlated with bladder pressure and bladder contractions. Hysteresis was observed in the firing rate responses of these afferents to bladder pressure. Analysis of contractions during slow fill and isovolumetric trials yielded Hmax values ranging from 0.38 to 1.00 and Hmean values ranging from 0.09 to 0.76 (Fig. 4). Hysteresis values did not differ across different volume ranges (Figs. 5, 6). In general, afferents that correlated better with bladder pressure tended to have lower Hmax values (Fig. 7b). However, there was no relationship observed between Hmean and how well an afferent correlated with bladder pressure (Fig. 7c). On closer examination, the inverse relationship seen between Hmax and correlation coefficient for the slow fill trials was observed in two of the five animals (experiments 2 and 3). For two other animals (experiments 1 and 5), this was actually a positive relationship. An inverse relationship between Hmax and correlation coefficient was also observed in both experiments (3 and 4) for the isovolumetric trials. Intuitively, it would seem to follow that there would be less hysteresis where there is a more linear firing rate-pressure relationship, but this was not the case for all animals. These individual differences in responses may partly be explained by how Hmax is calculated. Hmax is determined by only the maximum difference in firing rate, compared to Hmean which gives the average difference in firing rate over the whole loop. Hence, Hmax is more susceptible to outlier firing rate differences over a cycle as it simply gives the largest relative firing rate difference. Hmean, on the other hand, is the mean relative firing rate difference over the pressure range and provides a more representative measure of the hysteresis over all the cycles and seems to be less affected by large variability in firing rate over a cycle and from one cycle to the next (Fig. 2c). Hmean in general is lower if there is more variability in the amount of hysteresis across individual loops and if there is less hysteresis. Furthermore, if both Hmax and Hmean are large and the ratio of Hmax/Hmean is close to 1, this in general reflects more hysteresis in the system and larger, more uniform loops. Higher ratios of Hmax/ Hmean indicate more peculiar loops [32] and more variability in the loops between cycles.
Several factors may have contributed to variability in hysteresis loops (e.g. Fig. 2c ii) beyond natural irregular action potential firing, including inconsistencies in spike isolation and sorting and the firing rate calculation. We sought to minimize these sources of variability through several approaches. We utilized consistent unit thresholding and a fixed cutoff for determining bladder units (correlation coefficient ≥ 0.2). All sorted bladder units were manually verified by the same person. We also performed a preliminary hysteresis analysis comparing different firing rate calculations. We did not see a difference in computed hysteresis values between different methods (Fig. 8), and proceeded with a consistent firing rate calculation. Nevertheless, we cannot rule out potential contributions of these and other sources of error in the calculated hysteresis values. Still bladder afferent units had higher Hmax and Hmean values (on average across both data sets higher than 0.7 and 0.4, respectively), particularly when compared to non-bladder units which generally had a very small Hmean (close to 0). We applied the same equations for calculating Hmax and Hmean to neural activity-bladder pressure (or tension) plots in published literature ( Table 1). The range of Hmax and Hmean values from those plots were comparable to those values calculated from our data (Figs. 4,5,and 6). Note that the data in Table 1 was collected from a variety of animal models (cat, pig, and rat), fill rates, and signal sources (dorsal root ganglia and pelvic nerves). An advantage of the hysteresis calculation method used in this study is that it does not depend on the actual units of the input-output signals. Further testing could be done to rigorously compare effects of different types of signals and/or signal sources and different fill rates (for example, comparing rapid injection of fluid versus more physiological fill rates). Though our results were generally consistent, there was variability in the data sets. Some slow fill trials had many contractions (up to 16) whereas others had only 3 contractions. Thus we chose a set number of contractions per volume range. In addition, not all volume ranges within a fill cycle met our criteria for contraction count and size. Therefore, all volume ranges were not represented in all trials and animals. Some animals also had more bladder units than others. Another limitation of this study was the small sample size for the isovolumetric data (only 2 animals), though the results did overlap with the values from the slow fill trials. Another potential confounding effect was the difference in the infused saline temperature for experiment 5 compared to the other experiments. In experiment 5, saline was infused into the bladder at body temperature compared to room temperature for the other experiments. Room temperature saline infusion is a standard procedure in bladder studies [6,14]; however, functional differences can occur depending on the temperature of the saline [36]. We did see a higher average Hmean value for the "warmer" infusion compared to 3 out of the other 4 "cooler" infusions which may be due to this difference in temperature, though there was no difference observed for Hmax.
The origins of the pressure-firing rate hysteresis could be both myogenic and neurogenic. This may be due to the intrinsic mechanical properties of the bladder [37]. The bladder wall consists of smooth muscle, elastin, and collagen and has non-linear elastic, viscous, and plastic properties [38,39]. With repeated filling and emptying of the bladder, length-tension curves of the bladder produces characteristic hysteresis loops [40]. Even in the absence of neural innervation, outside of the body, the bladder muscle displays this hysteretic property, further demonstrating that hysteresis is inherent in the muscle. Subjecting bladder muscle strips to sequential stretching and relaxing resulted in shifts of the tension-time and length-tension curves in a manner typical for materials with hysteresis [41]. This non-linear behavior is also demonstrated in stress-strain curves after cyclic loading of the bladder muscle [38,42].
This hysteretic relationship may not only be explained by myogenic factors. Neuronal mechanisms may also play a role, though the particulars of this are not clear. Hysteresis could result from the membrane properties of the afferent fibers themselves and the mechanosensitive receptors. The primary mechanosensitive receptor in the bladder is thought to be the transient receptor potential vanilloid 4 (TRPV4) channel [43]. It is not known if the TRPV4 receptors demonstrates hysteresis, though another family member, the TRPV3, a temperaturesensitive receptor [44], has been shown to demonstrate hysteresis. Furthermore, other mechanosensory systems such as muscle spindles and joint receptors in cats have been reported to demonstrate hysteresis [21,45]. Many mechanoreceptors have also been reported to show adaptation with decreases in frequency and/or amplitude with increasing input stimuli which may be attributed to the viscoelastic nature of the receptors [46] and some bladder afferent responses may also demonstrate adaptation [47]. However, more detailed studies would be required to confirm neurogenic causes of this non-linear behavior or if it was solely due to the biomechanics of the muscle itself, which is beyond the scope of this study.
The rate of bladder filling has been reported to affect the bladder's intrinsic mechanical properties, with plastic behavior dominating during natural fill rates and viscous forces coming into play during rapid filling. Bladder compliance also decreases with increasing fill rates [48]. We infused saline at a near-physiological rate, in the medium-high range [31], so both plastic and viscous forces likely contributed to our results. Additionally afferent activity has been reported to be lower for a given pressure at higher fill rates (compared to more natural fill rates) [9]. Even though this may have occurred, the hysteretic relationship is still evident and the amount of hysteresis we see is comparable to values calculated from previous studies. Those studies used different fill rates and one of them, specifically looked at tension versus afferent activity (Table 1), whose relationship is reported to be independent of fill rates [9]. Ultimately, our goal is to develop a real-time closed-loop neuroprosthesis to determine bladder pressure from neural activity; hence, we are focusing on relationship between bladder pressure and bladder afferent activity.
Given the hysteretic relationship between bladder afferent activity and bladder pressure, a simple linear equation does not adequately describe the relationship between these two. There are closed-loop bladder control approaches that do not take into account hysteresis, which have demonstrated minimal success so far. For example, some studies have focused on estimated bladder volume instead of pressure [49][50][51], though pressure may be a more reliable estimator than volume [52]. Another study has shown a proof of concept closed-loop model using direct pressure monitoring as the feedback signal [53], but this requires a reliable way to directly measure bladder pressure that comes with its own challenges [54]. A closed-loop model for bladder control that takes into account hysteresis may result in a more realistic estimation of bladder pressure from afferent activity than using a simple linear regression model [23]. One approach could be using piecewise polynomials to model tonic increase in bladder pressure and contractions, and to use estimated hysteresis to adjust for the changing pressure-firing rate relationship on the rising and falling edges of contractions. Another method could be incorporating derivatives to model the hysteresis loop, similar to models proposed for dynamic hysteresis loops [55]. As we did not see differences in hysteresis values across different volume ranges, one hysteretic value may be applicable at different bladder volumes. Further experiments and algorithm development are needed to determine an optimal approach, which will also need to account for the effects of movement that will occur in a behaving implant recipient.

Conclusions
We did a quantitative assessment of the hysteretic relationship between bladder pressure and afferent activity, which has been previously commented on but not described quantitatively. The amount of hysteresis is similar across slow fill and isovolumetric trials and at different volumes and is comparable to values we calculated from published examples. These results provide more information on the relationship between bladder pressure and afferent activity and could be utilized in a closed-loop model for a bladder neuroprosthesis.
Abbreviations DRG: dorsal root ganglion (or ganglia); LUT: lower urinary tract; S1: sacral level 1; S2: sacral level 2; FR: firing rate; P: pressure; Hmax: maximum hysteresis ratio between the firing rate difference at each pressure and the overall firing rate range; Hmean: mean hysteresis ratio between the firing rate difference at each pressure and the overall firing rate range.