Introduction

Micturition is regulated by a dynamic reflexive process that switches between periods of storage and release of urine held in the bladder. This cycle is mediated by a complex network of autonomic and somatic neurons in the central and peripheral nervous systems [13].

The tibial nerve has a significant influence on micturition [4, 5], because projections from the tibial nerve to the spine inhibit afferent sensory neurons [68] and potentially several midbrain nuclei, including the Periaquaductal Grey (PAG) and Pontine Micturition Centre (PMC) [911]. This inhibition has been shown to act through both classical and opioidergic mechanisms [12]. This differential mechanism of action likely explains the ability of the tibial nerve to inhibit both the nociceptive and non-nociceptive activity of the bladder [13].

The ability of the tibial nerve to down-regulate bladder activity has led to the formal adoption of tibial nerve stimulation (TNS) as a clinical intervention for the management of symptoms of the lower urinary tract [13, 14]. Non-invasive transcutaneous TNS (known as TTNS) is used as a clinically approved treatment for the management of incontinence and overactive bladder [15]. However, despite the widespread adoption of this technique in clinical practice, key aspects of the dose-related efficacy of TTNS remain unexplored. Specifically, preliminary rodent and feline evidence suggests that there may be a frequency-dependent effect of low-frequency TTNS intervention driving a paradoxical up-regulation of bladder activity [1618].

We aimed to explore whether the preliminary evidence for a paradoxical frequency dependence manifested as a functional effect in a human population and potentially shed light on a possible mechanistic explanation for the effect. We sought to achieve these goals with a two-pronged approach: conducting a pilot study in a healthy human population and developing a detailed computational model that can explain such frequency dependence analytically.

Therefore, we conducted a single-blind randomized control trial in a healthy adult population. In addition, we included a runoff period to study if any observed effects linger in nature. With our computational model of the lower urinary tract and its controlling neural circuits, we simulated the bladder control network.

It is hypothesized that by varying only the frequency of applied TTNS, it will be possible to selectively up- or down-regulate bladder behavior in healthy adult humans. Moreover, we hypothesize that modeling the system computationally will allow elucidating the specific mechanism behind the effect.

We found that a frequency-dependent up- or down-regulation of bladder activity could be induced in a healthy human population. Furthermore, our computational findings corroborated these results and indicated that the effect may be mediated by brainstem-specific tibial projections and supported by spinal inhibition.

We propose a mechanistic framework for this effect based on the idea of filtering afferent sensory activity. These results reveal for the first time the presence of a frequency-dependent effect of TTNS in humans and provide a considerable foundation for future clinical research pertaining to the treatment of several disorders of the lower urinary tract.

Methods

Ethical Approval

The study was approved by the local ethics committee of The University of Edinburgh (ID: 977504). Forty-eight people were recruited under three experimental conditions. Before participation, they gave their full informed consent in writing.

Study design

A graphical overview of the study can be seen in Fig. 1. In summary, the study was a single-blind design between subjects. Participants were pseudo-allocated to one of three groups such that the sample size across groups remained equally distributed (see below). Participants were not informed of this decision to ensure blinding was maintained.

Experimental Study Methodology

Shown is a high-level overview of the study used to validate the frequency-dependent effects of TTNS predicted by the computational model. Participants were healthy adults (18+) asked to abstain from any nicotine/caffeine for 12 hours, and any fluids for 2 hours before the study. Upon arrival they were asked to ingest 750 ml water after emptying their bladder. A 30 min digestion period was employed before stimulation. Participants were pseudo-randomly allocated between groups, and were blind to the condition. Created with BioRender.com/p01p196.

To ensure a relatively stable bladder baseline across the sample, several measures were taken. First, participants were asked to restrict their fluid intake for two hours and caffeine or nicotine for 12 hours before the beginning of the experiment -caffeine or nicotine affect urination [19, 20]. Second, participants were asked to empty their bladder as much as they were reasonably able upon arrival, before ingesting 750 ml of water.

Participants underwent a 30 minute "digestion period" during which they were able to undertake any activities they wanted as long as they remained seated. During this period, they were prepped for neurostimulation: Any leg hair was shaved from the ankle and lower leg using a sterile disposable razor, and the stimulation site was disinfected using a 70% isopropyl alcohol solution to remove contaminants.

To ensure that the objective outcome of the study could be recorded in all cases, after the digestion period, participants were told to inform researchers when they first felt the urge to urinate. Electrodes (Med-Fit 50mmx50mm Hydrogel Electrodes) were applied to the right leg 1 cm posterior to and 10 cm superior to the medial malleolus in accordance with established TTNS protocols [13]. Participants received stimulation (Digitimer DS7A High Voltage Constant Current Stimulator) until they felt the urge to urinate. The specific parameters of the stimulation differed according to group allocation (see fig. 1):

  • Group A: 1 Hz, 200 µs pulse width, motor threshold.

  • Group B: Control group, without stimulation, though they still underwent skin preparation.

  • Group C: 20 Hz, 200 µs pulse width, motor threshold.

Upon reporting the urge to urinate, the time-elapsed was recorded and the participants were asked to rate how intense they felt their urge according to a validated questionnaire [21] (see Supplementary Materials).

Runoff Period

To analyze the presence of lingering effects and to explore whether the hypothesized excitatory effect was reversible, an additional runoff period was included in the study methodology, as shown in Fig. 1). After reporting the urge to urinate, the participants underwent an optional 10 minute extension to the study where:

  • Group A received 20 Hz stimulation, 200 µs pulse width, motor threshold.

  • Groups B & C were asked to remain seated for 10 minutes and no stimulation was applied.

After this runoff period, a final urge intensity score was obtained before the end of the study.

Model Topology

The simulated network was composed of a set of interconnected neuronal units each comprising 100 individual neurons modeled in terms of conductance, or Poisson point processes that generated simulated action potentials (see Neuronal Model). The topology of the model was adapted from previous work that simulated reflex voiding [22]. The circuit remained broadly unmodified, with the exception of several key modifications, including simulation of the tibial nerve and its projections (see Fig. 2). We focus on normal bladder function; therefore, the nociceptive bladder afferents were not included in our model. To simulate the effects of opioidergic and classical inhibition, neurons and synapses were represented in terms of ionic conductance. The biophysical model used to represent the bladder was adapted from [23](see Bladder Model).

Overview of the Computational Bladder Control Model

A. Shown is a block diagram of the simulated neuronal circuit and bladder model. The model used a modified biophysical representation of the bladder produced previously (dashed box, [23]) which used the firing rates of the pelvic (blue), hypogastric (orange), pudendal (green) efferents to calculate bladder state. Tibial input used to modulate the circuit shown in purple. PMC: Pontine Micturition Centre, PAG: Periaquaductal Grey, PGN: Preganglionic Bladder Neurons, ASC: Ascending Interneurons, PUD: Pudendal Afferent, Onuf: Onuf’s Nucleus, Pel: Pelvic Efferent, Hyp: Hypogastric Efferent. Scissors represent possible severed projections. Labels 1 - 5 represent key regions of tibial modulation using opioidergic (1 - 4) or classical (5) inhibitory mechanisms. Created with BioRender.com/z61j143.

The tibial nerve was modeled to project to second-order sensory afferents (via classical inhibitory synapses) and all neuronal units within the PAG and PMC (via opioidergic inhibitory synapses) as shown in Fig. 2. To provide fine control over the level of tibial modulation applied to the system, the firing rate of the neuron and its projections were clamped and, if required, could be severed by disconnecting synaptic junctions. This ensured that while in a severed state, the neurons though still modeled had no impact on post-synaptic targets.

Bladder Model

The model used a biophysical representation of the bladder adapted from the work published by Lister et al. [23] modified to allow integration with the conductance-based representation of the micturition control circuit used in the present work. In summary, for each step of the simulation (t), the firing rates of three key efferents: Pelvic, Hypogastric, and Pudendal were used to calculate the internal activation constants - ωe/i/s -which governed detrusor excitation/inhibition, and sphincter activation respectively. From these values, the bladder state could be calculated and the internal detrusor pressure (PB(t)) used to update the afferent neuronal activity (the main input for the control circuit; see Fig. 2). To ensure computational efficiency, the sampling rate of the model was set at 50 Hz, with each step of the simulation (Δt) representing a 20-ms increment. A detailed description of the model can be found in [23].

Animal Model Surgery and Preparation

To fit the synaptic weights of the model (see Supplementary Materials), we utilized real bladder pressure and neural data from a previously published work [24] (see Fig. 3). During the constant infusion of the bladder, simultaneous recording of neural signals from the spinal cord was carried out. To measure bladder pressure and volume, a catheter was inserted into the bladder wall from the dome area and secured using a purse-string suture in such a way that one catheter head was connected to the pressure transducer and the second one to the infusion pump. To measure the voided volume, a digital scale was placed under the rat model and the residual volume was calculated by the difference between the infusion volume and the voided amount. A single-channel recording was performed using epoxylite-insulated tungsten electrodes with resistance of 300-500 kΩ and shank diameter of 75 µm (FHC Inc., Bowdin, ME USA). To facilitate three-dimensional control of the electrode tip, a micromanipulator (SM-15, Narishige Group Product, Japan) was used, providing the possibility of locating the electrode in the best anatomical position. To reduce the effect of noise on the recorded signals, a custom-made Faraday cage was employed to envelope the entire experimental setup.

Methodological approach used to obtain animal data

A. Experimental setup and recording site used to obtain bladder data and associated afferent neural signal. B. Bandpass-filtered raw neural signal (top trace), corresponding bladder pressure (middle trace), and smoothed continuous neural firing rate (1st order Butterworth filter with a 1 Hz cutoff frequency, bottom trace) obtained from the experimental setup. Created with BioRender.com/g49j208.

Neuronal Model

The model was built using the Brian2 package for Python [25]. The relationship between the pelvic afferent firing rate and bladder pressure was calculated according to an empirically derived relationship determined previously [26] where the firing rate of the primary bladder afferent at each moment (F(PB(t))) was defined as:

All neuronal units, with the exception of primary bladder afferents and the tibial nerve, were simulated using a Conductance-Based Adaptive Exponential Integrate-and-Fire (CAdEx) model [27] modified to include opioidergic inhibition. Through the inclusion of an additional input current (Iap) neurons could be made tonically active where required. This was particularly important in the simulation of PAG and PMC, where inhibitory feedback mechanisms reliant on tonic activity were present.

The tibial nerve and primary bladder afferents were simulated as Poisson point processes which generated simulated neuronal spikes at a pre-specified rate. Doing so allowed the firing rate of the pelvic afferent to be set according to the aforementioned mathematical relation and the firing rate of the tibial nerve to be clamped to a specific value (representing external stimulation of the nerve).

Each unit within the network was composed of 100 simulated neurons connected via random synapses (where no autapses were permitted) to introduce noise into the system.

Membrane potential (v) was defined as

where Cm was the membrane capacitance; gL the leakage current & EL its associated reversal potential; gA was the adaption current and EA its reversal potential; Iap the input to tonically firing neurons; and Isyn the input from synaptic connections. Also included were terms that govern the action potential (vth and Δt). Like the original model, gA was defined as:

where the rate of change of the current was determined by a rate constant (τA), the subthreshold adaption parameters (vA and ΔA), and the maximum subthreshold adaption conductance (A). The model also included a post-spike reset mechanism, which too remained unchanged.

Upon spiking, the membrane potential was reset to its resting state (vR), and the adaption current increased (δ gA). The purpose of the adaption current was to allow modeling of the spike-frequency adaption. Tonically active neurons did not display spike-frequency adaption to ensure that their activity remained consistent.

Tibial Nerve Stimulation Analysis

To computationally analyze the impact of tibial modulation, simulations were performed in a range of modulatory states. Tibial modulation was applied to a 500 second simulation run at 21 different frequencies (0-20 Hz, in 1 Hz increments, Nrepeats = 10 in each condition). Any aspect of simulated bladder function or circuit activity could be recorded during these simulations. Where neuronal firing rate activity was analyzed, instantaneous firing rates for each recorded time step were smoothed (1st-order Butterworth filter, with a 1 Hz cutoff). To analyze the effect of specific tibial projections, projections to regions of interest could be severed as necessary to prevent the formation of synaptic connections.

Results

Participants Display a Frequency-Dependent Sensory Response to TTNS

Before analysis was conducted, several datapoints had to be removed from the dataset due to:

  • Failure to abstain from caffeine/nicotine before participation (n = 1)

  • Ongoing treatment for relevant co-morbidities (n = 1)

  • Participants reporting the urge to urinate before the end of the digestion period, preventing administration of neuromodu-lation (n = 3)

  • Participant reported a greater than normal water intake under typical conditions (n = 1)

  • Reported urge duration was classified as an outlier (z-score > 2) (n = 1)

The removal of these datapoints reduced the statistical power of the analysis and led to a final sample distribution of NA = 12, NB = 14, NC = 15. The abnormality of this distribution required that, where possible, non-parametric analyses were employed. Despite this, the results of the study were nonetheless promising. Analyzing the effects of the intervention, it was clear that there was a differential effect of stimulation frequency on urge-onset time. A Kruskal–Wallis test indicated that there was a significant effect of treatment condition on the time to onset of the urge (H(2) = 18.660, p < 0.001). Non-parametric post-hoc testing (Mann-Whitney U) confirmed that 20 Hz TTNS stimulation significantly delayed the onset of the first urge compared to control (p < 0.01) (see Fig. 4A). In contrast, compared to control, 1 Hz TTNS stimulation was found to significantly decrease the onset of the first urge (p = 0.03, see fig. 4A). A power analysis reported that, despite the reduced sample size of the final analysis, the population was of sufficient size to achieve a power score of at least 0.80.

Experimental Study Results

A: effects of low (1 Hz), high (20 Hz), and placebo (0Hz) TTNS stimulation on the time elapsed to first sensation of urge. P-values obtaine from post-hoc Mann-Whitney U tests. B: effects of the runoff period on self-reported urge intensity. Error-bars = 95% Confidence Interval. Urge intensity self-reported on a scale of 0–4 (see methods), runoff period was 10 minutes in duration. C: bivariate and univariate kernel-density estimate plot displaying the relationship between the time elapsed before participants reported the urge to urinate, and the self-reported intensity of the urge at this stage (i.e., before the additional 10-minute runoff period).

TTNS Has Little Frequency-Dependent Effect on Perceived Urge Intensity

A factorial ANOVA suggested there was a significant main effect of the runoff period (F(1,75) = 16.350, p = 0.0002) but not group allocation (F(2,75) = 0.931, p = 0.40) on urge intensity. In addition, there was no significant interaction between the factors reported (F(2,75) = 1.05, p = 0.35). In all cases, there was an increase in the average urge intensity reported by the participants before and after the runoff period (see Fig. 4B). At the end of the runoff period, the three conditions finished with a similar intensity score.

However, there was a relationship between urge intensity and sensory response (i.e., time elapsed to urge sensation). While urge intensity at the first reporting point (i.e., before the runoff period) was similar across groups, those in the excitatory condition reported the sensation earlier compared to those in the inhibitory condition (see Fig. 4C).

Normal Bladder Activity May be Computationally Modeled

Under baseline (unmodulated) conditions, the behavior of the model was promising. The model accurately simulated normal bladder function and the associated efferent nervous activity. Specifically, it was able to reproduce the expected filling and voiding cycle (see Fig. 5 top trace). During the filling phase, the efferent projections associated with urine storage (hypogastric and pudendal) were tonically active, while the pelvic projections associated with voiding were silent.

Behavior of the Simulated Bladder Model.

Shown is the bladder behavior (top trace) and associated efferent neuronal activity (middle trace) recorded over a 1000 second simulation run. Blowout boxes contain 5 second windows of the full recorded activity (highlighted in red) during filling (left box) and voiding (right box). Each raster trace was obtained from a randomly selected neuron within the Pudendal, Hypogastric, and Pelvic units (Nneurons = 100 in each case).

In contrast, during void events, bursts of pelvic efferent activity silenced the hypogastric and pudendal projections. This change in activity was associated with a decrease in bladder volume, indicative of a switch in the system from a storage state to voiding (see fig. 5). These results suggest that the model was able to mimic the complex switching behaviors required to mediate the behavior of the lower urinary tract, allowing it to serve as an accurate baseline upon which to test the impact of modulation.

Computational Modeling Confirms the Frequency-Dependent Effects of Tibial Nerve Stimulation

Analyzing the effects of TNS on the computational model revealed an intriguing trend. Under baseline conditions, the model produced a cycle of storage and voiding in accordance with expectations from previously obtained animal data (see [24], Fig. 6A). Applying high-frequency (20 Hz) stimulation to the system completely inhibited the voiding action (according to the expected effect of TNS on bladder behavior). However, if low-frequency (1 Hz) TNS was applied to the system, voiding behavior occurred earlier, and with a greater intensity -causing larger drops in stored volume for a given voiding event (see fig. 6A). These results suggest that there is indeed a frequency-dependent effect of TNS on bladder function, in line with the findings of the study. They may also explain the relatively small effect size observed in the study data. The temporal change of voiding occurrence was relatively small, while the magnitude of the change was much larger. In other words, excitatory TNS did not cause the system to void much earlier, but it did increase the intensity of the voiding events themselves.

Computationally Modeling TTNS.

A. Effects of low (1 Hz) and high (20 Hz) frequency TTNS on simulated bladder function compared to baseline (unmodulated) conditions. Total elapsed sim-cycles equate to 1000s total simulation time. B. Frequency dependent effects of TTNS on bladder contraction. Shown is the average total bladder contraction duration for a 500s simulation period under 21 different stimulation frequencies (0-20 Hz, in 1 Hz increments, Nrepeats = 10). Errorbar= 95%CI. C. Effects of disconnecting specific tibial-nerve projections on total contraction duration (for a 500s simulation period) under low-frequency (1 Hz, i) and high-frequency (10 Hz, ii) conditions. In both cases, baseline behavior represents the unmodulated behavior of the system. = p < 0.05, = p < 0.01, = p < 0.001, = p < 0.0001 Mann-Whitney-U test, Nrepeats = 10 in each condition.

To analyze why this may be the case, the research team explored specific aspects of bladder behavior during neuromodulation. In theory, to drive an increase in voiding activity without considerably changing the temporal nature of the cycle, one or more changes may be occurring:

  • The intensity of bladder contraction is increased

  • The duration of bladder contraction is increased

  • Mechanisms that drive storage are being inhibited to a greater extent.

Analyzing the behavior of the system under a range of different stimulation conditions, the research team discovered a frequency-dependent relationship between TNS and contraction duration. For a given simulation run (500 s total), at low frequencies the average total duration of each contraction was increased significantly, while at high frequencies they were eliminated completely (as voiding became impossible, see fig. 6B). These findings confirm the presence of a frequency-dependent effect of TNS on bladder behavior and indicate that the low-frequency effect may be driven by the extension of bladder contraction, causing a greater intensity of voiding for a given simulation run.

Brainstem Targeting Projections Mediate Low-Frequency Excitatory Effects

As the duration of contraction was highlighted as a measure of interest by the initial computational analysis, we sought to understand the nature of this frequency-dependent effect by analyzing the specific projections that may mediate it. To do so, an analysis of bladder contraction was performed in a range of connective states (where specific tibial nerve projections were severed) during low-and high-frequency stimulation.

A one-way ANOVA revealed a significant effect of the connective state on the duration of contraction during low-frequency (F(5,54) = 6.226, p < 0.001) and high-frequency conditions (F(5,54) = 16.688, p < 0.0001). Post-hoc Mann-Whitney-U tests revealed that significant increases in contraction duration were observed only when brainstem-targeting projections were present in some form (see Fig. 6Ci). However, PMC targeting projections alone were not sufficient to induce any form of change (see Fig. 6Ci). From these results, it is clear that brainstem, rather than spinal cord, targeting projections are necessary for the manifestation of the low-frequency effect, with the PAG likely mediating the central excitatory effect. Although not required to elicit excitatory effects, intriguingly it appears that the inclusion of spinal-targeting tibial projections has a stabilizing effect driving a more consistent up regulation of contraction duration (as the solely brainstem-targeting state had considerably more variance than any other condition, see Fig. 6Ci, brown bar).

High-Frequency Inhibition is Mediated by a Complex Network Interaction

In contrast, the combination of projections that underlie high frequency stimulation may be more complex. Post-hoc Mann-Whitney U testing revealed that three conditions extinguished or significantly inhibited bladder activity during high-frequency stimulation: Full connectivity (spine and brainstem), full brainstem connectivity (both PAG and PMC), and isolated spinal projections (see Fig. 6Cii). Together, these results suggest that a complex interaction between the brainstem and spinal connections underlies high-frequency inhibition, indicating that the effect may be more complex than initially thought. When spinal projections are included, the system becomes resistant to inhibition (since both regions of the brainstem must be targeted to extinguish bladder behavior). These findings require further analysis, but may support the hypothesis of spinal tibial projections as a stabilizing effect.

Discussion & Conclusion

Experimental Study Findings

The results of the present study clearly demonstrate that TTNS has a frequency-dependent effect on bladder function. They corroborate preliminary animal data within the literature [1618] and for the first time demonstrate the presence of frequency dependence in humans. The results suggest that the main effect of the intervention on the healthy population was sensation. Changes in the time taken for participants to experience the urge to urinate were clear, in contrast to the effects on the intensity of the urge that were not significant. Although 20 Hz stimulation appeared to increase the average intensity of the urge after treatment, this is unlikely to indicate the presence of any lingering effect. Rather, it is more likely that this difference compared to 1 Hz stimulation was due to more time elapsing before participants reported any feeling of urgency (see Fig. 4C). In other words, while a similar intensity was reported upon feeling the urge to urinate, the time taken to produce this feeling could be modulated -either reducing it (via 1 Hz stimulation), or increasing it (via 20 Hz stimulation).

More work is needed to assess the therapeutic value of the observed frequency dependence. In addition, the digestion period was probably slightly too long, as several individuals dropped out before receiving neurostimulation. Finally, the effect size of the frequency dependence was relatively small. Despite these limitations, the reported results present a clear picture of a frequency-dependent effect of TTNS in healthy participants and encourage further computational analysis to explore the specific mechanism behind the effect for future clinical studies.

Computational Findings

The results of the computational modeling work provide considerable mechanistic context to the results of the experimental study. The unmodulated behavior of the model was promising, suggesting that it was able to accurately mimic the complex switching behavior employed by the micturition control system[2]. The system produced tonic storage-related efferent activity indicative of a guarding reflex during periods of lower bladder volume which transitioned to void-promoting pelvic activation and silence of storage-related activity at appropriate times. The switch in the circuit state aligned with voiding events in the data. These results suggest the model is able to capture both aspects of micturition (upper circuit activity and lower-bladder behavior) providing an accurate overview of the system for further research. Although the average bladder capacity of the model was lower than would be expected in vivo (simulated capacity was approximately 100-150ml compared to a physiological range of 300-400ml [28]) the accuracy of the high-level behavior of the system allowed the effects of TTNS to be studied in silico.

The impact of TTNS on the system was clear: By varying solely the frequency of stimulation, it was possible to up- or down-regulate simulated bladder activity. This was in line with the hypothesis stated in the article, confirming the presence of a frequency-dependent effect highlighted by the experimental study in silico. An additional benefit of the computational approach was the ability to explain the relatively small effect size observed in the experimental study. As briefly mentioned in the results, the temporal change in voiding behavior was relatively small in the excitatory condition, with voiding occurring now much earlier in response to 1 Hz stimulation. This is in line with the results of the experimental study, where 1 Hz stimulation was found to induce the urge slightly earlier, but not at a considerably greater sensory magnitude. Exploring the implications of the modeling results, it is possible that the main excitatory effect is not mediated by a change in afferent sensation, but rather a change in efferent activity, elongated bladder contraction at the point of void.

The clinical implications of this frequency dependence are particularly important when considering the treatment of non-obstructive urinary retention (NOUR). Currently, the treatment of NOUR remains rather limited. Invasive sacral stimulation is the only approved neuromodulatory treatment option available to patients [29]. Previous attempts to use TTNS to treat NOUR have returned mixed results [30, 31]. However, this may be due to the use of inappropriate stimulation parameters. Previous research did not alter the stimulation parameters of the standard interventions used in the treatment of incontinence and, as such, may have in fact had a deleterious effect on NOUR symptomatology.

Moreover, at present, TTNS is typically applied to treat incontinence prophylactically, with weekly sessions promoting long-term improvement in symptomatology [32, 33]. However, these findings suggest that if TTNS is to be applied in the treatment of urinary retention, it may be more effective to administer the treatment in an acute manner, i.e., at the point of voiding. Doing so could up-regulate the bladder and facilitate unassisted voiding. If shown to be effective when implemented in this way, it could serve as an effective means of reducing the dependency of urinary retention patients on catheters, reducing the rates of urinary tract infections and improving overall quality of life [34, 35].

In addition to providing evidence for therapeutic value in the treatment of retention, the present work may also shed light on the mechanistic underpinnings of this frequency-dependent effect. Our results indicate that low-frequency excitation is mediated by brain stem-specific projections. However, this is not to say that spinal projections are not involved in the process. Spinal cord projections may act as a stabilizing mechanism, ensuring sustained elongation of bladder contraction during low-frequency stimulation. It may be possible to explain these findings in terms of filtering.

We characterize the role of the brainstem as a filter. The inhibitory feedback mechanisms within the PAG and PMC (see Fig. 2 labels 1–4) allow the modules to act as a set of high-pass filters connected in series, preventing efferent void signaling until a sufficient afferent firing rate is reached. Low-frequency TNS can induce an excitatory effect by inhibiting feedback mechanisms within the PAG, reducing the threshold of the filters, thus allowing a smaller afferent signal to pass through to the PMC (the primary micturion control region [2]) triggering void-promoting efferent activity. Inhibition of the PMC alone does not result in excitation, as the first filter in series (the PAG) remains fully operational and likely eliminates any afferent activity before it can reach the PMC. To this end, the inhibitory state of the PMC matters little in this condition. The spinal cord may serve as an additional stabilizing measure, inhibiting the effects of spinal feedback mechanisms that facilitate the guarding reflex.

Although this filtering hypothesis effectively explains the relationship between spinal and brainstem projections during low-frequency TTNS, it cannot explain the complex interaction that may be at play during high-frequency TTNS. It is unclear if this high-frequency inhibition occurs at the spine level, silencing afferent activity before it reaches the filters (Fig. 2, label 5), or at the level of the brainstem, silencing communication between the modules that allows the production of void-promoting efferent activity (Fig. 2, labels 2–4) or the efferent activity itself (Fig. 2, label 1). However, the complex interaction between spinal and brainstem-target projections seen here may explain the paradoxical findings within the literature, suggesting that although TTNS does have an inhibitory effect on spinal activity [8], this alone is insufficient to induce bladder changes in cases of spinal cord injury [4]. Although this may be due to a rewiring of the bladder control system, as there is evidence that early application of TTNS has some effect [36]. To this end, further research is required to explore how the spinal and supraspinal aspects of the system interact to inhibit bladder function during high-frequency TTNS.

Future research should also focus on exploring the clinical efficacy of low-frequency TTNS as a treatment option for urinary retention. Specifically, it would be prudent to determine if administration of the intervention in an acute manner, at a lower frequency, may remedy the mixed results seen previously. Additionally, these findings may contribute to the ongoing development of a wearable neuromodulatory treatment device which may provide a user-friendly method of administering the intervention [37].

Conclusion

Despite some limitations surrounding the methodology of the experimental study and the need for further research in this direction, the present findings nonetheless provide vital evidence for a frequency-dependent effect of TTNS in human beings. Moreover, the suggest a complex and not-necessarily solely brainstem-specific role of the tibial nerve in bringing about this paradoxical excitatory effect. We propose that these findings would lay the groundwork for a number of future research efforts which may elucidate the precise neural mechanisms underlying the complex effects of TNS on bladder function in both healthy and diseased states.

Acknowledgements

This project is funded in part by the University of Edinburgh and the Engineering and Physical Sciences Research Council (EPSRC) grant number EP/W031493/1.

Additional information

Code availability

The custom code used for the analysis in this study is available upon reasonable request from the corresponding author.

Author Contributions

A.M.T and K.N. conceived the study, A.M.T. constructed the computational model, and analyzed the results; M.J. and A.E. collected and processed the bladder data used to fit the computational model and run simulations in the present study; W.J. assisted with the recruitment and running of the study; E.L. created and provided the biophysical bladder simulation used in the computational model; A.M.T and M.J. produced the figures used in the publication; K.N. and S.M. provided research guidance and advice; All authors reviewed the manuscript.

Supplementary Materials

Bladder Urgency Survey

As part of the human study, participants were asked to rate how intense they felt the urge to urinate. This rating was adapted from a validated questionnaire typically used to rank urgency perception as a measure of urological dysfunction [21] (see Supplementary Fig. 1).

Description of the Simulation Pipeline

Each run of the simulation began by specifying a total runtime (in seconds). At t0, to introduce noise, the initial membrane potential of all neuronal units was randomized between the reversal potential for the leak-current (EL, -80mV) and the threshold for action potential (Vth, -50mV). The bladder was assumed to be empty (VB = 0) with an initial inflow rate determined via a stochastic process. At each time step (t), the bladder state was updated and the firing rates of the primary pressure sensitive afferents (F(PB(t))) calculated. This firing rate was then used to update the status of the micturition control circuit and its efferent projections, which would be used to calculate the bladder state for (t + 1). The target values for the time step were then saved and the simulation increased to (t + 1). Upon reaching the final time step of the specified duration (tend) any cached data were collated and exported, and the simulation stopped.

Synaptic Model

Synapses were modeled in terms of conductance using established formulae[38], modified to include opioidergic transmission. The magnitude of the synaptic current (Isyn) was defined as:

where gex/in represents the conductance of classical excitatory/inhibitory currents, and Eex/in their reversal potentials; and gop/Eop representing the action of opioid currents. Modeled synaptic conductance decayed exponentially according to

Synaptic strength was weighted, so that upon presynaptic spiking, the postsynaptic conductance of either excitatory, inhibitory, or opioidergic currents was incremented according to:

where w represents the weighting factor, and δgex/in/op the baseline increase in conductance in response to a presynaptic spike. The weights of all synaptic connections were determined by parameter fitting and constrained within a specific range to allow controlled plasticity within the network.

Computational Model Optimization

Bayesian optimization was used to fit the parameters of the model [39]. To do so, the input to the micturition control circuit had to be carefully controlled. To this end, for the optimization stage, rather than implement a noisy biophysical model, the bladder was represented as a mathematical vector that contained experimentally derived data on internal bladder pressure (PB(t)) and volume (VB(t)). The sample rate of the vector was 50Hz, matching the model. The bladder data used to do so was obtained from previously published work[24] (see Animal Model Surgery and Preparation).

During optimization, the activity of second-order sensory afferents was extracted, downsampled to 10 Hz, and smoothed (1st-order Butterworth Filter, 1 Hz cutoff frequency). These data were compared to the ground-truth neural data associated with the input bladder data (see Animal Model Surgery and Preparation), producing a normalized root mean square error (NRMSE) value representing the difference between the simulated activity and expected activity. The algorithm adapted the weights and parameters of the model such that the NRMSE was minimized (acquisition function - Gaussian hedging, acquisition optimizer - LMBFGS, Ncalls = 200). After running the algorithm, specific weights were manually adjusted to further optimize the efferent activity of the model, ensuring that it was accurate to typical bladder behavior.

Bayesian optimization of model weights and parameters considerably improved the accuracy of simulated afferent activity, bringing it closer to ground truth neural data (NRMSEbaseline=3.367; NRMSEoptimised =1.454). While subsequent manual tweaking did improve the efferent behavior of the model it did lead to a slight increase in NRMSE suggesting it did so at the cost of afferent accuracy (NRMSEfitted = 1.599). Despite the minor decrease in global function fit, this methodology nevertheless improved the efferent behavior of the model and allowed the simulation to better capture the peaking pattern of neuronal firing rates in the afferent arm of the circuit (see Supplementary Fig. 1) providing an accurate foundation for analyzing the effects of neuromodulation. Notably, the optimization algorithm reached its minimum NRMSE plateau in a reasonable number of calls suggesting the global function minimum was true (see Supplementary Fig. 2).

Although the optimized model did not capture every peak in the firing rate or the physiological noise present in the data, the current level of accuracy is more than sufficient to serve as a foundation for exploring TNS, as noise omission would likely only affect the model’s ability to capture very specific network behaviors irrelevant to the present study.

Model Parameters

Where possible these were matched to the parameters described in the original papers describing the neuronal, synaptic, and lower-urinary tract models[23, 27, 38]. Key parameters that were fit included the current mediating tonic activity (Iap), the magnitude of excitatory and inhibitory conductance (gex/in/op), and the upper magnitude of spike-frequency adaption (A). See Supplementary Table 1 for the parameters used in the simulation.

Example urge-intensity survey.

Shown is an example of the survey which was given to participants in printed form.

Model Performance During Fitting Process.

Shown is the overlap between ground truth (blue) and simulated (orange) afferent neural activity for a subset of bladder pressure data. Model performance is shown in a random unoptimized state (left), after 200 rounds of Bayesian optimization (center), and after final manual optimization of model weights (right). The magnitude of the overlap between simulated and ground truth data at each stage shown as Normalized Root Mean Square Error (NRMSE). Smoothed firing rate (1st-order Butterworth filter, 1 Hz cutoff frequency) was normalized against the maximum recorded value in each dataset.

Bayesian Optimization Convergence Plot.

Shown is a convergence plot detailing the Normalized Root Mean Square Error (NRMSE) of the ground truth vs. simulated afferent activity data at each phase of the optimization process.

Final Neuronal Model Parameters.

Where possible parameters were matched to the original specifications of the neuronal/synaptic model. Parameters that were altered by the model fitting process are marked as *.