A biomimetic electrical stimulation strategy to induce asynchronous stochastic neural activity

Objective. Electrical stimulation is an effective method for artificially modulating the activity of the nervous system. However, current stimulation paradigms fail to reproduce the stochastic and asynchronous properties of natural neural activity. Here, we introduce a novel biomimetic stimulation (BioS) strategy that overcomes these limitations. Approach. We hypothesized that high-frequency amplitude-modulated bursts of stimulation could induce asynchronous neural firings by distributing recruitment over the duration of a burst, without sacrificing the ability to precisely control neural activity. We tested this hypothesis using computer simulations and ex vivo experiments. Main results. We found that BioS bursts induce asynchronous, stochastic, yet controllable, neural activity. We established that varying the amplitude, duration, and repetition frequency of a BioS burst enables graded modulation of the number of recruited fibers, their firing rate, and the synchronicity of their responses. Significance. These results demonstrate an unprecedented level of control over artificially induced neural activity, enabling the design of next-generation BioS paradigms with potentially profound consequences for the field of neurostimulation.


Introduction
Electrical stimulation is a safe and effective method to interact with the nervous system, which has been successfully used in neuroprosthetics to treat brain disorders [1][2][3], control cardiac activity [4], alleviate chronic pain [5], restore lost sensory and motor functions [6][7][8][9][10][11], and even control autonomic functions [12]. During electrical stimulation, electric currents are exchanged between a stimulator and a pair of electrodes in contact with the body. This process involves the accumulation of charges and electrochemical reactions (oxidation and reduction of ionic species) at the electrode/tissue interface, as reviewed in detail by Merrill and colleagues [13]. Most commonly, charge-balanced rectangular biphasic pulses are used to activate neural structures, whereby a cathodic current-controlled rectangular pulse is delivered first, followed by an anodic pulse of equal charge. This design achieves a good tradeoff between effectiveness and safety [13,14]. Single biphasic pulses are the fundamental units of stimulation and are generally delivered in trains to modulate neural circuit dynamics over time. In this manuscript, we refer to this common electrical stimulation paradigm as 'pulsed stimulation' , whereby a train of cathodic-first biphasic pulses is delivered at frequencies lower than the maximal firing rate of the target neural structures.
Three parameters are used to adjust the effects of pulsed stimulation, namely stimulation amplitude, pulse width, and frequency. Both amplitude and pulse width control the number of neurons recruited by each pulse of stimulation, while stimulation frequency controls the rate at which those neurons are recruited. These parameters are tuned on an application-specific basis [8,15], often with the goal of recruiting neural structures in a way that most closely matches their natural idiosyncratic activity [8,9,[16][17][18][19]. However, pulsed stimulation activates the entire neural population in synchrony with each stimulation pulse [20], with a maximum temporal spread between recruited fibers dictated by the pulse width, which typically ranges from a few hundreds of microseconds to a few milliseconds at most (i.e. different fibers will be recruited at different times, but recruitment will be confined to the duration of a pulse). Such synchronous activity is not commonly observed in vivo, where neural activity is often largely asynchronous, driven in part by the probabilistic nature of action potential (AP) generation in sensory organs, such as muscle spindles [21] or retinal cells [22], and in part by the stochastic nature of synaptic transmission [23].
The high synchronicity of neural activity resulting from pulsed stimulation can have negative consequences in several neurostimulation applications. For instance, in the case of functional electrical stimulation (FES), where stimulation is used to induce muscle contractions, the synchronicity of the response has been shown to induce more jerky movements and higher muscle fatigue [9,24,25]. Similarly, in the context of tactile feedback restoration for amputees, it has been hypothesized that the synchronous neural activity induced in populations of cutaneous afferents accounts for the difficulty in producing natural tactile percepts [16,20].
Generating biomimetic, asynchronous neural activity with electrical stimulation remains an unsolved challenge. A recent study proposed replacing the cathodic and anodic phases of each stimulation pulse with a high-frequency burst of shorter pulses (approximately 100 µs, delivered at more than 7 kHz), with a longer total duration [26]. Using this approach, Zheng et al showed that it is possible to spread motor-axon recruitment over a window of approximately 2 ms, thus inducing motor responses that are somewhat less synchronous compared to pulsed stimulation (see table 1 of Zheng et al [26]). However, while this study demonstrates a first step toward inducing asynchronous and controllable neural activity, natural firing patterns remain significantly synchronous, and the proposed approach is unlikely to enable spreading neural activity over periods of time much longer than a few milliseconds, since delaying the delivery of the anodic phase for too long will increase the chances of tissue damage [13].
Here, we first studied how neural activity induced during pulsed stimulation differs from that generated during natural stochastic processes. Using computer simulations, we emphasized that electrically induced activity is not only more synchronous, but also less variable. We then proposed a biomimetic stimulation (BioS) strategy in which each individual stimulation pulse is replaced by a burst of amplitudemodulated high-frequency stimulation. We reasoned that by slowly increasing the stimulation amplitude, a burst of high-frequency pulses would first recruit the most excitable neurons (large-diameter fibers and neurons near the electrode), and then progressively recruit less excitable cells. In addition, we hypothesized that by using high-frequency stimulation (in the kHz range), the recruited neurons would rarely fire more than once per burst, thus allowing the distribution of neural recruitment over indefinitely long bursts, while still being able to control the number of APs. Indeed, high-frequency stimulation waveforms have been shown to elicit transient neural activity under certain conditions, with neurons firing only once to a few times at the stimulation onset and then becoming unresponsive to further stimulation pulses until the stimulation is either stopped or restarted [27][28][29][30]. Neural responses to kilohertz stimulation are notoriously complex, with different types of response regimes for different sets of parameters, including conduction block and tonic firing. Here, our hypothesis was that there would be a sufficiently large parameter range for which only a limited number of APs are elicited to be of practical relevance. Combining computer simulations and ex vivo experiments in nerve fibers, we corroborated these hypotheses, and showed that this strategy allows for the generation of asynchronous, stochastic, yet controllable, neural activity. We argue that the ability to desynchronize neural activity during electrical stimulation could have far-reaching implications in several applications of neurostimulation.

Computational model of a bundle of nerve fibers during electrical stimulation
We simulated a bundle of a heterogeneous population of myelinated fibers. Fibers were modeled as simple myelinated axons with no cell body, and had two components: a biophysical compartmental model, which we used to evaluate the effect of the stimulation on the membrane potential; and a propagation model, which we used to simulate AP propagation.
The biophysical model was implemented using the NEURON simulation environment [31], and was based on the multicompartment cable model proposed by McIntyre and colleagues [32]. This model was shown to accurately predict the excitation properties of mammalian myelinated fibers to both pulsed [32] and high-frequency [29] extracellular electrical stimulation. Ten compartments were used to model the attachment (MYSA), paranode (FLUT), and internode (STIN) fiber segments between successive nodes of Ranvier ( figure 1(a)). Segment dimensions are based on realistic morphological measurements, which were obtained from nine fibers with diameters ranging from 5.7 to 16.7 µm [32]. To generate additional fibers, with diameters between those that were observed, we used a quadratic polynomial fit over the sampled dimensions. Independently of the diameter, fibers were modeled with 41 nodes of Ranvier (441 segments in total).
To simulate AP propagation along a given fiber length, we implemented a propagation model. When an AP was induced in the biophysical model, its propagation velocity was computed and used to estimate the time the same AP would take to travel along the imposed fiber length. In this way, even if the biophysical model was only a few centimeters long (which is enough to evaluate the effect of electrical stimulation on the membrane), we were able to estimate AP propagation along arbitrarily long fibers, without drastically increasing the computation time.
We used this fiber model to create a heterogeneous population of fibers with properties tuned to match the scenario of upper-limb peripheral nerve stimulation. Specifically, we created a population of 100 group-II afferents, with diameters following a normal distribution with mean of 10 µm and standard deviation of 2 µm, and length of 40 cm. Fibers were distributed uniformly in a circle of radius 100 µm, with a random longitudinal shift to avoid lining up all nodes of Ranvier.
To evaluate the effect of electrical stimulation, we considered a point source electrode positioned in the center of the modeled population, approximately at 7.25 mm from the beginning of the fiber, and an infinite, homogeneous, anisotropic medium. We then simulated an electrical stimulation waveform and computed the external potential at each model compartment, using the formula V e(x,y,z) = I 4π √ σyσzx 2 + σxσzy 2 + σxσyz 2 , where I denotes the injected current, σ the anisotropic conductivity of the medium, modeled to match the conductivity of the white matter [33] (σ x = σ longitudinal = 0.6 S m ; σ y = σ z = σ transversal = 0.083 S m ), and x, y, z the position of the considered compartment with respect to the electrode [34].
When evaluating the effect of pulsed stimulation (frequency: 40 Hz; pulse width: 50 µs; amplitude: −45 µA; figures 1(c)-(e)), BioS (repeating frequency: 40 Hz; burst frequency: 8000 Hz; pulse width: 50 µs; amplitudes: −27 µA; figures 2(c)-(e) and constantamplitude bursts (figure 3(c)), we ran simulations lasting 500 ms, in which the stimulation was turned on 15 ms after the start of the simulation; this ensured that all membrane potentials could reach their resting values. Stimulation amplitudes were set to recruit approximately 80% of the entire population with a single stimulation pulse (figures 1(c)-(e)) or with a single BioS burst (figures 2(c)-(e)). Using these amplitudes ensured that the 'stimulation volumes' addressed by BioS and pulsed stimulation were the same. In all BioS simulations, the starting burst amplitude was set close to recruitment threshold (−10 µA). To reduce computation time, simulations performed to characterize the effect of different BioS burst frequencies (range: from 1000 to 8000 Hz in steps of 500 Hz) and amplitudes (range: −10 to −55 µA in steps of 5 µA, figures 2(b) and 3(a)) or to evaluate the effect of constant-amplitude bursts (figure 3(b)) were performed over 200 ms and 50 fibers. Simulations in figures 1-3 were performed at body temperature (37 • C); simulations in supplementary figure S5 (available online at stacks.iop.org/JNE/17/046019/mmedia), were performed at room temperature (23 • C).

Computational model of a bundle of nerve fibers exhibiting natural neural activity
Natural neural activity-which we used to evaluate the ability of a stimulation waveform to produce biomimetic firings-was modeled using homogeneous Poisson processes. For each modeled fiber, the interval between consecutive APs followed a negative exponential distribution with a fixed mean. This widely used approach provides a convenient approximation of the stochastic properties of natural neural firings [35]. When comparing electrically induced and natural firings, we simulated natural neural activity with a mean firing rate matching the one induced by the tested electrical stimulation condition (pulsed stimulation in figures 1(c)-(e), and BioS in figures 2(c)-(e)).

Analysis of the simulated neural activity
The population neural activity was estimated by convolving the induced APs-recorded at the end of the fibers and expressed as trains of Dirac delta functions-with a Gaussian wavelet (duration: 1 ms; standard deviation: 0.33 ms), and by subsequently adding up the results (figures 1(c) and 2(c)). The neural activity was then binned in windows of 25 ms (figures 1(d) and 2(d)), starting 6 ms after the stimulation onset (estimated minimum period of time required for an AP induced below the electrode to propagate along the whole length of a simulated fiber). Three statistics were computed over the extracted responses: the MAV, the peak to peak amplitude, and the mean standard deviation. In addition to these three statistics, we also extracted the average number of times a recruited fiber fired and the overall number of recruited fibers, to assess fiber recruitment.
The MAV measured the magnitude of the response, i.e. the number of APs generated in the entire population (greater values in this statistic indicate higher levels of neural activity, supplementary figure S6). For each extracted window of neural activity, as defined above, we computed the MAV according to the following: Each fiber is composed of two components: a multicompartmental model that describes the physics of the fiber's nodes of Ranvier and of the myelinated segments, which we used to evaluate the effect of an external electrical field on the fiber potential; and a propagation model that simulates the propagation of the action potentials generated in the multicompartmental model along the whole length of the fiber (40 cm). We modeled the effect of electrical stimulation by considering a point source electrode located in the center of the population, and assuming an infinite, homogeneous, anisotropic space with conductivity properties matching those of the white matter. (b) Distribution of the modeled fibers' diameter. (c) Effect of pulsed stimulation (40 Hz, 50 µA, and 50 µs pulse width) on the modeled population, and comparison with simulated natural firings. Panels report, from top to bottom: the simulated electrical stimulation profile; the effect of electrical stimulation on the membrane potential of a representative fiber at the closest node to the electrode; a raster plot reporting the induced spiking activity reaching the end of the modeled fibers (black dots); and the estimated population activity (black curve). In the raster plot, fibers are sorted by diameter, from smallest to largest. The aquamarine dots and curve represent the simulated natural spiking and population activity, respectively. (d) Representative spiking activity induced electrically (black) or naturally (aquamarine) during a binned stimulation period (top), and average (±standard deviation) population activity over all the simulated stimulation periods (bottom). (e) Bar plots representing the mean and standard deviation of the statistics used to compare the electrically induced activity from the natural one, which we computed over the activity binned in stimulation intervals (n = 18 responses, for each condition). From left to right, the statistics are: the mean absolute value (MAV), which we adopted as a measure of the overall level of neural activity; the peak to peak amplitude, measuring the synchronicity of the neural activity (synchronous responses have higher peak to peak amplitudes); the mean standard deviation of the neural responses (no standard deviation is reported); the number of recruited fibers; and the average number of APs induced in the recruited fibers.
where N s is the total number of samples in a window and x i is the ith sample.
The peak to peak amplitude measured the response synchronicity; assuming that a similar number of APs is induced, the summation of multiple synchronous responses would cause higher maximum values (supplementary figure S6). For each extracted window of neural activity, the peak to peak amplitude was measured as the vertical distance between the highest and lowest amplitude points, as follows: Peak to peak = max (X) − min (X) where X = {x 1 , . . . , x Ns } and N s is the total number of samples in a window.
Finally, the mean standard deviation measured the response variability. This was obtained by computing the standard deviation of each sample within the extracted window of neural activity across all recorded responses, and taking the mean across all samples, as follows: ith samples across all N w windows, and N s is the number of samples within a single window.

Rat dorsal rootlets extraction
All procedures were approved by the Veterinary Office of the canton of Geneva in Switzerland and were in compliance with all relevant ethical regulations. Dorsal rootlets were extracted from four adult male Lewis rats, following the procedures described in [36]. Briefly, rat spinal cords were exposed, and nerve roots were cut at the exit foramen. These were then dissected under a microscope to obtain a nerve strand, which was further teased into rootlets approximately 60 µm in diameter (supplementary figure S1). A total of 16 dorsal rootlets were extracted: ten for the

Staining and rootlets imaging
Staining analyses (hematoxylin) were performed over two of the extracted rootlets to evaluate the heterogeneity of the sampled fiber population (supplementary figure S1). Total fiber count and estimation of diameters were performed using ImageJ64. Contour of each axons were drawn using the freehand selection tool. The area of each contour was measured in software. Each axon diameter was calculated from the measured area, approximated as a perfect disk.

Nerve-on-a-chip platform
A nerve-on-a-chip platform was used to record and stimulate the extracted rat dorsal rootlets, as described in the work by Gribi and colleagues [36]. The platform consisted of two stimulation electrodes and eight recording electrodes arranged in a microchannel ( figure 4(a)). The entire device measured 4 × 2.5 cm 2 . The PDMS microchannels had a cross-section of 100 × 100 µm 2 . Recording and stimulating electrodes were encapsulated in PDMS microchannels (10 mm long). Recording electrodes had contact areas measuring 100 × 300 µm 2 , while stimulation contacts were larger, measuring 100 × 600 µm 2 . Wires were soldered to the chip for connection to an external stimulation (IZ2H stimulator controlled via an RZ2 processor, Tucker Davis Technologies) and recording device (AM system differential AC amplifier, bioamp processor, Sequim, USA and an RZ2, Tucker Davis Technologies). To insert the rootlets into the platform, a suture was tied to one end of each rootlet, allowing it to be gently pulled through the microchannel. All recordings were performed with the nerve-on-a-chip platform located in a faraday cage at room temperature (∼23 • C), using the recording electrode number 6 to maximize the signal to noise ratio [36], and a sampling rate of 48828 Hz.

Data filtering
All data were filtered using a 6th-order, band-pass, digital Butterworth filter with cut-off frequencies of 60-2000 Hz. This range was chosen to remove both the power line noise and the stimulation artefacts.
Indeed, since all the tested stimulation conditions used a pulse-width of 50 µs, the duration of the stimulation artefacts was sufficiently short to be efficiently removed without affecting the neural signals (figure 4(b)).

Experimental stimulation parameters
When evaluating the effects of single stimulation pulses, BioS bursts, or constant-amplitude bursts (figures 4-6), a repetition frequency of 0.5 Hz was used. This ensured that consecutive neural responses were not affected by effects of a previous stimulation. Each tested stimulation condition was recorded for a minimum of 40 s (i.e. n ≥ 20 repetitions). To evaluate the different predicted response regimes induced by BioS bursts, we tested three different burst frequencies: 1000, 2000 Hz, and 8000 Hz. 1000 Hz was chosen to evaluate whether using low frequencies would cause individual fibers to fire repeatedly throughout the burst (figures 2(b) and 3). 8000 Hz was chosen to test the hypothesis that fibers would rarely fire more than once at such high frequencies, and whether summation effects would recruit more fibers than in the case of a single pulse. 2000 Hz was chosen as an additional data point to evaluate bursts with intermediate properties.
To define the minimum stimulation amplitude for a BioS burst, we searched for the lowest current necessary to induce an observable neural response with a single stimulation pulse. The exact value was different for each rootlet but remained in the range of 10-15 µA. Similarly, we set the maximum amplitude by increasing the stimulation current until further increments did not elicit visible increases in the neural response. Across the tested rootlets, the maximum current remained in the range of 30-40 µA. In experiments involving single pulse stimulation and constant-amplitudes bursts, these were always delivered at the maximum amplitudes used for the BioS bursts. During amplitude modulation experiments, we tested different maximum BioS burst amplitudes (figure 5). These various amplitudes were always between threshold and saturation, specifically ranging from threshold amplitude to 200% of threshold amplitude.

Experiments using single stimulation pulses, BioS bursts, or constant-amplitude bursts
Neural responses were extracted using a time-window starting 5 ms before each stimulation pulse/burst and lasting 30 ms. As in the case of the simulated signals, to analyze the induced responses we computed the MAV, the peak to peak amplitude, and the mean standard deviation (supplementary figure S6). MAV and peak to peak amplitudes were computed over the whole extracted signals (30 ms long), for all the stimulation conditions. Instead, the mean standard deviation was only computed over the actual duration of the induced response: 0-5 ms after the stimulation onset, in the case of a single pulse; or 0-25 ms in the case of a burst. To allow comparisons across different rootlets, each statistic was normalized with the median value obtained with a single pulse of stimulation.
When assessing the effect of amplitude modulation, we bootstrapped (n iterations = 10 000) the modulation depth of each statistic. This was defined as the difference between the highest and lowest median statistics obtained over the tested amplitudes.  figure 3(c)), and by constant-amplitude high-frequency bursts delivered at 1000, 2000, and 8000 Hz, respectively from left to right, in a representative extracted rat dorsal rootlet (burst duration: 20 ms; pulse width: 50 µs). Dark-gray shaded areas highlight the initial synchronous neural response (−1 to 5 ms from the stimulation onset); light-gray shaded areas highlight subsequent asynchronous responses (5 to 25 ms from the stimulation onset). (b) Scatterplots reporting, from left to right, the median of the peak to peak amplitude, mean absolute value, and mean standard deviation of the neural responses induced by the tested burst conditions, for all (n = 10) the tested rootlets. For each rootlet, statistics are normalized with respect to the median value of single pulse stimulation (highlighted by the dashed lines). Lines connecting different burst conditions indicate the results for different rootlets. Bars represent the median of each statistic across the tested rootlets. Black stars indicate a significant difference between the starred burst condition and 1 (i.e. the median value of single pulse stimulation); p < 0.05, two-sided Wilcoxon signed rank test with Bonferroni correction for multiple comparisons. Cyan stars indicate a significant difference between the starred burst condition and BioS bursts delivered at the same frequency (figure 4(d)); violet asterisks indicate a significant difference between two burst conditions; p < 0.05, two-sided, Wilcoxon signed rank test with Bonferroni correction for multiple comparisons.

Experiments using continuous BioS or pulsed stimulation
Single neural responses induced during continuous BioS or pulsed stimulation (figures 7 and 8) were extracted using time-windows starting at the onset of each BioS burst or stimulation pulse and lasting an entire stimulation period. The extracted responses were analyzed as previously described. Each stimulation condition was tested for 10 s.
Temporal analysis of the MAV (figure 7(b)) was performed by grouping the induced neural responses into 2 s windows. The bootstrapped modulation over time (n iterations = 10 000, figure 7(c)) was computed as the difference between the medians of the MAVs of the neural responses induced during the last 2 and first 2 s of recordings.

Statistics and data analysis
Data acquisition and analysis were not performed blind to the experimental conditions. To avoid contaminating the data with sporadic movement artefacts during continuous stimulation experiments, statistics that were larger or lower than the median of the sampled distribution by four times the standard deviation were excluded from the analysis. No more than 3.4% of the recorded samples were excluded per trial. No additional data were excluded from the analyses. Statistical significance was analyzed using the two-sided, paired, Wilcoxon signed rank test, the two-sided Wilcoxon rank sum test, or a bootstrap test based on the Monte Carlo algorithm resampling scheme (n iterations = 10 000). In all cases, Bonferroni corrections were used to correct for multiple comparisons. No assumptions on data distributions were performed. Additional details about the number of repetitions and the tests used for each experiment are reported in the corresponding figure legends.

Code and data availability
Recorded data and software routines used for data analysis are available from the corresponding authors upon reasonable request. The code for the computational model is available at https://github.com/FormentoEmanuele/BioS.

Pulsed stimulation induces synchronous and regular neural activity
We investigated the ability of pulsed electrical stimulation to induce biomimetic patterns of neural activity in a population of myelinated fibers. We first developed a biophysical model of a bundle (b) Violin plots reporting, from left to right, the bootstrapped (n iterations = 10 000) modulation depth of the peak to peak amplitude, MAV, and mean standard deviation of the neural responses induced by the tested stimulation conditions (computed as the difference between the maximum and the minimum of the medians calculated over the tested stimulation amplitudes), for all the tested rootlets. Values are normalized with respect to the baseline of single pulse stimulation, i.e. the median value of the lowest tested amplitude. Dark grey plots report the statistics for single pulse stimulation, light grey plots for BioS. Orange stars (below the violin plots) indicate a significant modulation (modulation > 0); * * p < 0.005; one-sided bootstrap test with Bonferroni correction for multiple comparisons. Black stars indicate a significant difference between single pulse and BioS; * * p < 0.005; * p < 0.05; two-sided, bootstrap test with Bonferroni correction for multiple comparisons. Violin plots reporting the statistics for each tested stimulation amplitude are reported in supplementary figure S3.
of fibers-composed of a heterogeneous population of group-II fibers, 40 cm in length in an anisotropic medium-and evaluated the effect of electrical stimulation on the fibers' membrane potential (figures 1(a) and (b)). We simulated a train of biphasic, charge-balanced, symmetric stimulation pulses at a frequency of 40 Hz and amplitude sufficient to recruit approximately 80% of the modeled fibers ( figure 1(c)). Each stimulation pulse instantaneously induced an AP in the recruited fibers, at the node closest to the electrode. Since conduction velocity depends on fiber diameter, the heterogenous distribution of diameters (10 ± 2 µm) within the modeled fiber population caused a weak, progressive desynchronization as the APs traveled away from the site of stimulation. After traveling along the entire length of the fiber, the time difference between the earliest and the latest APs reached approximately 12 ms (figure 1(c)). We next evaluated the difference between electrically induced and natural neural responses at a given mean firing rate (figures 1(c) and (d)), simulated using Poisson processes. We estimated nerve activity during electrical stimulation and during natural firing by convolving each AP with a Gaussian wavelet (figures 1(c) and (d)). We then binned the extracted activity into stimulation intervals (25 ms, figure 1(d)) and computed three statistics: the MAV, as a measure of the overall level of neural activity; the peak-to-peak amplitude, as a measure of the synchronicity of the response; and the mean standard deviation, as a measure of the response variability. In addition, we computed the average number of APs per recruited fiber and the overall number of fibers recruited by the stimulation (figure 1(e)). Since we imposed equal mean firing rate, the MAV of electrically induced and natural neural activity was similar ( figure 1(e)). However, the peak-to-peak amplitude was, on average, 251% greater during electrical stimulation, suggesting that stimulationinduced neural responses were much more synchronous than natural activity ( figure 1(e)). Electrically induced neural responses were also characterized by lower variance (figure 1(e)). This was caused by the repeatability of fiber recruitment at each pulse of stimulation, which contrasts with the stochasticity of natural neural activity. These results highlight a sharp discrepancy between naturally occurring neural activity and the neural activity triggered by pulsed electrical stimulation.

BioS: a stimulation strategy to induce biomimetic patterns of neural activity
Next, we designed an electrical stimulation strategynamed BioS, for Biomimetic Stimulation-to induce neural patterns that are more similar to naturally occurring activity. This strategy consists in replacing each stimulation pulse with an amplitudemodulated high-frequency (1-10 kHz) burst of charge-balanced biphasic pulses ( figure 2(a)). Two hypotheses serve as the theoretical basis for this strategy. First, we reasoned that by progressively increasing the stimulation amplitude of each pulse within a burst, fiber recruitment would be distributed over time: large-diameter fibers, or fibers close to the electrode, would be recruited toward the beginning of the burst, as their recruitment threshold is lower; small-diameter fibers, or fibers further from the electrode, would be recruited toward the end of the burst, as their recruitment thresholds is higher ( figure 2(a)). Second, since neuronal membranes cannot follow high-frequency electric field oscillations, but can exhibit transient responses [27][28][29][30], we hypothesized that the recruited fibers would rarely fire more than once per burst. Therefore, each burst would approximately trigger approximately the same number of APs within the target population as a single pulse delivered at the highest amplitude within the burst, while desynchronizing the induced activity by recruiting fibers throughout its entire duration. Finally, we hypothesized that adjusting the burst amplitude, duration, and repetition frequency would enable precise control over the induced neural activity.
To evaluate the proposed stimulation strategy, we first used the computational model described above to characterize the effect of different burst frequencies and amplitudes on fiber recruitment ( figure 2(b)). When using burst frequencies above 2 kHz, simulations indicated that recruited fibers rarely fire more than once. Increasing the maximum burst amplitude had limited influence on the number of times a fiber was recruited but increased the overall number of recruited fibers ( figure 2(b)). Interestingly, increasing the frequency above approximately 5 kHz also increased the number of recruited fibers. At these frequencies, BioS bursts recruited more fibers than a single stimulation pulse delivered at the same amplitude-a phenomenon caused by the build-up of neuronal membrane depolarizations from individual pulses throughout the burst [37,38]. Since the effects of individual sub-threshold pulses are short lived (due to passive membrane properties), only two stimulation pulses occurring in quick succession can interact with each other to produce facilitation. At lower burst frequencies (<5 kHz), such interactions were less significant.
We then assessed the similarity between the neural activity induced during three different conditions: a BioS burst repeated at 40 Hz with amplitude set to recruit approximately 80% of the modeled fibers (burst duration: 20 ms, burst frequency: 8 kHz), simulated natural neural activity with a mean firing rate set to be equal to the BioS condition, and pulsed stimulation (figures 2(c)-(e)). As hypothesized, BioSinduced neural activity displayed lower synchronization compared to pulsed stimulation, with APs distributed throughout the entire burst duration (figures 2(c) and (d)). Indeed, the average peak-topeak amplitude decreased by 40% with respect to pulsed stimulation and was only 16% greater compared to natural activity ( figure 2(e)). In addition, only a slight increase in MAV was observed. Although the total number of recruited fibers was fixed, each recruited fiber fired 1.37 times on average, explaining the small increase in overall neural activity. Finally, BioS bursts induced-responses displayed a markedly higher variability compared to pulsed stimulation, with a mean standard deviation approaching the value obtained using simulated natural neural activity (figure 2(e)).
We next characterized the ability of BioS to induce asynchronous neural activity over a range of burst frequencies and amplitudes. For each tested parameter we compared the peak to peak amplitudes of the induced neural responses to those computed when recruiting the same fibers with pulsed stimulation ( figure 3(a)). We found that BioS induced neural activity that was largely more asynchronous compared with pulsed stimulation for a wide range of burst frequencies (above 2000 Hz) and amplitudes (below −10 µA). BioS parameters can thus be extensively tuned to control the number of recruited fibers ( figure 2(b)), without affecting the ability to induce asynchronous responses.
We modeled the effect of constant-amplitude high-frequency bursts to test whether the amplitudemodulation which characterizes BioS is necessary to induce asynchronous neural activity (figures 3(b)-(c)). We found that constant amplitude bursts induced neural responses with overall similar peak to peak amplitude compared with pulsed stimulation, suggesting that constantamplitude bursts fail to induce the asynchronous neural activity generated by BioS. In particular, we found that each burst of constant amplitude stimulation induced a synchronous response in all recruited fibers within the first few pulses of stimulation ( figure 3(c)). In addition, for some stimulation parameters, a subset of fibers responded to constant-amplitude bursts with additional asynchronous APs after the first synchronous response. This caused the simulated neural response to have larger peak to peak amplitudes compared with pulsed stimulation.
These results support the theoretical foundations of our hypotheses and suggest that BioS stimulation could be used to artificially induce controllable neural activity which is both highly asynchronous and highly variable. These two characteristics are essential for eliciting biomimetic patterns of neural activity.

BioS bursts induce asynchronous stochastic neural activity
To validate the predictions obtained with our model, we designed an ex vivo experiment using explanted rat dorsal rootlets and a nerve-on-a-chip platform [36]. This platform allowed us to electrically stimulate the small explanted population of heterogeneous fibers (100-150 fibers, supplementary figure S1) and to measure the induced activity 1.3 cm away from the location of stimulation (figures 4(a) and (b)).
We stimulated ten dorsal rootlets with either single pulses, or 20 ms long BioS bursts with frequencies of 1000 Hz, 2000 Hz, and 8000 Hz. We chose these burst frequencies to sample the different response regimes induced by high-frequency stimulation (figures 2(b) and 3(a)). We used an interval of 2 s between consecutive repetitions. Figure 3(c) reports the average (±standard deviation) response of one dorsal rootlet to these four stimulation conditions. Single pulses of stimulation activated the entire population of recruited fibers simultaneously, producing a large compound AP in the measured signal. In contrast, BioS bursts recruited fibers throughout their entire 20 ms duration, without inducing an initial synchronized response. To analyze these responses, we computed the same statistics we adopted in our simulations, namely peak-to-peak amplitude, MAV, and mean standard deviation ( figure 3(d), and supplementary figure S2). Results across the ten tested rootlets show that the peak to peak amplitudes of the responses obtained with BioS bursts were significantly lower (median 47% lower) than those obtained with single pulse stimulation, for all tested burst frequencies (figure 3(d), * p < 0.05, two-sided Wilcoxon signed rank test with Bonferroni correction for multiple comparison). Regarding MAV, we found no significant difference between the responses of BioS bursts delivered at 2000 Hz and those induced by single pulses of stimulation (yet, the median was 29% higher). However, BioS bursts delivered at 1000 and 8000 Hz led to responses with significantly higher MAV. Finally, BioS bursts induced responses with drastically higher mean standard deviations ( figure 3(e)).
These results validate our model predictions, suggesting that each BioS burst induces asynchronous and stochastic neural activity, rarely recruiting individual fibers more than once.

BioS burst amplitude-modulation is essential for desynchronizing neural responses
To further tests our model predictions, we evaluated whether amplitude-modulation of a high frequency burst is necessary to desynchronize neural activity. We stimulated the extracted dorsal rootlets with constant-amplitude high-frequency bursts of 20 ms at 1000, 2000, and 8000 Hz, and evaluated the neural responses ( figure 5).
In agreement with our simulations, constantamplitude bursts induced a strong, synchronous response during the first few pulses of stimulation, which resembled the compound AP induced by a single pulse of stimulation ( figure 5(a)). Indeed, we did not observe any significant differences between the peak-to-peak amplitudes of the responses induced by a constant amplitude burst and single pulse stimulation (figure 5(b), * p < 0.05, two-sided Wilcoxon signed rank test with Bonferroni correction for multiple comparison). Bursts at 8000 Hz induced responses with higher peak-to-peak amplitudes compared with those induced at 2000 Hz. This confirmed our modeling results, suggesting that increasing the stimulation frequency above a certain threshold recruits a larger number of fibers due to summation effects. In further agreement with our simulations, the strong, synchronous response induced at the beginning of each constant amplitude burst was followed by multiple weak, asynchronous responses generated throughout the entire burst duration ( figure 5(a)). Consequently, the MAV and the mean standard deviation of the induced responses were significantly higher than those obtained with single pulse stimulation ( figure 5(b)).
These results suggest that constant-amplitude high-frequency bursts induce responses that are only partially desynchronized, supporting our hypothesis that modulating the stimulation amplitude of each burst is essential for inducing an asynchronous, stochastic, and controllable neural response.

BioS burst maximum amplitude controls the number of recruited fibers
We explored whether modulating the maximum amplitude reached within a BioS burst would allow precise control over the number of recruited fibers, the same way modulating pulse amplitude does in the case of pulsed stimulation. We recorded the neural responses of three dorsal rootlets to BioS bursts and single pulses delivered at stimulation amplitudes ranging from 120% to 200% of the threshold amplitude ( figure 6(a)). The first pulse within a BioS burst was also set to the threshold amplitude.
We found that the MAV of the responses to both BioS bursts and single pulses progressively increased with increasing stimulation amplitude (figure 6(b) and supplementary figure S3). Results in figure 6(b) are reported as violin plots, which show summary statistics similar to a box plot (center circle reports the median, thicker bar extends from the 25th to the 75th percentile), while also showing the shape of the distribution as a grey filled area. We found similar average modulation depths-computed as the difference between the maximum and the minimum of the medians calculated over the tested stimulation amplitudes-for the two stimulation conditions ( figure 6(b), across rootlets average depth of 209% and 203% for BioS and single pulse stimulation, respectively). Yet we found statistically significant differences within rootlets (figure 6(b)). In the case of single pulses, we observed a significant increase in peak to peak amplitude in response to higher stimulation amplitudes, consistent with the fact that an increasing number of fibers were being recruited synchronously (figures 5(b), (c) and supplementary figure S3). On the other hand, BioS bursts demonstrated only modest increases in peak to peak amplitude when using higher stimulation amplitudes. Finally, in the case of BioS bursts, the mean standard deviation of the response robustly increased with higher stimulation amplitudes, while only modestly increased in the case of single pulse stimulation where the increase was statistically significant in only 2 rootlets (p < 0.005).
These results demonstrate that BioS bursts can be used to recruit a controllable number of fibers, while maintaining strong desynchronization, and that the higher the number of recruited fibers, the larger the variability of the induced response.

Tradeoff between synchronicity and firing stability during continuous stimulation
Given the limited ability of axons to keep up with high frequency stimulation profiles [27][28][29][30], we reasoned that repeatedly delivering BioS bursts (continuous stimulation) might lead to rapid depression of the neural responses. Based on this assumption, we further hypothesized that interleaving consecutive bursts with a period of no-stimulation (duty cycling) might be necessary to overcome this limitation. To test these hypotheses and evaluate whether BioS can induce sustained asynchronous neural activity, we stimulated 3 dorsal rootlets with 10 s long, trains of BioS bursts with duty cycles of 50%, 70% and 90% and burst repetition frequency of 20 Hz ( figure 7(a)). For comparison, we also stimulated the extracted rootlets with pulsed stimulation delivered at 20 Hz.
All the tested BioS conditions induced neural responses which were significantly less synchronized than pulsed stimulation (lower peak to peak amplitude, p < 0.0001, two-sided Wilcoxon rank sum test with Bonferroni correction for multiple comparisons, figure 7(b)). Using a duty cycle of 50% consistently induced the responses with the highest peak to peak amplitude. This result is consistent with the hypothesis that BioS bursts distribute neural recruitment over their entire duration thereby evoking more synchronized responses when their duration is shortened. This also suggests that neural synchronicity can be modulated by changing the burst duration.
Regarding the MAV, we observed a decrease over time during BioS, which was absent during pulsed stimulation (figures 7(c) and (d)). For all the tested duty cycles and rootlets, this decrease reached a plateau after 8 s at most (figure 7(c), no significant decrease between samples collected 6-8 s after the stimulation onset and those collected 8-10 s after). Stimulation profiles with higher duty cycles generally produced stronger responses at the beginning of the trial compared with lower duty cycle profiles (figure 7(c)), but also suffered from a stronger depression of the responses (figure 7(d)). When BioS was delivered with a duty cycle of 50%, the maximum reduction of the response over the 3 tested rootlets was of 24% (bootstrapped median, n iterations = 10 000), while this increased to 35% and 42% when the duty cycle was incremented to 70% and 90%, respectively ( figure 7(d)). This finding not only suggests that high duty cycle BioS profiles induce stronger neural fatigue, but also that longer bursts initially evoke a greater number of APs in the recruited population of fibers.
Taken together these results show that BioS can induce sustained asynchronous neural activity, while also highlighting a tradeoff between the level of desynchronization that can be achieved and the temporal stability of the evoked responses. Specifically, while high duty cycle BioS bursts will cause the most desynchronization (since the bursts are longer), they will also more strongly attenuate the neural response over time.

BioS repetition frequency regulates the firing rate of the recruited fibers
Finally, we explored whether modulating the BioS repetition frequency would enable control over the firing rate of the recruited fibers, analogously to the case of pulsed stimulation. For this, we tested the effect of BioS and pulsed stimulation at frequencies ranging between 10 and 80 Hz (n = 3 rootlets). For all the tested frequencies, we delivered BioS with a 70% duty cycle and recorded the responses for 10 s. Increments in stimulation frequency induced a progressive increase in the MAV of the responses evoked by both BioS and pulsed stimulation, in all the tested rootlets (supplementary figure S4, * p < 0.05, * * p < 0.005, * * * p < 0.0001, two-sided Wilcoxon rank sum test with Bonferroni correction for multiple comparisons). This modulation, however, was significantly larger in the case of pulsed stimulation, reaching a maximum modulation depth of 422%, as opposed to 224% in the case of BioS ( figure 8(b), * * p < 0.005, bootstrap test with Bonferroni correction for multiple comparisons). Increasing the frequency also led to a steady, yet weak, decrease in the peak-to-peak amplitude of the responses generated by pulsed stimulation, but not by BioS ( figure 8(b) and  supplementary figure S4). Finally, higher stimulation frequencies led to responses with progressively greater standard deviation (supplementary figure S4).
Overall, these results confirm that adjusting the repetition frequency of continuous BioS profiles enables control over the firing rate of the recruited fibers. Furthermore, these results also highlight differences between BioS and pulsed stimulation in the level of modulation that can be achieved.

Discussion
We developed a biomimetic electrical stimulation strategy (BioS) that induces controllable, strongly asynchronous, and stochastic neural activity, thus overcoming the limitations of pulsed stimulation. Instead of single pulses, the proposed strategy uses high-frequency amplitude-modulated bursts as fundamental stimulation units ( figure 2(a)). Two hypotheses underpinned this design. The first was that ramping up stimulation amplitude within a burst would ensure the progressive and desynchronized recruitment of target fibers. Initially, low amplitude pulses would recruit large diameter fibers and fibers very close to the stimulating electrode. Then, as the amplitude of subsequent pulses increases, smaller diameter fibers or fibers more distant from the electrode would be recruited as well. Finally, the last few pulses in the burst, with the highest current amplitudes, would activate the smallest and most distant fibers in the target population. The second hypothesis was that the limited ability of neurons to follow highfrequency stimulation profiles [27][28][29][30] would prevent fibers from firing repeatedly within each burst. Therefore, a BioS burst would induce a similar number of APs as a single stimulation pulse, while spreading recruitment throughout its duration. Here we discuss the significance of our results, some key BioS parameters and how to tune them, and the implications of this novel stimulation strategy for the field of neuromodulation.

BioS induces asynchronous, stochastic, yet controllable, neural activity
We leveraged a computational model and ex vivo experiments to validate the hypotheses supporting the proposed stimulation strategy. We obtained compelling evidence suggesting that BioS bursts induce asynchronous neural activity by distributing fiber recruitment throughout their duration and found BioS configurations that led to similar overall levels of neural activity compared to pulsed stimulation. These findings corroborated our fundamental hypotheses.
Both computational modeling and ex vivo experiments demonstrated that BioS induces stochastic neural activity. We propose that this phenomenon arises from the presence of high-frequency pulses within a BioS burst. Indeed, when using burst frequencies above approximately 2 kHz, neuronal membranes have insufficient time to return to their resting potential between consecutive pulses in a burst. This implies that each stimulation pulse induces a membrane depolarization that affects the fiber's response to the next pulse (as demonstrated in figures 2(b) and 4). Consequently, the exact moment a fiber is recruited depends on its membrane potential at the beginning of a burst, which can vary slightly from burst to burst due to natural membrane potential fluctuations, causing the response of an individual fiber to BioS bursts to be variable.
We showed that both the number of recruited fibers and the induced firing rate can be controlled by modulating BioS amplitude and repetition frequency, analogously to the case of pulsed stimulation. Interestingly, although increasing the stimulation amplitude used for BioS bursts resulted in progressive increases in average neural activity, this was also accompanied by small increases in peak to peak amplitude ( figure 6(b) and supplementary figure S3). This was likely caused by the fact that when burst amplitude increases, the modulation range likewise increases. For instance, at 120% of the threshold amplitude, a BioS burst is modulated from threshold to 20% above threshold in a fixed amount of time. If the amplitude is increased to 180% of the threshold, the burst is modulated from threshold to 80% above threshold in the same amount of time. Since the number of individual pulses within a BioS burst only depends on the intra-burst frequency (fixed at 2 kHz in the ex vivo experiments) and the burst duration, it is unaffected by amplitude modulation. Consequently, the difference in amplitude between consecutive pulses within a burst becomes larger as the modulation range increases, which in turn makes it more likely for two or more fibers to be recruited by the same pulse. This is consistent with the small increase in peak to peak amplitude observed at larger BioS amplitudes.
During frequency modulation, we observed that BioS induced progressively higher average neural activity, confirming that adjusting the BioS repetition frequency allows fine-tuning of the recruited fibers' firing rate ( figure 8). However, the magnitude of this modulation was significantly lower compared with pulsed stimulation. We propose two concurrent mechanisms to explain this observation. First, as seen in figure 7, shorter bursts induce weaker neural activity compared to longer bursts (lower MAV). Therefore, since an increased stimulation frequency is associated with a shortened burst duration, a lower number of APs will be generated during a burst repeated at high frequencies. This is also consistent with the nearly absent modulation in the responses' peak-topeak amplitude observed when increasing the BioS repetition frequency. Indeed, if the same number of APs were induced in a shorter period of time (i.e. when high BioS repetition frequencies are used), a higher peak-to-peak amplitude would be expected, as discussed in the case of amplitude modulation. However, since this was only the case in one of the three tested rootlets (supplementary figure S4), it is likely that shorter bursts induced weaker overall neural activity. We believe that this phenomenon is caused by a lower number of APs in the recruited fibers, rather than a reduced number of recruited fibers. Second, the MAV likely underestimates the overall neural activity when multiple APs are asynchronously induced in a brief period of time. Indeed, when using short burst durations, the probability that the electrical signals originating from two APs occurring in quick succession will interfere with each other increases (e.g. depolarization and hyperpolarization of two fibers canceling each other out at the recording site). As the probability of such interferences increases, the magnitude to which the MAV of the measured potential underestimates the overall induced neural activity likewise increases. Consequently, we hypothesize that the observed MAV modulation during BioS is underestimated, and that this effect is strongest for very high repetition frequencies.
Interestingly, increasing single pulse stimulation frequency led to a weak, but significant decrease in the responses' peak-to-peak amplitude. This result likely stems from the refractory period of the target fibers. Indeed, when increasing the stimulation frequency above a certain threshold, it is likely that the period between two consecutive stimulation pulses becomes shorter than the refractory period of some of the fibers in the population, thus reducing the probability of all fibers responding to each stimulation pulse. This interpretation is consistent with the observed increase in mean standard deviation at higher frequencies, which suggest that consecutive pulses recruit a different number of fibers. These results also point out that stimulating with simple high frequency trains of pulses-at frequencies above the maximum firing rate of the target fibers, but below the kHz range to avoid the types of blocking effects exploited by BioS (e.g. at 300 Hz)-should also lead to asynchronous neural activity. Indeed, since fibers will not be able to fire at each stimulation pulse, they will slowly desynchronize due to their slightly different refractory periods. However, this strategy would elicit uncontrollable neural firings, since the population firing rate would be tied to the maximum firing rate of the fibers, making this approach unviable for most applications.
These results demonstrate that BioS can induce asynchronous, stochastic, yet controllable, neural activity. This level of control over the recruited fibers is essential for generating biomimetic patterns of neural activity and cannot be achieved with pulsed stimulation.

Key BioS parameters and design features
Some BioS conditions induced significantly higher levels of neural activity compared to pulsed stimulation ( figure 4(d)), suggesting that certain BioS parameters caused some fibers to respond more than once per burst. Our modeling results showed that fibers could fire more than once in response to a single BioS burst and suggested that these events largely depend on the burst frequency, with repeated responses almost exclusively observed when using burst frequencies below 2000 Hz ( figure 2(b)). These modeling results were not fully corroborated during ex vivo experiments, where bursts frequencies as low as 1000 Hz rarely induced consecutive firings in the same fibers. We hypothesize that this discrepancy was the result of the different temperatures used in both experiments, with the simulations performed at 37 degrees, and the ex vivo experiments performed at room temperature (approximately 23 degrees). Adjusting the temperature in the computational model brought simulations in line with the experimental results (supplementary figure S5). These results highlight the importance of properly adjusting the BioS burst frequency to minimize the number of times a fiber is recruited within a burst, which is essential to enable precise control over the induced neural activity. Indeed, changes in stimulation parameters during kilohertz stimulation are known to result in different neural response regimes, ranging from conduction block to tonic firing. Making sure BioS is applied with the appropriate parameters is therefore particularly important and is likely to require some application-specific tuning, since the neural population, tissue properties and electrode configuration can impact the exact range of parameters allowing for an optimal response. Nevertheless, it is important to note that even if a fiber is recruited multiple times within a burst, its firing rate would remain controllable (by adjusting the BioS repetition frequency) as long as it remains lower than the maximum firing rate allowed by the fiber's refractory period.
We demonstrated that amplitude-modulation of the high frequency carrier is essential for inducing asynchronous, controllable neural activity. In particular, we found that constant-amplitude highfrequency bursts induced neural activity consisting of both a synchronous response, of equal or greater magnitude to that induced by a single stimulation pulse delivered at the same amplitude, followed by several asynchronous responses. These results suggest that constant-amplitude bursts induce responses that are more synchronous and involve a higher number of APs in the recruited fibers, compared with BioS bursts. This finding is largely consistent with previous high-frequency studies, reporting similar transient neural activations when constantamplitude high-frequency stimulation profiles were used for nerve conduction block [28,39], but minimal activations when the stimulation amplitude was gradually ramped-up [30].
Although we only tested BioS bursts with linearly increasing amplitude, more complex amplitudemodulation profiles might further desynchronize the induced neural activity. To maximize neural desynchronization, each pulse within a burst should recruit the same number of fibers, so that fiber recruitment is evenly distributed throughout the entire burst duration. As can be seen in figure 4, a linear increase in stimulation amplitude leads to inhomogeneous fiber recruitment, with a distribution skewed towards the end of the burst. The shape of the fiber recruitment function depends on several factors, including the conductive properties of the stimulated tissue, the spatial arrangement of the targeted fiber population, and the distribution of fiber diameters. Assuming a homogenous space populated by equally distributed, identical fibers, and a point source electrode, a linear increase in stimulation amplitude during a BioS burst would lead to a cubic fiber recruitment function, i.e. consecutive stimulation pulses within a burst would recruit a cubically-increasing number of fibers (this phenomenon is simply caused by the spherical propagation of the electric field around the electrode). Therefore, to flatten this recruiting function and maximize neural desynchronization, the amplitude-modulation profile should follow a cubic root function. In realistic scenarios, the optimal amplitude-modulation profile should be computed as a function of the target tissue and electrode properties.
In addition to the stimulation amplitude modulation function, fine tuning the range of stimulation currents is also essential to ensure optimal neural desynchronization. Using a minimum amplitude lower than the smallest fiber threshold would reduce the period of time during which fiber recruitment can be desynchronized, as only a portion of the pulses in a burst would effectively recruit any fibers. A similar effect would occur if the maximum amplitude were larger than the highest fiber threshold, whereby a certain number of the highest amplitude pulses within the burst would not recruit any further fibers. Inversely, if a high minimum amplitude is used, the first pulse within a burst would recruit multiple fibers and lead to a synchronized response. Therefore, the minimum stimulation amplitude within a BioS burst should be calibrated to only recruit the neural structures with the lowest thresholds, while the maximum amplitude should be set at the lowest amplitude necessary to recruit the entire target population. An exception to this rule may be the case of very high BioS burst frequencies (e.g. 8 kHz), where significant summation of consecutive stimulation pulses occurs. In this scenario, starting with a subthreshold pulse may allow certain fibers to fire an AP below threshold amplitudes.
Finally, we showed that BioS can be used to induce continuous asynchronous neural activity, but that the stability of the induced response depends on the stimulation duty cycle. In particular, a high duty cycle might lead to strong depression of the induced responses. We hypothesize that the stronger depression of the responses observed when using high duty cycles (i.e. neural fatigue) may be caused by neurobiological effects, such as local imbalances in ionic species or non-linear ion-channel effects similar to those observed in high frequency nerve block. We also found that adjusting the BioS duty cycle allows control over the level of desynchronization, with high duty cycle stimulation inducing more asynchronous responses compared with low duty cycle stimulation.

Fiber diameters, conduction velocity and natural desynchronization
As can be seen in our modeling results ( figure 1(c)), a single pulse of stimulation will lead to progressive desynchronization of induced APs as they travel along a heterogenous population of axons, each with a different conduction velocity (in our model, APs travel 40 cm along fibers varying in diameter approximately from 6 to 15 µm). Depending on the precise application, electrical stimulation might recruit a much larger range of fiber diameters, such as in the case of peripheral nerve stimulation where a large number of different neuron types are represented. This would lead to larger 'natural' desynchronization than seen in our model, where the range of fiber diameters was chosen to reflect the range one might expect to find in a single population of fibers (e.g. type II afferents). Furthermore, application where stimulation happens further from the target (e.g. stimulating afferent fibers in distal leg nerves) would results in further desynchronization. However, although this 'natural' desynchronization may help mitigate the high synchronicity of electrically induced neural activity, it is important to consider two key problems. First, natural desynchronization resulting from differences in conduction velocity will lead to highly patterned and periodic neural activity, unlike BioS, which leads to more variable and stochastic neural activity. Second, while natural desynchronization effects may appear large when observing very heterogenous neural population (composed of wide ranges of fiber diameters), one is typically interested in activating a single population of neurons (e.g. type II cutaneous afferents). Therefore, since the range of diameters within a population will be much more limited, the contribution of this 'natural' desynchronization effect will likewise be limited. This is why we chose to model a limited range of fiber diameters, since the ability to desynchronize a specific population of neurons is likely to be the more relevant criteria for most practical applications.

Study limitations
Our model of the neural effects of stimulation should, in our opinion, be valued for its ability to offer qualitative rather than quantitative predictions which are directly comparable to experimental results. Indeed, we relied on this model to provide guidance about ranges of parameters and to understand the overall mechanisms underlying BioS stimulation. However, considering the simplifications of the model and the large number of available parameters, we do not expect, nor claim, that our model holds quantitative predictive power. In particular, considering the difficulties with predicting the effects of high frequency neural stimulation, it is likely that modest differences between our modeling results and experimental results exist. In fact, we observed and described a number of these differences in this work.
We did not characterize the biological safety of stimulating using BioS bursts. However, considering that the frequencies used in BioS are similar to those used in typical FDA approved neuromodulation therapies for chronic pain [40], we do not anticipate any obvious safety issues. Compared to the high frequency stimulation protocols used in chronic pain therapies, BioS delivers less current overall, since only a few pulses are delivered at maximum amplitude in each cycle. Nevertheless, a full safety characterization would be necessary before applying BioS in clinical settings.
The statistical metrics used to probe the synchronicity, level of activity and variability of the induced neural activity may sometimes fail to accurately capture specific aspects of the responses. Anticipating such 'failure modes' may help interpret these metrics more accurately (e.g. similar to how familiarity with the shortcomings of the mean, median and mode helps in the interpretation of these values). As described in the methods section, the overall level of activity in the target population was estimated using the MAV of the response. In a simple scenario (i.e. no interference), the MAV obtained with ten individual AP, and the MAV obtained with a compound AP made up of ten fibers would be the same. As such, measuring the MAV over a given time window is a direct measure of the total number of APs produced within that time period. However, when APs occur in quick succession, the likelihood of interference increases (where the positive and negative phases of an AP may partially or completely cancel each other out). This effect becomes more likely when a large number of fibers are recruited asynchronously in a short period of time (such as when using high BioS repeating frequencies). Under these conditions, MAV would tend to underestimate the level of activity induced by stimulation.
Synchronicity was gauged using the peak to peak amplitude of the response. This metric captures the height of the largest compound AP, which is directly proportional to the number of fibers firing synchronously. High values of peak to peak amplitude can only be obtained when a large number of APs occur near simultaneously. Indeed, even small temporal shifts would impact the size of the compound AP measured in the aggregate response. However, similarly to MAV, peak to peak amplitude is susceptible to interference effects. Indeed, two large synchronous responses could partially or totally cancel each other out, leading to an underestimation of the synchronicity of the response. Nevertheless, this is unlikely to have played a meaningful role in our results, since the timing required to cause interference prevents a fiber from interfering with itself, therefore making it difficult for any large compound AP to be canceled out (i.e. if two thirds of the fibers are involved in an AP, only the other third could play a role in any interference effect). However, when measuring highly asynchronous responses, it is plausible that interference may have led to underestimation of the level of synchronicity. In addition, an accurate comparison between the synchronicity of different responses using the peak to peak amplitude is only possible when a similar number of fibers is recruited between the tested conditions. As demonstrated by our computational model (figure 2(b)) and exvivo experiments (figure 4), 1000/2000 Hz BioS and constant-amplitude stimulation bursts recruited a similar number of fibers as those recruited by a single pulse of stimulation delivered at the burst maximum amplitude. Therefore, with the exception of the 8000 Hz condition, this condition was largely met.
Finally, the mean standard deviation was used to estimate the variability of neural responses, and therefore the level of stochasticity in the activity. Any changes in the response from one stimulation pulse to the next are captured by this metric. In this case, the most likely effect that may have played a role in our results is the presence of any additional source of variability that may not be associated with changes in neural response due to the stimulation (e.g. temperature fluctuation affecting neural dynamics). Since we assume this type of 'external' variability affects all measurements equally, we expect that this effect only accounted for a very small proportion of the measured variability, and therefore did not play any substantial role in the interpretation of our results.

BioS stimulation of the nervous system
Electrical stimulation of the nervous system is a valuable tool, allowing researchers to deliver information to the nervous system. However, as highlighted in this work, current strategies fail to generate signals that 'speak' the language of the nervous system: the neural activity induced by pulsed electrical stimulation is highly repeatable and synchronous, while naturally occurring activity is often stochastic and highly asynchronous. BioS has the potential to reduce this nonconformity, and thus improve our ability to communicate with the nervous system. In this section, we explore the impact of the proposed stimulation strategy on selected neurostimulation applications.
FES uses pulsed stimulation to control muscle activity. However, because neural responses induced by pulsed stimulation are strongly synchronous, FES induces jerky movements and higher muscle fatigue, compared with physiological recruitments [24,25]. During sustained physiological muscle contractions, motoneurons fire asynchronously at a rate of 6 to 8 impulses per second [9]. At these frequencies, asynchronous firing is essential to produce the smooth force profiles underlying natural movements, by ensuring the generation of a continuous force. Consequently, delivering FES with pulsed stimulation at 6 Hz would lead to uncontrollable spasmodic movements due to the synchronous contraction of motor units. To compensate for this, current FES paradigms use much higher frequencies (between 20 and 40 Hz). This enables the generation of sustained muscle contractions, but at the expenses of increased muscle fatigue [9,24,25]. Using BioS during FES would thus allow a decrease in muscle fatigue and much finer control over the induced muscle contractions, potentially increasing the benefits provided by this technique to patients with motor disabilities.
Another application in which BioS could be helpful is sensory feedback restoration [20]. Indeed, electrical stimulation of cutaneous afferents typically induces unnatural sensations. In fact, evoked sensations are often reported as vibration or fluttering, which would be consistent with the synchronicity of the restored sensory signals. Two recent studies demonstrated the importance of delivering biomimetic patterns of electrical stimulation to restore the sense of touch in trans-radial amputees during control of a myoelectric hand prosthesis [8,17]. Using peripheral nerve stimulation, these studies showed that modulating the frequency and amplitude of pulsed stimulation over time to mimic the natural firing patterns of touch induces more natural sensations. However, the ability to induce highly biomimetic patterns of neural activity was likely hampered by the non-physiological, synchronous nature of the elicited neural responses. Indeed, even if the population firing rate closely matched naturally occurring firing patterns, the induced activity remained time-locked to each stimulation pulse. Therefore, it is possible that using BioS to further improve the level of biomimicry of the induced neural activity would evoke even more natural sensations.
An important consideration for practical applications of BioS is the increased number of stimulation artefacts associated with higher frequency stimulation. Such artefacts can typically be adequately removed (as we did in this study) by using filtering, blanking or multiplexing approaches. Nevertheless, it is possible that certain applications where stimulation is applied very close to recording sites, or where the recordings are particularly sensitive, might be more difficult to perform with BioS. We do not envision significant issues with stimulation artefacts for the types of applications discussed in this section (i.e. invasive sensory feedback and FES).
In this work, we focused on nerve fibers. Nevertheless, BioS would likely elicit the same effects in any type of neuron. Indeed, the mechanisms underlying the proposed strategy only rely on fundamental membrane properties, which are common to most neural structures. Thus, applications involving cortical or subcortical stimulation may also benefit from BioS [41][42][43][44]. However, it is important to note that some applications likely benefit from synchronicity. For instance, during spinal cord stimulation for gait restoration, the efficacy of certain electrical stimulation protocols depends on the magnitude of the postsynaptic potentials induced in targeted motoneurons by the recruitment of proprioceptive fibers [45]. It is to be expected that higher synchronicity in the afferent activity would cause larger post-synaptic changes, making BioS a poor candidate for this specific application.
In conclusion, BioS is an addition to the neurostimulation toolkit which enables researchers to control a new parameter when recruiting neural structures: synchronicity. Inducing synchronous neural responses will no longer be the default choice. Instead, researchers will be able to choose the strategy that best fits their needs and application, thus opening up interesting new opportunities.