Interleaved intramuscular stimulation with minimally overlapping electrodes evokes smooth and fatigue resistant forces

Objective. It is known that multi-site interleaved stimulation generates less muscle fatigue compared to single-site synchronous stimulation. However, in the limited number of studies in which intramuscular electrodes were used, the fatigue reduction associated with interleaved stimulation could not consistently be achieved. We hypothesize that this could be due to the inability to place the intramuscular electrodes used in interleaved stimulation in locations that minimize overlap amongst the motor units activated by the electrodes. Our objective in the present study was to use independent intramuscular electrodes to compare fatigue induced by interleaved stimulation with that generated by synchronous stimulation at the same initial force and ripple. Approach. In the medial gastrocnemius muscle of an anesthetized rabbit (n = 3), ten intramuscular hook wire electrodes were inserted at different distances from the nerve entry. Overlap was measured using the refractory technique and only three electrodes were found to be highly independent. After ensuring that forces obtained by both stimulation modalities had the same ripple and magnitude, fatigue induced during interleaved stimulation across three independent distal electrodes was compared to that obtained by synchronously delivering pulses to a single proximal electrode. Main results. Contractions evoked by interleaved stimulation exhibited less fatigue than those evoked by synchronous stimulation. Twitch force recruitment curves collected from each of the ten intramuscular electrodes showed frequent intermediate plateaus and the force value at these plateaus decreased as the distance between the electrode and nerve entry increased. Significance. The results indicate that interleaved intramuscular stimulation is preferred over synchronous intramuscular stimulation when fatigue-resistant and smooth forces are desired. In addition, the results suggest that the large muscle compartments innervated by the primary intramuscular nerve branches give rise to progressively smaller independent compartments in subsequent nerve divisions.


Introduction
Neuromuscular electrical stimulation is a general term that refers to the use of electric currents to evoke muscle contractions by inducing action potentials in the motor neurons or axons that innervate the muscles [1]. When it is used with the purpose of eliciting useful and functional motor behaviors like walking or standing in paralyzed individuals, it is referred to as functional electrical stimulation (FES) [2]. So far, FES has found limited clinical use as a means to restore movement to paralyzed muscles due to its inability to produce contractions that can be sustained for a long period of time [3]. Due to the rapid fatigue onset, conventional FES allows achieving motor functions of limited duration that are insufficient to help paralyzed patients regain independent activities in daily life.
One of the hypothesized causes of the rapid onset of muscle fatigue associated with conventional FES is that motor units are recruited in a way that is the reverse of how they are activated during a natural contraction [4]. Physiologically, slow motor units (fatigue resistant) are recruited before the fast (fatigable) ones. Yet, during FES, the natural order of recruitment is reversed because fast motor units are more easily activated by an external electric field than the slow ones due to the larger diameter of the axons that innervate them [2,5,6].
Another hypothesized cause of the rapid onset of fatigue is related to the stimulation frequency. In a physiological contraction, multiple independent motor units are activated asynchronously at low frequencies, with each one of them generating an unfused response. However, the temporal summation of multiple out of phase unfused responses gives rise to an overall fused and fatigue-resistant force [7]. In contrast, during conventional FES, a group of motor units is stimulated synchronously via a single electrode (commonly placed near the nerve entry site, also known as the motor point). The use of a single stimulation site during FES dictates that stimuli is delivered at a relatively high frequency, known as the fusion frequency, in order to evoke a fused contraction (e.g. 50 Hz) [4], which is well above the natural firing frequency of the motor units (10)(11)(12)(13)(14)(15)(16)(17)(18)(19)(20) [8,9] leading to early fatigue of the activated muscle fibers [10][11][12].
Fatigue and ripple are two clinically important properties of artificially evoked forces that are influenced by stimulation frequency. Force ripple is defined as the degree to which a muscle contraction is not smooth. Lack of contractile smoothness arises when the delivered electrical pulses are widely spaced in time that the force twitches evoked by the individual pulses do not fully overlap or when individual twitches with unequal peak amplitudes summate [13]. High stimulation frequencies generate low ripple contractions that fatigue very quickly, while low stimulation frequencies yield fatigue-resistant contractions that contain high ripple [12,14]. Because improvement in one property has to be traded for worsening in the other, clinically desirable forces that are low in ripple and low in fatigue are difficult to obtain with conventional synchronous stimulation.
Interleaved stimulation is an electrical stimulation modality that attempts to combine the fatigueresistance of low frequency stimulation with the low ripple of high frequency stimulation. In interleaved stimulation, multiple electrodes are used to activate independent motor unit groups within a single muscle or a number of agonist muscles. A phase-shifted stimulus is presented to each motor unit group at a frequency that is lower than its fusion frequency. In this way, each motor unit group is stimulated at a low frequency resulting in low fatigue while the net frequency across all motor unit groups is kept high so that low ripple is obtained. It is believed that the low stimulation frequency, the reduced duty cycle of each motor unit group or muscle, and the improved blood flow are mechanisms by which interleaved stimulation delays muscle fatigue [15]. Fatigue benefit of interleaved stimulation over synchronous stimulation has been demonstrated using surface electrodes in humans [16][17][18][19], spinal cord [20,21] and peripheral nerve electrodes [22][23][24][25][26][27] in cats, and epimysial electrodes in frogs [28]. Yet, the enhanced fatigue-resistance often observed with interleaved stimulation could not be consistently shown with intramuscular electrodes. Intramuscularly, Lau et al [29] and Peckham et al [30] demonstrated that less fatigue occurred during interleaved stimulation than synchronous stimulation, while Buckmire et al [31] recently showed that both interleaved stimulation and synchronous stimulation induced the same amount of fatigue. We hypothesize that this inconsistency was due to excessive overlap between the motor unit groups activated by the electrodes within the interleaved protocol [23][24][25]32]. That is, in the study by Buckmire et al [31] it would be possible that almost the same populations of motor units were stimulated by the electrodes, thus causing the individual muscle fibers to be activated at a high frequency rather than at a low frequency. We consider that in order to expose the effect of interleaved stimulation on fatigue, the overlap must be controlled.
While fatigue was used by many authors as the main outcome measure for comparing the benefits of interleaved and synchronous stimulation, ripple was used only by a few. Authors that quantified ripple showed that, for the same frequency, the amount of ripple obtained with multisite interleaved stimulation was always higher than when a single electrode stimulated at a high frequency was used [22,33]. Downey et al [34] showed that with surface electrodes, interleaved stimulation generally produced high ripple, but when optimal stimulation frequency was used, ripple as low as that obtainable with single site high frequency stimulation was achieved. Hughes et al [28] demonstrated that epimysial multisite interleaved stimulation could not attain the same level of force smoothness achievable with single site high frequency stimulation. Because of the tradeoff that exists between ripple and fatigue when stimulation frequency is varied, it was unknown whether or not increasing the stimulation frequency used duringinterleaved stimulation to obtain a force with similar ripple as that developed during synchronous stimulation would undermine its fatigue resistance.
In this work, our objective was to fairly compare the effect of interleaved and synchronous intramuscular stimulation on fatigue by ensuring that forces obtained from both protocols were ripple matched. A subjacent objective was to validate the applicability of a technology we are currently developing in our lab. Achieving multiple intramuscular stimulation sites with standard wired and bulky electrodes would result in highly invasive and complicated implantation surgeries. To address this problem, we are developing very thin (∅ < 1 mm) wireless flexible thread-like implantable micro-stimulators called eAXONs (for electronic AXON) which will allow implantation to be performed in a minimally invasive fashion via injection avoiding complicated surgeries [35][36][37].

Animal preparation and electrode implantation
The in vivo procedure was approved by the local Ethical Committee for Animal Research (CEEA-PRBB) and by the Catalan Government (application number: BPM-18-0028AE (project number: 10 109)). Three adult male White New Zealand rabbits weighing 4.1-4.5 kg were used in three experimental sessions, one animal in each experimental session.
The animal was anesthetized by administering dexmedetomidine (0.1 mg kg −1 ), ketamine (10 mg kg −1 ), and buprenorphine (0.05 mg kg −1 ) intramuscularly. Hair covering the right shank and thigh was clipped using an electric shaver. The animal was then intubated and anesthesia was maintained by 3%-4% sevoflurane. The cephalic vein was cannulated and lactated Ringer solution was administered intravenously. The animal was monitored during anesthesia via pulse oximetry and capnography while the plane of anesthesia was assessed by checking for corneal reflexes. A 5 cm incision was made on the posterior surface of the shank from the popliteal fossa all the way down to the calcaneus. The hamstring fascia was then cut to gain access to the medial gastrocnemius (MG) muscle. The MG was the muscle of choice because: 1) its intramuscular innervation pattern is known [38], 2) it has a relatively simple as well as easily accessible extra-muscular nerve supply, and 3) it is a mixed fiber type muscle. Multisite stimulation in a heterogeneous muscle, like the MG, is expected to produce more ripple than if a homogeneous muscle was used, such as the soleus. So the MG muscle model provided a worst case scenario for testing our hypothesis.
Carefully, the head of the gastrocnemius muscle was pried away using blunt dissection to reveal the site where the MG nerve entered the muscle. Once the motor point of the MG muscle was visualized as depicted in figure 1(a), a hook wire intramuscular electrode (Catalog number 003400160 by SGM d.o.o., Split, Croatia) was inserted approximately 10 mm deep into the muscle parallel to the projected path of the nerve. This electrode will be referred to hereafter as the motor point (MP) electrode. Subsequently, nine additional intramuscular hook wire electrodes (hereafter collectively referred to as non-MP electrodes) were implanted at different predetermined locations away from the MP according to the scheme shown in figure 1(b). To do that, a punctured and wetted 5 cm × 2 cm paper grid with circular markings at 5 mm intervals was overlaid on the exposed surface of the muscle to guide the implantation as shown in figure 1(c). The hook wire electrodes which were housed in a 1.3 cm long 26 G needle were inserted at an approximate 30 • angle with respect to the muscle surface to approximately the same depth (∼1 cm). The nine non-MP electrodes were given an ID based on their position from the MP; proximal-medial (PM), proximal-central (PC), proximal-lateral (PL), central-medial (CM), central-central (CC), centrallateral (CL), distal-medial (DM), distal-central (DC), and distal-lateral (DL). This electrode arrangement, in which a separation distance of 5 mm was maintained between electrodes in rows and 10 mm between those in columns, was chosen because of the oblong shape of the muscle and our attempt to maximize the area of the muscle that could be explored with the nine electrodes that were available. The lower third of the muscle was spared as stimulation in this region evoked very weak force responses. Two small skin incisions about 1 cm apart were made on the thigh, creating a tunnel through which a short piece of silver wire (∅ ∼ 1 mm) was passed and served as a reference electrode. Once implantation of electrodes was completed, the muscle was covered with glycerin to prevent dehydration before the animal´s foot was mounted on the isometric force measurement system.
Upon postmortem examination, the tips of all the implanted hook wire electrodes were found to be entirely located within the MG muscle and not piercing the underlying musculature.

Isometric force measurement bench and stimulator/recorder
A system was implemented to measure isometric forces evoked from hind limb muscles of rabbits in response to electrical stimulation. A schematic drawing of the isometric force measurement setup is shown in figure 1(d). It consists of a force measurement bench and a stimulator/recorder (i.e. electronics).
The force measurement bench consists of: 1) an acrylic baseplate on top of which a 1.5 mm thick iron sheet was screwed in with bolts to provide a strong ferromagnetic surface for the attachment of the rest of the components, 2) a height adjustable ankle steel clamp equipped with a magnetic base and atraumatic 'cushioned' caps that are intended to go over and clamp the tuberosity of the tibia to immobilize it, and 3) a metallic stand with a magnetic base for mounting the force transducer. The idea of using magnets for fixing the components was taken from [39]. The force transducer (STC-10kgAL-S by Vishay Precision Group, Malvern, PA, USA) was securely bolted to the metallic stand which in turn was affixed to the ferromagnetic surface via magnets. The animal was laid on its side with the leg to be tested being on the top. The leg was held above the surface of and parallel to the baseplate via the ankle clamp and then the foot was fastened via cable ties to a rigid plastic footplate which in turn was mechanically connected to the sensing element of the force transducer by a threaded rod at its center. The foot plate served as an interface between the animal's foot and the force transducer. The knee and ankle joint were flexed at approximately 90 • . This ensured that the MG muscle was at the optimum length that produced the strongest response when it is electrically stimulated. The positioning of the foot resulted in a baseline force of approximately 1 N which was later subtracted from forces evoked by electrical stimulation. The distance between the ankle joint and the center of the footplate was approximately 8 cm.
In essence, the stimulator/recorder consisted in a computer controlled 16-channel current source able to generate monophasic pulses with a fixed duration of 250 µs and a programmable magnitude of up to 4 mA and a resolution of about 0.15 µA. Structurally, the core of the stimulator/recorder was a realtime controller (cRIO-9063 by National Instruments Corp., Austin, Texas, USA) that executed real-time algorithms uploaded from a PC laptop (Windows 10, model 532U3C Ultrabook by Samsung Electronics Co., Ltd, Suwon, South Korea) connected to it via USB. Two input-output modules were plugged into the real-time controller: a signal conditioning and acquisition module for bridge sensors (NI-9237 by National Instruments Corp.) and a 16-channel voltage output module (NI-9264 by National Instruments Corp.). The input module NI-9237 was connected directly to the force transducer referred above. The output module NI-9264 was connected to a custom made 16-channel voltage to current converter. The real-time algorithms, which were specific for each one of the stimulation assays described below, were created in the instrumentation software platform LabVIEW (version 2017, by National Instruments Corp.). The LabVIEW platform was also used to record, pre-process and visualize the signal from the force transducer. The signal from the force transducer was acquired at a sampling rate of 1000 samples per second and low pass filtered by an 8th order Bessel filter with a cutoff frequency of 100 Hz (unless otherwise specified in the assays described below).
In all assays, stimulation pulses were monitored with an oscilloscope (TPS2014 by Tektronix, Inc. Beaverton, Oregon, USA); not shown in figure 1.

Electrical stimulation assays 2.3.1. Recruitment curves.
Individual muscle twitches were produced by delivering 250 µs cathodic monophasic pulses with varying amplitude that increased from 0 to 4 mA in steps of 0.1 mA through each of the ten implanted electrodes. These pulses were automatically delivered in series by the stimulator/recorder at a rate of 1 pulse per second. The peak forces of the evoked muscle twitches were automatically determined for each pulse magnitude and plotted to construct the recruitment curve for each stimulation site.
Intermediate plateaus in the recruitment curves were identified whenever twitch peak forces increased by less than 5% over two or more increments in current for MP electrodes and three or more increments in current for the non-MP electrodes [31]. Three parameters were extracted from the recruitment curves collected for each electrode; 1) maximal force defined as the force generated by the muscle in response to maximum stimulation allowed by our generator (4 mA), 2) plateau force defined as the mean force value at the first intermediate plateau region of the recruitment curve, 3) plateau width defined as the difference between the current value at the end of the first linear segment and the beginning of the second linear segment of the recruitment curve. As later explained in the discussion section, we conjecture that the plateau width and force are parameters that can be used to understand the organization of the motor unit groups within the muscle.
A peak twitch force of 0.3 N (approximately 10% of twitch force produced with maximum stimulation through MP electrode) was chosen as an adequate force level at which both overlap and fatigue were tested. This force level was selected because it was attainable even with stimulation of distalmost sites that generated relatively weak responses. Before performing the stimulation overlap test described below, the current amplitude was adjusted for all nine non-MP intramuscular electrodes such that the magnitudes of their evoked twitch forces were matched to the predefined target value (0.3 N). Such adjustment was performed with the help of another algorithm that consisted of sequentially, at a rate of 4 pulses per second, delivering pulses through each electrode. The peak magnitudes of the resulting twitches were graphically displayed and the pulse amplitude on each electrode was manually programmed whenever the peak magnitude of its corresponding twitch significantly deviated from the desired target value. This step was repeated on a trial and error basis until the peak twitch forces obtained from all the electrodes were approximately equal.

Stimulation overlap.
Non-MP electrodes that selectively accessed independent motor unit groups were identified by quantifying the stimulation overlap amongst all of them using the refractory technique [25,26,40]. This technique measures overlap between an electrode pair by comparing the sum of the twitch forces evoked when the two electrodes are stimulated separately with the twitch force generated when both electrodes fire with an interpulse interval of 750 µs between them . This time delay ensures that the second stimulus is delivered during the refractory period of the motoneurons excited by the first stimulus. If the two electrodes were independent, the forces evoked in response to stimulating one electrode immediately right after the other should be the same as the sum of the forces evoked when the two electrodes are stimulated individually. Contrastingly, if the two electrodes were overlapping, then forces evoked when both electrodes are stimulated in rapid succession would be lower than the sum of their individual responses. This is explained as a fraction of the motoneurons depolarized by the stimulus delivered through the first electrode would still be refractory and would be unable to respond to the stimulus passed through the second electrode. The degree to which the magnitude of the integrated force deviates from the sum of the forces obtained when both electrodes are stimulated separately is proportional to the size of the overlapped region. Overlap was expressed as a percentage with values ranging from 0% to 100% whereby 0% corresponded to no overlap and a 100% indicated full overlap. Percentage overlap was calculated for each electrode pair using the following equation: To analyze the stimulation overlap amongst all nine non-MP electrodes in the array, 81 individual sets of pairwise electrode comparisons were performed, 72 of which were a result of testing each electrode against all other electrodes twice (stimulation was delivered in one order during the first test and the sequence was reversed in the second test) and 9 of which are due to testing each electrode against itself (self-overlap). Thus, the overlap analysis generated a 9 × 9 matrix of overlap values which was further analyzed by a k-means clustering algorithm implemented in MATLAB (MathWorks, Natick, MA, USA) to determine the groups of electrodes that were stimulating the same muscle region (high overlap) and those activating the different muscle regions (low overlap) [41].
To our knowledge, there were no previous reports that showed how overlap amongst multiple intramuscular electrodes can be influenced by their location in the muscle and their arrangement with respect to each other. To investigate this, we analyzed the overlap for each row and each column of electrodes in the 3 × 3 electrode matrix (as shown in figure 1(b)). Since each row and each column of electrodes (when considered in isolation) represented a unique electrode placement, we could thus readily demonstrate the effect that electrode positioning has on overlap. The row and column electrodes are referred to herein as horizontal and vertical triads (a combination of three adjacent electrodes) respectively.
The mean overlap for each triad was calculated by first averaging across the six overlap numbers generated from the six different pairwise test permutations possible with three electrodes within a triad (six permutations per three electrodes). At the end of this step, there was a single overlap value for each triad per rabbit. Then, the overlap values for the same triads in the different rabbits were pooled and averaged. Overlap amongst the electrodes used in interleaved stimulation was calculated in the same way.

Fatigue test.
For the three rabbits, the clustering algorithm indicated that the nine non-MP electrodes were grouped into three sets of highly overlapping electrodes. Out of each of the three sets identified, a single optimal electrode that ensured the lowest overlap when combined with electrodes from the other two sets was selected and was subsequently employed in the interleaved protocol.
Two procedures were performed before the commencement of the actual fatigue test. The goal of the first procedure was to deliver single pulses to check if equal twitch forces could be generated from the three selected electrodes when their corresponding current amplitudes-those previously determined from the overlap test-were applied. If the twitch forces from the three electrodes had different magnitudes, current amplitude was adjusted (only fine adjustments were needed). This careful force matching was crucial for minimizing the ripple in the force produced with interleaved stimulation (combining twitches of unequal amplitude results in a tetanic force with high ripple).
The aim of the second procedure was to determine the current amplitude required for the synchronous protocol to elicit a force similar to that produced by the interleaved protocol. First, the force produced with the interleaved protocol was determined by measuring the force response to delivering a 20 s pulse train of interleaved stimulation using the currents amplitudes determined in the first procedure (check following paragraph for stimulation features). Once the force generated with interleaved protocol was known, a series of 20 s pulse trains of 20 Hz synchronous stimulation was applied to the MP electrode and used to adjust the current intensity until the force achieved with the synchronous protocol was very similar to (but slightly lower than 1 ) the one produced by the interleaved protocol. 1 Since it was very difficult to perfectly match the synchronous force to the interleaved force, we decided to intentionally generate a synchronous force that was very similar to, but slightly lower than, the interleaved force. This ensured that not only were the forces in both protocols matched in terms of absolute ripple, but also that the relative ripple for synchronous protocol was slightly higher compared to the interleaved protocol as a way to set up a worst case scenario for testing our hypothesis.
During these pulse trains, force signal was acquired at 500 samples per second and low pass filtered by an 8th order Bessel filter with a cutoff frequency of 50 Hz. The average force magnitude and the RMS force ripple over the last 5 s of the pulse trains were computed automatically by the algorithm. (It was not necessary to adjust the frequency of the synchronous protocol to match ripple as in preliminary assays it was found that synchronous stimulation at 20 Hz consistently produced similar ripple, or higher, than that produced by the interleaved protocol at 30 Hz).
The actual fatigue tests had a duration of 5 min. Muscle fatigue produced during interleaved and synchronous stimulation was assessed by examining the evolution of a quasi-tetanic contraction (i.e. almost fully fused twitches) over time. To produce a quasitetanic contraction with interleaved stimulation, 250 µs monophasic pulses were delivered through the three independent electrodes at 10 Hz each, with a small time shift (33.33 ms) amongst them, resulting in a composite frequency of 30 Hz (Int 30 Hz). The current value determined in the first procedure for each electrode was used. To obtain a contraction with synchronous stimulation that is equivalent (in terms of magnitude and ripple) to that achieved with interleaved stimulation, a continuous train of 250 µs monophasic pulses was applied on the MP electrode at 20 Hz (Syn 20 Hz) using the current amplitude determined in the second procedure. All the fatigue tests were separated by a 10 min period of rest to prevent fatigue/potentiation from influencing results from the subsequent tests. Force signal was acquired at 50 samples per second and low pass filtered by a 4th order Bessel filter with a cutoff frequency of 5 Hz. The order of testing of the different stimulation protocols was randomized. The difference in the patterns of stimulation between interleaved and synchronous protocols is illustrated in figure 2. Electrical stimulation activated only the MG muscle and no signs of spread of stimulation to other muscles were observed.
To determine the type of the motor units activated during the synchronous and interleaved stimulation, the durations of the twitches, recorded with the synchronous and each of the interleaved electrodes for creating the recruitment curves, were analyzed. Individual twitches evoked over current steps between threshold (current strength at which twitch force was first detected) and the strength used to fatigue the muscle were analyzed for each electrode, The twitches were normalized to their peak values, and their duration defined as the time from the outset of the twitch until the force decayed to 10% of its peak value were measured [6]. Muscle fatigability was assessed by measuring the fatigue index (FI) and fatigue time (FT). FI was defined as the ratio between the peak force and the force at 100 s into the test, while FT was defined as the time it took for peak force to drop to 80% of its value. The peak force was defined as the maximum force value reached during the first 100 s into the test. This specific time window (from 0-100 s) was used to determine the peak force because, in one of the fatigue tests, force was recorded up to a 100 s and this allowed us to be consistent in our analysis of data obtained from different tests.

Sihler stain
Following the protocol outlined in [38], at the conclusion of the experiment, rabbits were euthanized with an overdose of sodium pentobarbital (120 mg kg −1 ) and the MG muscle along with the intramuscular electrodes remaining in situ were excised and immediately fixed in 10% formalin for one month. Thereafter, muscles were then depigmented by maceration in a solution of 3% potassium hydroxide and 3% hydrogen peroxide (0.04 ml of H 2 O 2 per 100 ml of KOH) for four weeks. The muscles were decalcified in Sihler's solution I (one part glacial acetic acid, one part glycerin, and six parts 1% aqueous chloral hydrate) for three weeks till the muscles were transparent. The nerves were stained by immersing the muscle in Sihler's solution II (two part Mayer's hematoxylin, one part glycerin, and six parts 1% aqueous chloral hydrate) for two weeks, followed by destaining step (in Sihler's solution I) lasted for 6 h to remove excess stain. The stained muscles were neutralized in 0.05% lithium carbonate solution for 1 h and cleared in glycerin at increasing concentrations (40%, 60%, 80%, and 100%; for one week in each concentration). After clearing, the muscles were stored in pure glycerin in the dark.

Statistical analysis
Two-way analysis of variance without interaction (treatment and animal) was performed to compare the twitch durations obtained from the different stimulation sites. The animal was used as a blocking factor to account for any variability caused by inter-animal differences [42]. A similar approach was used to test for statistically significant difference in mean overlap amongst the different triads as well as to analyze FI, FT and RMS ripple for the synchronous and interleaved protocols. Linear regression analysis was used to analyze how plateau force and width changed as a function of electrode-MP distance. All results are presented as the mean ± standard error of the mean (SEM). All statistical tests were carried out at a significance level of α = 0.05 using SPSS statistical package (IBM, Chicago, IL, USA).  Figures 3(b)-(d) shows the maximum force elicited from the MG muscle in response to maximum stimulation (4 mA) via each electrode whose location in the muscle is displayed in figure 3(a). The bold horizontal lines in the bars correspond to the mean force value at the plateau region(s) of the recruitment curve obtained from each electrode. Stimulation with the MP electrodes evoked the strongest responses out of all the implanted electrodes. Peak twitch forces produced with a 4 mA stimulus pulse delivered through the MP electrodes were 3.89 N, 3.92 N, and 4.31 N in rabbits 1, 2, and 3 respectively. Responses became weaker as the distance between the electrode and the MP increased. Figure 3(e) shows that the plateau force significantly decreased as the distance separating the electrode from the MP increased (p < 0.05), while figure 3(f) shows that the plateau width significantly increased as the distance between the electrode and the nerve entry point increased (p < 0.05). Figure  3(g) shows the intramuscular nerve branching pattern of the MG muscle from rabbit 3 obtained by Sihler staining. In the three rabbits, there were five primary nerve branches that continued to progressively divide to give rise to gradually thinning daughter branches.

Selectivity of stimulation with intramuscular electrodes
The ability of three adjacent intramuscular electrodes (triad) to selectively activate independent regions of the muscle was assessed as a function of their location in the intramuscular nerve pathway and their arrangement with respect to each other. The light gray, gray, and dark gray horizontal triads shown in figure 4(a) occupied the proximal, central and distal regions of the intramuscular nerve pathway respectively and the stimulation overlap amongst electrode members of each horizontal triad is depicted in figure 4(b). Stimulation overlap produced by the proximal triad (59.1% ± 8.1%, mean ± SEM) was significantly higher than the overlap generated by the central (21% ± 5.5%,) and distal triad (28% ± 7%). The fact that the same horizontal triad produced lower stimulation overlap when placed more distally than proximally highly correlated to the wider separation of the primary intramuscular trunks due to  . The light gray, gray, and dark gray vertical triads shown in figure 4(d) were located in the medial, central and lateral regions of the intramuscular nerve pathway respectively and the stimulation overlap amongst the electrode members of each vertical triad is presented in figure 4(e). Post-hoc analysis revealed that the overlap amongst electrodes in lateral triad was significantly higher (63.4% ± 7) than overlap amongst electrodes in the medial triad (38.5% ± 6.5%), but was not statistically different compared to overlap amongst electrodes in the central triad (47.3% ± 9%). The substantially high overlap measured for the lateral triad could be because, according to the Sihler stain, electrode members of this triad lied in the same region as and in alignment with a single primary nerve branch, which was not the case for other vertical triads that happened to align with different nerve branches. A tracing of the branching pattern shown in figure 4(c) is displayed in figure 4(f) where the approximate locations of the electrodes with respect to the five main intramuscular branches are shown. Consistently, stimulation overlap produced by a triad was the lowest whenever its electrodes were aligned with different nerve branches.

Fatigue induced by interleaved vs synchronous stimulation
The results of applying the clustering algorithm to the overlap data obtained from each rabbit are shown in figure 5. Electrodes that were found to be highly independent and were employed in interleaved stimulation were electrodes PM, DL, and DC for rabbits 1 and 2, and electrodes PL, CL, and CM for rabbit 3. The average overlap amongst the interleaved electrodes was 20.5%, 24.6%, and 4.0% in rabbits 1, 2, and 3 respectively. Electrode positions, currents and frequencies used in the synchronous and interleaved stimulation protocols for each rabbit are summarized in table 1. Figure 6(a) shows the time course of force responses produced from a single trial performed in the MG muscle of rabbit 1 during both interleaved and synchronous stimulation protocols. Of the 14 fatigue trials performed across 3 rabbits, 7 trials were Int 30 Hz and the rest were Syn 20 Hz. Force traces obtained from all 14 trials are provided in the supplementary material. It can be clearly seen that force evoked with interleaved stimulation could be maintained for a longer period of time compared to that generated by synchronous stimulation which decayed very rapidly soon after the beginning of the test. Figure 6(b) displays the mean FT calculated for each stimulation protocol. Mean FT for interleaved stimulation (111 s ± 17.7 s) was significantly longer than that measured for the synchronous protocol (12.6 s ± 3.7 s; p < 0.05). Furthermore, the mean FI of the interleaved protocol was significantly larger (0.86 ± 0.03) than that of the synchronous protocol (0.3 ± 0.06) as shown in figure 6(c). Figure 6(d) shows that the mean of the RMS ripple of forces produced during both stimulation protocols were not statistically different (p = 0.96) as expected.
Twitch force recordings for different stimulation sites exhibited different timing features. In figure 7(a), it is shown, for rabbit 3, the average of seven normalized twitches obtained at varying current amplitudes with the synchronous electrode and with the electrodes 1, 2, and 3 employed in the interleaved protocol. The durations of twitch waveforms collected from the synchronous electrode were longer than those obtained with the interleaved electrodes. Figure 7(b) displays a boxplot of all the twitch durations at different current amplitudes obtained from the synchronous and interleaved electrodes across all rabbits. The duration of the twitches was significantly longer for the synchronous electrode (156 ms ± 1 ms) than interleaved electrode 1 (82 ms ± 2 ms), 2 (100 ms ± 2 ms) and 3 (75 ms ± 3 ms). Amongst the interleaved stimulation electrodes, electrode 2 had a significantly slower twitch than did electrodes 1 and 3.

Recruitment characteristics as a function of electrode position
Recruitment properties were characterized for each intramuscular electrode by measuring three aspects of the generated recruitment curves: the maximum force, the plateau force, and the plateau width. We found that the maximum force decreased as the distance between the electrode and the nerve entry site increased. This can be explained by a reduction in the number of motor units that can be recruited as the stimulation electrode moves away from the main nerve branches. Sihler stain indeed revealed that nerve branches were progressively becoming thinner as they travelled distally thus innervating fewer motor units. Plateau width and force, on the other hand, are parameters that we conjecture can be used to understand the organization of the motor unit groups within the muscle. It has been reported previously that the hind limb muscles of mammals are divided by their primary (main) intramuscular nerve branches into discrete motor unit groups called neuromuscular compartments (NMCs) which can be activated independently of the other [43][44][45]. The organization of muscles into NMCs occurs because the axon populations in the different primary nerve branches are independent of each other. This is believed to occur because when a nerve trunk splits into different primary branches, each daughter (primary) nerve branch receives a unique subset of the axons present inside the nerve trunk.  As a consequence, when one of these primary nerve branches is activated at low stimulation levels, only that nerve branch is recruited. Since the axons in the activated nerve branch do not have collaterals in the other unstimulated nerve branches, there is no means by which action potential developed in the activated branch can antidromically invade the nerve branch point and excite the distant unstimulated nerve branches as explained in figure 8.
Since the location of a nerve branch with respect to an intramuscular electrode affects its recruitment threshold, and because of the widespread distribution of the independent nerve branches, their recruitment with intramuscular stimulation occurs sequentially according to their distance from the electrode as Figure 5. Selection of the intramuscular electrodes that were used in interleaved stimulation. (a),(b), and (c) results obtained by applying the clustering algorithm to overlap data collected from rabbit 1, 2, and 3 respectively. The clustering algorithm identified three different sets of highly overlapping intramuscular electrodes denoted as set 1, set 2 and set 3 for each rabbit. Each black circle represents an individual set and inside each circle are the IDs of the intramuscular electrodes that are members of that specific set. Overlap number inside each back circle indicates the mean % overlap amongst its electrode members (i.e. intra-set % overlap) while the numbers outside the circle below the horizontal brackets indicate the mean inter-set % overlap. Note that for rabbit 2, PC electrode was excluded because it was an outlier (i.e. it highly overlapped with electrodes in both sets 1 and 2). This could be due to the possibility that it lied near a main nerve that branched out to a large region of the muscle and resulted in high overlap with several electrodes. Electrode ID encircled in red in each black circle indicates the electrode that was selected from each set and used in independent multisite interleaved stimulation. current amplitude of the delivered pulses is increased. This piecewise recruitment of muscles with intramuscular stimulation has been linked to compartmentalization. The link is explained as follows: at moderately low current amplitudes, axons within the nerve branch closest to the electrode are the first to be activated. When the current amplitude is increased further, more axons within the same nerve branch become recruited causing force to rise linearly. Up until a certain current value beyond which any further increase in current does not produce additional force (i.e. plateau). At this point, the 'first' nerve branch is fully recruited and the current value is still below the activation threshold of the axons of a more distant neighboring branch. When current is further increased, fresh axons from the neighboring branch along with those in the firstly recruited branch are activated producing a second linear rise in force. Intermediate plateaus appeared in 60% of the recruitment curves collected in this study and corresponded to the consecutive activation of several independent spatially distributed NMCs. The step-wise activation of the MG muscle reported here is consistent with previous studies in which intramuscular electrodes were used to stimulate the soleus [46], ankle extensors [47], MG [48] and biceps femoris [49] muscles of the cat, deltoid muscle of Rhesus monkeys [31], and the anterior tibialis muscle of humans [50]. Furthermore, the fact that our overlap test showed that we accessed three minimally overlapping regions further corroborates the compartmentalization hypothesis of muscle innervation.
The width of the intermediate plateaus in the recruitment curves represents the range of current values between the full activation of the 'firstly' recruited nerve branch and the threshold activation of 'secondly' recruited nerve branch. In other words, the wider the separation between two linear segments in a recruitment curve, the greater the degree to which one nerve branch can be activated in isolation without spillover to the other (i.e. selectivity of stimulation) [47].
The primary nerve branches (controlling the different NMCs) were closely packed together proximally but they became more widely spaced distally. And so, as the distance between the IM electrode and MP increased, (i.e. the branches became more separated), the difference between the recruitment thresholds of neighboring nerve branches increased and so did the plateau width of the recruitment curves collected. Our results indicate that stimulation at sites near the MP evoked strong responses but resulted in recruitment with narrow dynamic ranges (current range over which stimulation remained selective was narrow), while at more distal sites dynamic range of recruitment was broadened (high selectivity) yet lower forces were generated.
Only motor axons in the different primary nerve branches are known to be neurally segregated from each other and can be uniquely activated, and whether or not distal higher order nerve branches can be distinctly stimulated in a similar fashion is unknown [51][52][53]. For that reason, it is presumed that the smallest portion of a muscle that can be selectively activated is the one innervated by a primary nerve branch. In this study, we used the plateau force as a measure of the maximum force produced by the 'firstly' recruited nerve branch and quantified it for each electrode position. Remarkably, a  are shown for simplification) at the nerve branch points inside the intramuscular nerve pathway determines the compartmentalization pattern of the muscle. Primary, secondary, and tertiary labels on each segment of the intramuscular nerve tree indicate the nerve branch order. The yellow highlight over an axon indicates that it is active. Two possible axonal branching scenarios are illustrated: (a) Motor axon 1, 2,3, and 4 enter the same primary nerve branch but then axons 1 and 2 enter one secondary nerve branch while axons 3 and 4 enter the other secondary branch until eventually each axon enters a different tertiary nerve branch resulting in four independent tertiary nerve branches. Low intensity stimulation by an intramuscular electrode placed near the 'red' tertiary branch will lead to highly restricted muscle fiber activation. Action potential developed in the 'red' tertiary nerve branch cannot invade the other tertiary nerve branches antidromically since the axons in activated branch are independent of axons in other nerve branches. This axon branching pattern results in a muscle with compartments about tertiary nerve branches. (b) Motor axon 1 and 2 exclusively enter one of the primary motor nerves but then they branch at each successive nerve branch point and send collaterals into all the secondary and tertiary nerve branches non-selectively. Low intensity stimulation by an intramuscular electrode placed near any tertiary nerve branch will lead to a widely spread muscle fiber activation. Action potential developed in the activated branch will spread to all the tertiary nerve branches antidromically since the axons in the activated nerve branch also have collaterals in the other tertiary nerve branches (i.e. they are dependent). This axon branching pattern results in a muscle with compartments about the primary nerve branches. linear regression analysis revealed that the plateau force decreased as the distance between the electrode and the MP increased. This suggests that progressively smaller motor unit groups (NMCs) were selectively activated by the more distal intramuscular electrodes. An interpretation of this is that NMCs innervated by primary nerve branches are further subdivided by the daughter nerve branches into smaller sub-compartments. If indeed the MG muscle is highly compartmentalized, this would mean that a large number of intramuscular electrodes (larger than the number of the primary nerve branches in the muscle) can be implanted without causing a significant loss in their independence which could be beneficial for producing finer, smoother and more fatigue-resistant contractions [21,24].
Evidence of higher compartmentalization has never been reported possibly because the available techniques used to study compartmentalization of axons like glycogen depletion or EMG recording are suitable for extramuscular nerves but not for intramuscular ones. This is because, while extramuscular nerves are readily accessible and can be individually stimulated, access to intramuscular nerves to stimulate them requires surgical dissection which can be so extensive that it could ultimately lead to the destruction of the muscle. We believe that the analysis of plateaus offers a non-destructive approach to infer the segregation of axons especially in the nerves located in the distal regions of the intramuscular nerve pathway.

Stimulation selectivity with intramuscular electrodes
The number and the locations of intramuscular electrodes required to perform multisite stimulation should be determined based on the number and the locations of the NMCs found in a muscle in order to achieve multisite stimulation with minimal overlap. However, knowledge of the compartmentspecific nerve branches and their respective locations is almost impossible to obtain since there is no current imaging technique available that would allow visualization of intramuscular nerves in vivo. This problem is further complicated by the high variability in the branching patterns of intramuscular nerves and possibly compartmentalization patterns of axons inside the nerves for the same muscle across individuals. Thus, in an effort to help guide the selection of sites for the implantation of multiple IM electrodes to obtain relatively selective recruitment, we studied the effect of varying the placements of the electrodes in the intramuscular nerve pathway on the selectivity of stimulation.
Overlap measured for the different triads highly correlated to the distribution of the intramuscular nerves and where the electrode members of the triad lied with respect to them. The horizontal central triad situated near the mid-belly exhibited far less stimulation interaction than did the same triad at a proximal location. As demonstrated via the Sihler stain, the nerve branches were more widely spaced at the mid-belly region than they were proximally and thus nerve branches were stimulated with less interaction. Likewise, the high overlap produced by the lateral triad strongly correlated to the presence of a large vertical primary nerve branch in the lateral region of the muscle. It seems that selectivity of stimulation was poor because members of this triad were aligned with the same nerve.
The results reported here, emphasize the roles of both the spatial arrangement of intramuscular electrodes with respect to each other and their locations relative to the intramuscular innervation in determining the selectivity of stimulation.
According to what we have exposed above in the first section of the discussion (i.e. possible existence of a large number of independent NMCs, well above what the primary nerve branches suggest) we would expect a large number of independent stimulation sites. Different reasons could explain why our overlap analysis did not turn over a high number of independent sites (the clustering algorithm only identified three independent groups of electrodes). One reason could be because of the arrangement of the nine electrodes used. We speculate that the proximal electrodes excited the parent nerves that gave off the daughter nerve branches which were being stimulated by the more distal electrodes leading possibly to a high degree of stimulation overlap. A single horizontal array of adequately spaced intramuscular electrodes spanning the mid-belly portion could have proven to be a more optimal configuration resulting in less stimulation overlap. Another reason could be that because the pulse amplitudes were adjusted for a specific target twitch force (0.3 N) these amplitudes were either too high or too low to selectively stimulate minor NMCs (i.e. above or below the plateau forces).

Fatigue induced by interleaved vs synchronous stimulation
Interleaved stimulation has been shown to induce less muscle fatigue than synchronous stimulation with most types of electrode configurations in both human and animal subjects [15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30]. However, when intramuscular electrodes were used, fatigue reduction commonly seen with interleaved stimulation could not always be achieved. We reasoned that in those studies, overlap caused by the suboptimal placement of the intramuscular electrodes was another factor aside from interleaved stimulation that was simultaneously affecting fatigue. And so, to assess the effect of interleaved stimulation on fatigue without any bias, overlap had to be minimized. Moreover, it is established that although interleaved stimulation produces higher fatigue-resistance compared to synchronous stimulation, ripple was always worse off. This ripple is not just because of the low stimulus frequencies used during interleaved protocol, it is also due to the difficulty in obtaining equal force tension from different stimulation sites in the muscle or different muscles [13,33]. In this study, we compared the fatigue of interleaved and synchronous stimulation after overlap in the former was minimized. At the same time, we wanted to see if the use of a higher total stimulation frequency during interleaved protocol to achieve a force with the same initial ripple as that obtained with the synchronous protocol would come at the cost of greater fatigue. The results presented here as well as some preliminary results obtained from rats (unpublished work due to methodological limitations) showed that, for contractions of the same smoothness and magnitude, interleaved stimulation with minimally overlapping electrodes induced less fatigue compared to single site synchronous stimulation.
We conjecture that Buckmire et al could not achieve fatigue reduction with interleaved stimulation because the electrodes that were employed highly overlapped (although overlap was not quantified, authors estimated up to 50% overlap between the two intramuscular electrodes) [31]. On the other hand, when a priori knowledge of the compartmentalization pattern of the muscle was used to guide the selection of sites for implanting the intramuscular electrodes, as was done by Liu et al [54], independent compartments were successfully accessed. When those electrodes were employed in interleaved stimulation, fatigue was significantly reduced compared to synchronous stimulation [29].
During the first few seconds of the fatigue tests an enhancement (potentiation) in the force response was observed despite the delivery of pulses of fixed current amplitude, as reported by other authors in [22,55] [56]. As expected, the low frequency stimulation used in this study to fatigue the muscle produced potentiation [56]. Moreover, the low frequency fatigue induced by the interleaved and synchronous fatigue protocols might have been an additional factor that further contributed to potentiation [57,58]. However, the magnitude of this potentiation did not vary between the fatigue trials and so it did not significantly interfere with our ability to match the magnitude of force to the same value at the beginning of the synchronous and interleaved protocols.
The axons of motoneurons may be compartmentalized in the intramuscular nerve pathway not just based on their target compartment but also on their histochemical type. When motor units belonging to the same histochemical class spatially group in a localized region of a muscle, the muscle is referred to as exhibiting 'regionalization' [59][60][61]. The MG of the rabbit is a heterogeneous muscle made up of different proportions of slow and fast motor units with the former being more fatigue resistant than the latter [62]. Previous reports showed a tendency for slow and fast fibers to be spatially segregated in the MG muscle of the rabbit [63] and thus we were concerned that this would cause the electrodes located in different sites to produce contractions with different fatigue properties which could have confounded our fatigue assessment. Generally, the fatigue rate of contractions elicited with stimulation through an intramuscular electrode will depend on the fatigue properties of the predominant type of motor units activated by it. To confirm that the increased fatigue resistance obtained with interleaved stimulation was caused by interleaving stimulation per se, and not because the electrodes employed during interleaved stimulation were favorably inserted in regions consisting predominantly of slow (fatigue resistant) motor units, we measured and compared the contraction properties of motor units activated by the electrodes employed in both synchronous and interleaved protocols across all rabbits. The duration of the twitches elicited from stimulation at different current values with synchronous electrode were significantly longer than those evoked with interleaved electrodes (1, 2 & 3). This indicates that the synchronous electrode recruited more slow (fatigue resistant) motor units than the interleaved electrodes that activated mostly fast (fatigable) ones. Despite the spatial bias created by the nonuniform distribution of slow and fast motor units which have caused the interleaved electrodes selectively activate fast (fatigable motor units), the fatigue resistance of interleaved stimulation was minimally affected. This suggests that the fatigue rate of an electrically evoked contraction depends not only phenotypic type of motor units, but also on the frequency of the stimulation pulses. Our finding is in agreement with results from Vroman et al who showed that muscle fiber composition does not significantly affect fatigue [64]. Similarly, Bickel et al observed no difference between the fatigue induced by repeatedly stimulating a slow and a fast muscle [65]. Evidence of the strong dependency of the rate of muscle fatigue on stimulation frequency has been reported before [11,6666 ] and is reinforced by the force-time traces provided in the supplementary materials of the present study.
We had to increase the stimulation frequency of the interleaved protocol in order to achieve the same ripple as that of the synchronous protocol (30 Hz vs. 20 Hz). This did not adversely affect fatigue resistance. We believe that both the use of lower stimulation frequencies, minimally overlapping electrodes, and access to multiple independent motor unit populations during interleaved stimulation protocol were crucial to achieving such fatigue reduction [50].
Some of the limitations of this study were that the fatigue tests were conducted at a relatively low force level and not over a wide range of forces. Other limitations were the low number of rabbits studied and that the maximum capacity of our current generator was only 4 mA.

Conclusion
We showed that interleaved intramuscular stimulation with minimally overlapping electrodes produces less fatigue than single site synchronous stimulation for the same initial force and ripple. We also showed that overlap amongst the intramuscular electrodes appears to be related to the spatial distribution of the intramuscular nerves. A horizontal row of three electrodes in the central region of the intramuscular nerve pathway was the optimal placement that produced the least amount of overlap because nerve branches were most widely spread apart. Furthermore, we inferred that the MG muscle may be more extensively compartmentalized than previously thought in that motor axons are segregated not just in primary nerve branches but also in higher order nerve branches. The finding that the preferential activation of fast (fatigable) motor units during interleaved stimulation did not hinder its fatigue resistance lends evidence to support that the FES-induced muscle fatigue is caused primarily by the high frequency of stimulation used and to a lesser extent by the type of motor units recruited.