Spatially selective stimulation of the pig vagus nerve to modulate target effect versus side effect

Electrical stimulation of the cervical vagus nerve using implanted electrodes (VNS) is FDA-approved for the treatment of drug-resistant epilepsy, treatment-resistant depression, and most recently, chronic ischemic stroke rehabilitation. However, VNS is critically limited by the unwanted stimulation of nearby neck muscles—a result of non-specific stimulation activating motor nerve fibers within the vagus. Prior studies suggested that precise placement of small epineural electrodes can modify VNS therapeutic effects, such as cardiac responses. However, it remains unclear if placement can alter the balance between intended effect and limiting side effect. We used an FDA investigational device exemption approved six-contact epineural cuff to deliver VNS in pigs and quantified how epineural electrode location impacts on- and off-target VNS activation. Detailed post-mortem histology was conducted to understand how the underlying neuroanatomy impacts observed functional responses. Here we report the discovery and characterization of clear neuroanatomy-dependent differences in threshold and saturation for responses related to both effect (change in heart rate) and side effect (neck muscle contractions). The histological and electrophysiological data were used to develop and validate subject-specific computation models of VNS, creating a well-grounded quantitative framework to optimize electrode location-specific activation of nerve fibers governing intended effect versus unwanted side effect.


Introduction
Vagus nerve stimulation (VNS) is FDA-approved to treat multiple indications and presents a significant clinical and market opportunity for future expansion (Groves andBrown 2005, Beekwilder andBeems 2010). VNS is currently approved for the treatment of epilepsy, depression, obesity, and stroke rehabilitation (U.S. Food & Drug Administration 1997, 2015. Clinical trials are evaluating VNS for other indications including heart failure, diabetes, and rheumatoid arthritis (Boston Scientific Corporation 2021, Drewes 2021, SetPoint Medical Corporation 2021. Unsurprisingly, given its potential in multiple therapeutic applications, there is a substantial market for VNS. In 2020, the global market for VNS was estimated at $8 billion, with projected growth to nearly $17.9 billion by 2028 (Verified Market Research 2022).
The vagus nerve (VN) has such a wide range of clinical applications due in part to its complex organization consisting of over 100 000 nerve fibers in humans (Hoffman and Schnitzlein 1961). Notably, these fibers are responsible for multiple sensorimotor functions of the neck and throat, as well as autonomic functions of the cardiac system, respiratory system, and most organs within the abdomen. Sensory afferent fibers from the visceral organs ultimately influence noradrenergic, serotonergic, and cholinergic inputs to the cortex (Groves and Brown 2005, Dorr and Debonnel 2006, Pavlov and Tracey 2012. Somatic motor efferent fibers run through the pharyngeal, superior laryngeal (SL), and recurrent laryngeal (RL) branches of the VN, which innervate the pharyngeal, cricothyroid (CT), and cricoarytenoid (CA) muscles of the neck (Rea 2014, Câmara andGriessenauer 2015).
Motor activation of the neck muscles is the most common side effect in VNS implant recipients (FDA 2017). Contemporary VNS cuff electrodes wrap ∼340 • (dependent on nerve diameter) around the VN, and they indiscriminately stimulate the entire cross section of the VN trunk, including both sensory and motor nerve fibers. Electrical stimulation activates large-diameter motor efferent fibers at lower thresholds than smaller diameter and unmyelinated fibers (Grill 2004, Horch andKipke 2017). Therefore, VNS often elicits therapy-limiting side effects (mediated by large-diameter motor nerve fibers) at lower stimulation amplitudes than required to achieve intended therapeutic effects (mediated by medium-and small-diameter fibers) (Yoo et al 2013). Activation of large-diameter motor nerve fibers is responsible for cough, throat pain, voice alteration, and dyspnea reported in up to 66% of patients (The Vagus Nerve Stimulation Study Group 1995, Handforth et al 1998, De Ferrari et al 2017, Nicolai et al 2020. In a clinical study of VNS to treat heart failure, desired VNS evoked heart rate (HR) responses indicative of cardiac engagement were achieved in only 12% of measurements taken at the 6-and 12 month end points using stimulation parameters (i.e. signal frequency, pulse duration, and output current) similar to those used in the treatment of epilepsy and depression; off-target stimulation of the neck muscles was the principle therapy-limiting-side effect (De Ferrari et al 2017). Although VNS is a promising therapy for treating various disorders, its efficacy can be limited by side effects induced by activation of offtarget pathways.
Recent studies of VNS suggest that smaller, multicontact electrodes may provide better selectivity for isolating desired therapeutic effects from off-target side effects (Tosato et al 2007, Ordelman et al 2013, Ahmed et al 2020, Fitchett et al 2021, Shulgach et al 2021. Specifically, the studies demonstrated that the magnitude of HR changes may be altered by changing the electrode location around the circumference of the VN. While encouraging, the studies focused only on cardiac response and did not quantify the more common side effects (i.e. motor efferent activation of neck muscles). Further, the studies were phenomenological in design and did not document how the underlying organization of nerve fibers within the VN, coined 'vagotopy,' related to the functional outcomes of stimulation.
In this study, we addressed this gap by quantifying the relationship between the anatomical organization of fiber types in the pig cervical VN and the functional outcomes of electrical stimulation, including recordings of evoked neural activity, activation of deep neck muscles, and changes in HR. The pig is the animal model most resembling the human in terms of diameter and fascicular complexity of the cervical vagus, as well as vagal innervation of deep neck muscles (Nicolai et al 2020, Pelot et al 2020, Settell et al 2020. VNS was performed using the six-contact ImThera cuff electrode (LivaNova, London, UK); the ImThera electrode was awarded an FDA Investigational Device Exemption to conduct clinical studies to stimulate the hypoglossal nerve to treat obstructive sleep apnea, indicating that promising results in application of the ImThera device to the VN could be quickly translated to future human studies. We conducted histology on the VN from each animal to trace fascicles from the nodose ganglion (i.e. a location of known functional organization) to the cervical nerve trunk (i.e. where the cuff electrode was placed); we could thus identify functional groupings of fascicles within the cuff electrode that predominantly contained either motor and parasympathetic efferents fibers or sensory afferent fibers from the visceral organs. These physiological and histological data were used to implement and validate individual-specific, anatomically realistic, biophysical computational models of contact location-specific activation. These computational models provide a conceptual and mechanistic framework for quantifying the opportunities and limits of locationspecific fiber activation that govern desired changes in sympathetic/parasympathetic tone versus unwanted activation of efferent motor nerves innervating deep neck muscles.

Materials and methods
The Institutional Animal Care and Use Committee approved all animal care and procedures at the University of Wisconsin-Madison. Six University of Wisconsin minipigs were included in the study weighing 45 ± 3.8 kg (mean ± SD). All procedures described were performed on the right cervical VN in acute terminal studies under isoflurane anesthesia and lasted 10-12 h. We recorded HR, blood pressure (BP), saturation of peripheral oxygen (SpO 2 ), blood glucose levels (BGL), end-tidal carbon dioxide (CO 2 ) via capnography, electroneurogram (ENG) of the VN via longitudinal intrafascicular electrodes (LIFE), and electromyogram (EMG) of deep neck muscles. These recordings were obtained both as data for analyses and to ensure the welfare of the animal.  (Malagodi et al 1989, Nicolai et al 2020. The recording electrodes were cleaned and reused between experiments, with repairs made, or electrodes remanufactured, as needed.

Materials
A DIN 42802 touch-proof connector (Digitimer, Welwyn Garden City, United Kingdom) was soldered to a stainless-steel (SS) wire (Inside Diameter (ID): 0.001 in, Outside Diameter (OD): 0.0055 in, 10-15 cm length), which was soldered to a single strand platinum iridium (PtIr) wire (OD: 0.004 in, ID: 0.002 in) using special flux (Worthington Stainless Steel Flux, Mfr. Model #331929, Columbus, OH). Heat shrink tubing was applied to the connections and silicone was applied to the ends of the heat shrink tubing to protect the solder joints from adverse environmental effects and to stabilize the position of the protective tubing. A recording window was created in the PtIr wire by removing approximately 1 mm of insulation using a thermal wire stripper (Patco Services Inc. PTS-30S). This recording window was approximately 2 cm distal to the SS/PtIr solder joint. Using 7-0 silk suture, a surgeon's knot was tied around the PtIr wire 2 mm from the window (toward the SS/PtIr solder joint). The suture was covered in Quick-Sil (WPI, Sarasota, Florida) to serve as a mechanical stop and ensure that the exposed electrode window was within the nerve.

Stimulating device: ImThera
The ImThera is a medical device that features six independent stimulating electrodes or 'contacts,' initially designed to stimulate the hypoglossal nerve (figure 1). To avoid ambiguity, the ImThera will be referenced as a 'device' or 'cuff ' throughout this paper, while the six independent stimulating electrodes will be referred to as 'electrodes' or 'stimulating contacts.' When implanted, the device forms a cuff around the nerve, allowing its contacts to have axial and longitudinal separation along the nerve (figure 1).

Validating electrochemical properties of recording and stimulating devices
As outlined by Wilks et al (2017), voltage transient (VT) analysis, electrochemical impedance spectroscopy (EIS), and cyclic voltammetry (CV) were conducted to document the properties of the electrodes (see supplementary material: Electrochemical Testing and Validation). Briefly, EIS was performed on all devices prior to each surgery for quality control. LIFE devices with an impedance below 10 kΩ at 1 kHz were deemed suitable and used in the experiment. All ImThera EIS measurements were stereotypical of platinum-iridium electrodes (see supplementary material: figure 1). Additionally, VT and CV analyses were conducted upon receipt of the ImThera to establish safe stimulation limits (oxidation/reduction and charge injection limits) (supplementary material: Electrochemical Testing and Validation).

Surgical methods: general overview
All animals (n = 6) were sedated using an intramuscular injection of Telazol (6 mg kg −1 ) and xylazine (2 mg kg −1 ). The animals were then intubated and anesthetized using 1%-2% isoflurane. Depth of anesthesia was routinely assessed (via palpebral reflex or nasal septum pinch) and drug dosages were appropriately adjusted to ensure proper sedation levels. Respiration was controlled via a ventilator set to 13-18 breaths per minute to keep end-tidal CO 2 levels within 34-45 mmHg. Normal saline was administered continuously via intravenous (IV) with analgesia (Fentanyl, initial 5 µg kg −1 bolus, 3-5 µg kg h −1 maintenance). HR was measured via EKG. Pulse rate and SpO 2 were measured via an infrared probe (a) Depiction of the surgical preparation in which the right cervical vagus nerve (VN) is exposed. (b) Depiction of pertinent anatomy. The superior laryngeal (SL) nerve branches into the internal superior laryngeal (ISL) nerve that innervates the cricoarytenoid (CA) muscle and the external superior laryngeal (ESL) nerve that innervates the cricothyroid (CT) muscle. The recurrent laryngeal (RL) nerve also innervates the CA. 'Direct' activation of the neck muscles occurs via activation of motor nerve fibers from the stimulating electrode, branching at the RL and innervating the CA. 'Indirect' activation of the neck muscles results from current escaping the stimulating cuff, traveling through fluid and surrounding tissue and activating motor nerve fibers in the SL, which causes contraction of the CA and CT. This figure is not drawn to scale, but the representation of a longer RL than SL is accurate and an essential factor affecting EMG latencies. (C) (Left), Transverse plane depiction of the six-contact ImThera device on the VN. The caudal direction points out of the page, towards the reader. (c) (Right) Additional perspective of the ImThera cuff, showing the axial and longitudinal orientation of stimulating contacts. placed on either the tongue or lip. BP was measured via an invasive probe (Millar) placed in either the left or right femoral artery. The final three animals additionally had their BGL values monitored every 30-60 min. Animals were euthanized using IV KCl (4 mEq kg −1 ).
An incision was made via electrocautery lateral (right) of midline, between the mandible and sternal notch, to access the right VN. In two animals, the omohyoid muscle had to be transected to reach the carotid sheath. Blunt dissection was used to isolate the VN from the carotid artery within the carotid sheath. The VN was exposed from the inferior (nodose) ganglion to approximately 10 cm caudal of this point. The length of the surgical window was intended to maximize the distance between stimulating and recording electrodes to allow separation between the stimulus artifact and the fast-traveling electrically evoked compound action potentials (eCAPs). Care was taken to avoid applying mechanical forces on the VN while other structures were isolated using surgical vessel loops. The surgical window was kept open throughout the procedure. Pooled liquid in the surgical window was removed periodically with manual wicking. Saline-soaked gauze was placed around exposed portions of the VN, and saline spray was applied to the surgical window as needed to keep all structures tissues hydrated.

Surgical methods: ImThera implantation
The proximal edge of the ImThera device was placed on the right VN caudal to the nodose ganglion at a distance of 1.1 ± 0.80 mm. A 7-0 silk suture was placed around the device to keep it stationary. Effort was made to minimize mechanical forces on the VN from the device.

Surgical methods: life implantation
Five LIFE devices were placed within the VN caudal to the ImThera cuff to record ENG, caudal to the ImThera cuff. These devices were implanted longitudinally (caudal to cranial) such that the entirety of the recording window was parallel to the course of the nerve. The most cranial LIFE served as a reference electrode for the remaining four LIFE devices. The recording electrodes were staggered in the cranial/caudal direction to discern spatiotemporal aspects of the ENG signals. Cranial/caudal staggering ensured signals were confirmed as neural in origin if temporal differences were observed among recording electrodes. Medial/lateral staggering enabled sampling and later analysis of ENG signals from various locations within the nerve. The distance between the stimulating electrode and the recording electrodes was maximized within the surgical window to reduce the contamination of the neural signals with stimulation artifact. This resulted in an average distance of 7.93 ± 1.75 cm from the caudal edge of the cuff to the most cranial LIFE and 0.79 ± 0.38 cm from the most cranial (which was used as the reference electrode) to the most caudal LIFE.

Surgical methods: EMG implantation
EMG needle electrode pairs were placed in the CT, CA, and pharyngeal constrictor muscles.

Experimental protocol
We conducted stimulation trials in the 'intact' condition, followed by multiple control conditions. In all cases we stimulated to approximate clinical parameters, using a symmetric biphasic rectangular pulse with 0 ms interphase delay, 0.2 ms per phase, cathodic phase first, and 0.4 ms total duration, mimicking common clinical parameters (supplementary material: figure 2) (Koo et al 2001). All stimulating pulses were delivered in a monopolar configuration to avoid potential interactions between simultaneously active cathode/anode electrode pairs along the nerve. We placed an 18-gauge, stainless steel needle (152 mm 2 surface area) within the thoracic wall to serve as the counter electrode. The stimulus was delivered through one of six contacts of the ImThera cuff at 25 Hz with random amplitudes from 50 to 5000 µA. We paused for 60 s between each amplitude to allow for physiological states to return to baseline. The range of randomized stimulation amplitudes was adjusted on a case-by-case basis based on physiological responses (e.g. in some animals, severe bradycardia was evoked during stimulation at high amplitudes, see section 3). Stimulation was delivered sequentially through electrode contacts 1 through 6. For each contact and stimulation amplitude, a total of 750 pulses (i.e. 30 s of stimulation) per amplitude per stimulation contact were delivered in the intact condition. Subsequently, animals were stimulated under the following conditions to confirm the source of recorded signals: (a) neuromuscular junction (NMJ) block (vecuronium, IV, 0.15 mg kg −1 bolus, 0.15 mg kg −1 h −1 maintenance rate), (b) transection of the SL branch, (c) transection of the RL branch, and (d) transection of the vagus trunk cranial and caudal to the stimulating cuff (double vagotomy). Multiple studies demonstrated that EMG signals can appear as artifacts in neural recordings (Yoo et al 2013, Nicolai et al 2020; therefore, we conducted both NMJ block and transections to assess the origin of recorded signals (supplementary material: Off-Target Activation). After stimulating in the NMJ condition, we waited for the paralytic agent to wear off before proceeding with the remaining conditions. In all these additional conditions (NMJ to double vagotomy), we delivered 25 pulses per amplitude (i.e. 1 s of stimulation) due to limitations in our protocolallotted surgical time.

Histological analysis
Histology of the nerve was conducted to relate vagotopy with stimulation-evoked physiological responses. To ensure that no movement had occurred during the procedure, the placement of the cuff on the nerve (longitudinal location and relative rotation) was recorded and measurements to known landmarks (e.g. nodose ganglion) were made before and after the procedure. Following all functional studies and euthanasia, animals underwent microdissection and histological procedures as previously described (Settell et al 2020). Briefly, the VN was exposed from the nodose ganglion to the RL bifurcation at the level of the subclavian artery. To match histology to relative in vivo locations, tissue dye was used to mark the ventral and lateral edges of the nerve (Bradley Products, Inc., Davidson Marking System, Bloomington, MN), and sections were removed for fixation and processing. All sections were placed in 10% neutral buffered formalin at 4 • C for 24-48 h. Samples then were embedded in paraffin wax and 5 µm thick slices, approximately every 40 µm, were collected and mounted on charged slides (Sakura Tissue-Tek VIP). Slices were stained with Gomori's trichrome and were imaged with a Motic Slide Scanner (Motic North America, Richmond, British Columbia). Measurements acquired during the surgery were used to superimpose the in vivo location and orientation of the cuff onto the histology.

Data analysis 2.2.8.1. Preprocessing ENG and EMG signals
The ENG and EMG data were digitally filtered using in-house developed software (Ludwig Lab 2021) to increase the signal-to-noise ratio and to remove longterm baseline drift without introducing filter ringing from a large stimulation artifact, which could otherwise be misconstrued as neural signals (Lindeberg 1990, Widmann andSchröger 2012). Specifically, we subtracted the median filtered signal (kernel size: 201) from the original signal (Jarske and Vainio 1993). Next, we applied a Gaussian filter (standard deviation for Gaussian kernel: 0.87) and a finite impulse response filter to reject common-mode noise and its harmonics (60, 120, and 180 Hz) with 1.0 Hz bandwidth. An additional median filter (kernel size: 11) was applied to EMG signals to eliminate intermittent spikes present only in EMG recordings.

eCAP magnitude calculation
After filtering, neural signals were averaged across pulse trains and segmented into five time-restricted windows based on known conduction velocities per fiber type (i.e. Aα, Aβ, AΥ , Aδ, and B) (supplementary material: table 1). The onset of stimulation was set as the time zero starting point. The boundaries of the time-restricted windows were calculated by dividing the distance from the stimulating device to the recording electrode by the conduction velocity of each fiber type. For neural signals uncontaminated by stimulation artifact, the root mean square of the measured voltage (V RMS ) within each timerestricted window was calculated to determine the magnitude of the eCAP generated by that nerve fiber type.
If neural signals contained noticeable eCAPs but were still contaminated by stimulation artifact despite the preprocessing steps, the signals were processed additionally to isolate the eCAP (supplementary material: figure 3). First, to find the onset and offset of the eCAP, the contaminated signals were filtered using a Savitzky-Golay filter to smooth the signal and eliminate false minima (Abraham and Golay 1964). Next, a peak-detection algorithm was applied to find the local maximum within the time window for each fiber type. The time points of the closest minima occurring before and after the maximum were classified as onset and offset of the eCAP. The later minimum (offset) could exceed the bounds of the timerestricted window by 25% to allow for the entirety of the signal to be captured. Once the onset and offset of the eCAP were found, these points were superimposed on the filtered signal (without Savitzky-Golay filter). Similar to other studies examining the amplitude of the evoked response in EMG recordings by calculating the integral of the measured signal (Haughton et al 1994, Kent-Braun 1999, Zory et al 2005, we calculated the definite integral of the eCAP. To account for the stimulation artifact, a simplified model of the stimulation artifact was created by forming a line between the onset and offset of the eCAP. The definite integral of this simplified model was subtracted from the definite integral of the signal to give an area measurement of the eCAP (supplementary material: figure 3).

Dose-response curve (DRC) calculation
DRCs were generated for neural, muscle, and physiological responses to randomized stimulation amplitudes. The 'intact' condition was used to analyze: (a) nerve fiber responses measured via ENG, (b) muscle response measured via EMG, and (c) changes in HR and BP from baseline measured via EKG and invasive arterial catheter, respectively.
As described above, the magnitude of the eCAPs was calculated either as V RMS or as area-under-thecurve following the peak-detection algorithm. The amplitude of the EMG response was calculated via V RMS ; as evoked EMG responses occurred outside of the window contaminated by the stimulus artifact (∼5 ms), no stimulus subtraction methods were required to measure the magnitude of the evoked signal. We calculated 95% confidence intervals (CI) for EMG and ENG responses by bootstrapping the sampled signal and calculating the amplitude of the measured eCAP based on the bootstrapped sample estimates (see supplementary material: Interpolation and Statistics). Changes in HR and BP were defined as the maximum change during a stimulation pulse train from the baseline calculated over the interval immediately prior to stimulation (i.e. t − 3 to t − 1 s before stimulation).

Response interpolation and statistical analyses
We examined trends in evoked measurements based on the relative location of the stimulating electrode to the efferent and afferent clusters within the VN. For each animal, the cross-sectional histology of the VN (at the location of the center of the stimulation cuff) was annotated with the locations of the six stimulating contacts (supplementary material: figure 4(a)). To ensure that the impact of electrode location with respect to the efferent mode was consistent across the cohort analyses, the histology/contact figure was rotated such that the efferent-mode fibers were in the top half of each image (12 o'clock position) and the afferent-mode fibers were in the bottom half (6 o'clock position) (supplementary material: figure 4(b)).
In each animal, for each of the six stimulating contacts, we calculated the maximum change in HR and the stimulation amplitudes that elicited 50% of the maximum response (ED 50 ) for Aα and EMG. To determine ED 50 values for Aα and EMG measurements, we fit the neural and EMG DRCs to the sigmoidal Hill equation (Hill 1910) using GraphPad (Prism, San Diego, CA). For each metric (Aα ED 50 , EMG ED 50 , maximal HR change) in each animal, we linearly interpolated 500 values between the neighboring contacts along the nerve circumference, which resulted in 3000 values (500 data points ×6 neighboring contact pairs) per metric, per animal. These 3000 values were plotted around the perimeter of the rotated micrograph (efferent-mode facing the 12 o'clock position) using Python.
To build an aggregate heatmap and calculate statistics across the cohort, a generalized VN cross-section was created utilizing the same efferent/afferent oriented gestalt, and the interpolated values were averaged, based on their relative location, across the cohort. For the cohort maximum change in HR, Aα ED 50 , and EMG ED 50 , we averaged the interpolated arrays across the animals for which we were able to obtain recordings with sufficient resolution (n = 6, n = 4, n = 5, respectively).
Additionally, we conducted a paired t-test comparing the Aα ED 50 values, across the cohort, for the three contacts closest to the efferent cluster versus the ED 50 values for the three contacts furthest from the efferent cluster.

Computational modeling
We used the ASCENT ('Automated Simulations to Characterize Electrical Nerve Thresholds') pipeline v1.1.0 (Musselman et al 2021) to model VNS for all six pigs with the ImThera cuff using individualspecific nerve morphology from histology (supplementary material: figures 5 and 6, table 2). For each pig, we segmented the epineurium and fascicles in the histological cross section at the level of the ImThera cuff using (NIS-Elements Ar software, v5.02.01, Build 1270, Nikon Instruments Inc.) as described previously (Pelot et al 2020). We modeled the threedimensional geometry of the ImThera cuff around the modeled nerve, a 100 µm layer of saline around the cuff, and the surrounding tissue as muscle (25 mm length, 10 mm diameter). ASCENT built, meshed, and solved the finite element models using COM-SOL Multiphysics v5.4 (Burlington, MA). For each of the solved COMSOL models-one for each of the six contacts for each of the six animals (total = 6 × 6 = 36 models)-the ASCENT pipeline sampled the resulting electric potentials along the length of the nerve at the centroid of each fascicle, i.e. along the paths of axons. Using NEURON v7.6 (Hines and Carnevale 1997), ASCENT applied the potentials to biophysically-realistic models of mammalian myelinated nerve fibers, corresponding to the largediameter (13 µm MRG) myelinated Aα and mediumdiameter (3 µm MRG) myelinated B fibers (McIntyre et al 2004, Musselman et al 2021. We used a binary search algorithm (1% resolution) to identify the activation threshold for each fiber in response to a charge-balanced, biphasic pulse (200 µs/phase with no interphase delay, as used in vivo). Data plotting and analyses were done in Python v3.7. Further details are provided in supplementary materials: Computational Modeling.
We used a mixed model analysis of variance (ANOVA) for each fiber type to assess factors that influenced the models' relative error to in vivo thresholds (fixed effects: full factorial of vagotopy (included in model or not), response level (Aα fibers only), distance from the efferent mode; random effect: animal). For significant terms, we performed paired student's t-tests on the models' errors with and without the fixed effect.

Results
The goal of this study was to understand how the underlying anatomy of the pig cervical VN relates to functional response to VNS evoked using an FDA investigational device exemption (IDE)approved six-contact cuff electrode, with a focus on developing a conceptual framework for differentiating putative therapeutic effect (i.e. HR reduction) from side effect (concomitant laryngeal muscle activation). Data from the recent NECTAR trial of VNS to treat heart failure suggested that the intended changes in sympathetic/parasympathetic tone were limited by unwanted neck muscle activation (De Ferrari et al 2017); therefore, we used the stimulationevoked changes in HR as a primary measure of therapeutic engagement. In the following sections, each animal's histology is presented first to provide context for evaluating the observed functional results. Next, the physiological results are presented in context of the underlying fascicular organization and the associated positions of the active electrode contacts. Finally, subject-specific computational modeling is used to develop a conceptual framework to understand how the underlying fiber type-specific fascicular organization relates to the observed functional results.

Histology to trace efferent and afferent fascicles from nodose to surgical window
The VN has distinct functional organization at specific points along its path connecting the brainstem to the visceral organs (Ellis 1964, Câmara andGriessenauer 2015). Somatic motor efferents responsible for deep neck muscle activation originate within the nucleus ambiguus in the medulla oblongata and more distally become the pharyngeal, SL, and RL branches, which innervate the pharyngeal, CT, and CA muscles of the deep neck (figure 1). Parasympathetic efferents originate in the dorsal motor nucleus of the vagus in the medulla oblongata and distally become vagal branches to visceral organs. Visceral afferents course within these same branches back to the cervical vagus; these afferents are pseudounipolar cells with cell bodies in the nodose ganglion and central axons terminating in the solitary nucleus (NTS).
Across the cohort, histology of the VN at the level of the cuff showed a functionally clustered organization of fascicles (supplementary material: figure 4). Cranial to the cuff, pseudounipolar cell bodies aggregate in the nodose ganglion (figure 2(b), slice 1) (Settell et al 2020). These afferent neurons of the nodose remained clustered in a distinct grouping of fascicles as they descended caudally to the level of the surgical window (figure 2(b), slices 2-6, 'Afferent Mode'). We traced this organization by taking histological slices every ∼40 µm. The remaining fascicles (figure 2(b), 'Efferent Mode') contained both motor and parasympathetic efferents 14 (Ellis 1964, Rea 2014, Câmara and Griessenauer 2015, Settell et al 2020. Thus, the cross section of the nerve at the level of the stimulation cuff contained two distinct functional groupings of fascicles: an efferent grouping and an afferent grouping (figure 2(b), slice 6). Further, the efferent grouping of fascicles is expected to preferentially contain the largest diameter myelinated motor nerve fibers while the afferent fascicle cluster is expected to contain smaller diameter myelinated and unmyelinated sensory afferents from the visceral organs; indeed, we observed this distinction in fiber 14 In previous studies, we conducted immunohistochemistry with an antibody against choline acetyltransferase on the pig cervical vagus nerve. We demonstrated that cholinergic fibers consistent with both motor and parasympathetic efferents are preferentially clustered in one of the two bi-modal groupings of fascicles (Settell et al 2020). diameters between fascicles in both modes to confirm the presence or relative absence of the largest diameter myelinated fibers (figure 2(d)).

Sources of electroneurography and electromyography signals
We conducted multiple controls to confirm the sources of the various components of the EMG and ENG signals. We first considered the relevant gross anatomy, including the neural pathways activated by VNS that evoke EMG responses: the RL and SL branches of the vagus ( figure 1). Prior studies demonstrated that electrical stimulation of the cervical vagus with an epineural cuff can activate motor nerve fibers within the cervical vagus that become the RL (DuBois and Foley 1936, Foley and DuBois 1937, Nicolai et al 2020. In addition, motor nerve fibers of the SL which are located near the cuff can be activated through current leakage (Nicolai et al 2020). Across the cohort, the motor nerve fibers in the RL traverse a long path (22.89 ± 1.03 cm) from the implanted cuff electrode to their innervation of Figure 3. Representative ENG signals measured by intrafascicular electrodes. (a) Signal measured via LIFE in the VN during the 'intact' condition. Candidate eCAPs are denoted by red arrows. Candidate eCAPs in the green shaded region correspond to the latencies of Aα and B fibers, while candidate eCAPs in the red shaded region correspond to the latencies of Aγ to B fibers. In this condition, it is difficult to determine which signals are true eCAPs and which signals are EMG contamination. Also, note the stimulus artifact at the onset of the stimulation (grey shaded region). (b) Upon applying a paralytic agent, the EMG is eliminated and only the neural signal remains in the recorded trace. (c) Without paralytic, RL and SL branches are transected, eliminating direct and indirect pathways to nearby muscles within the neck, which provides an additional confirmation of the lack of EMG contamination in the ENG signal. (d) Verification of neural signal sources via double transection of the vagus nerve, rostral and caudal to the stimulation cuff, after which no signal remains. Note, the 'intact' condition (panel (a)) is the result of averaging the recordings from 750 pulses, while the remaining conditions (panels (b)-(d)) are averaged across 25 pulses, resulting in slightly more noise in the traces. the CA. Conversely, the motor nerve fibers of the SL travel a much shorter path: the external superior laryngeal branch travels 4.7 ± 0.6 cm from the nodose to the cricoid cartilage, under which it inserts in the CT muscle, and the internal superior laryngeal branch travels 3.4 ± 0.4 cm to its insertion in the CA muscle. Activation of the long RL pathway causes muscle contractions, and the resulting EMG appears as artifact in ENG recordings that could easily be conflated with activation of fibers with lower conduction velocities, i.e. long latency components in the eCAP (figure 3 and supplementary material: figure 19) (Nicolai et al 2020). Activation of the SL pathway generates a shorter latency EMG response that could also appear as artifact in the ENG recordings (supplementary material: figure 19).
We confirmed the sources of the eCAP versus EMG artifact in the ENG signal through use of a paralytic agent vecuronium and multiple nerve transections. As shown in figure 3(a), there were multiple peaks in the eCAP signal in the intact condition in the latency ranges for the Aα and B fibers. Transection of the RL and SL branches eliminated the two slowest peaks (figure 3(a), two right most red arrows), thus indicating their source was EMG contamination. However, transections of the RL and SL alone do not account for the possibility of an unanticipated motor nerve fiber pathway being activated and evoking EMG that appears as artifact in the ENG. Therefore, eCAPs were also obtained following intravenous administration of the muscle blocker vecuronium to eliminate all other confounding muscle signals. eCAP components after RL and SL transection or under a general muscle block remained unchanged (figures 3(b) vs. (c)), indicating the RL and SL transections successfully eliminated all EMG contamination from the ENG recordings. Finally, we confirmed that remaining eCAP components were truly neural signals by transecting the vagus on both sides of the stimulating cuff, which eliminated all signal ( figure 3(d)).
In addition to controls with a paralytic agent and with transections, we also observed correlations between ENG activity and physiological responses by comparing the DRCs for Aα fibers with EMG across contact locations (figure 4). Since Aα fibers have a short latency eCAP response, Aα eCAP components could be reliably discriminated from the stimulation artifact in only four of the six animals. Activation of Aα motor nerve fibers within the cervical vagus directly cause contractions of the CA muscle via the RL (figure 1). As expected, there was strong agreement between Aα ENG and the EMG DRCs with respect to the onset threshold (i.e. stimulation amplitude producing the initial rise in signal), saturation threshold (i.e. minimum stimulation amplitude producing maximum signal), and differences in the responses across the six contacts (figure 4, supplementary material: figure 22).

Location-specific differences for evoking unwanted neck muscle activation
We quantified the effects of contact location around the nerve circumference on the recruitment of Aα fibers and the associated EMG responses, leveraging our mapping of the underlying nerve morphology. Figure 4(b) shows the Aα and EMG DRCs for two representative animals as a function of contact location; the data across the contacts are color-coded by the distance of the contact from the centroid of the grouping of fascicles containing somatic motor and parasympathetic efferent fibers. The DRCs have lower thresholds for the three contacts closest to the efferent mode and higher thresholds for the three contacts furthest from the efferent mode.
We rotated the histology for each animal such that the efferent and afferent groupings were oriented toward the top and bottom of the page, respectively, to better visualize and quantify the impact of distance from the stimulating contact to the efferent fascicles on the EMG responses across the cohort (supplementary material: figure 4). Next, we calculated the current amplitude required to elicit half-maximal Aα and EMG responses (ED 50 ) from the DRC for each contact for each animal. We mapped these amplitudes to the location of the stimulating contacts and interpolated ED 50 values around the nerve circumference ( figure 4(a)). We repeated this analysis for the data for other animals (4/6 for Aα, 5/6 for EMG); the Aα signal was contaminated by the stimulation artifact in two animals, and the EMG DRC had insufficient resolution to identify clearly the points of threshold and saturation in one animal. We calculated the mean ED 50 across animals for each location around the nerve circumference (figure 4(c)). As with the individual animal data ( figure 4(a)), this aggregate analysis shows consistency between electrode distance from the efferent cluster of fascicles and the ED 50 values for the Aα and EMG responses, with lower thresholds for contact locations closer to the efferent fibers. We conducted a paired t-test comparing the Aα ED 50 values for the three contacts closest to the efferent cluster versus the ED 50 values for the three contacts furthest. As expected, the required ED 50 for the three contacts closest to the efferent cluster (0.150 ± 0.028 mA) were lower than the ED 50 values for the three contacts furthest from the efferent cluster (0.245 ± 0.083 mA; t(5) = − 3.52, p = 0.005, paired t-test). Additionally, the Euclidian distance from the stimulating contacts to the evoked ED 50 was plotted (supplementary material, figure 24).
In two of the six animals, the EMG signal increased further in response to the highest stimulation amplitudes after initially reaching a point of apparent saturation (figure 5, supplementary material: Off-Target Activation). This additional increase in EMG response reflects the appearance of a new short-latency EMG response, which was eliminated by transection of the SL branch in proximity to the cuff (figure 1, supplementary material: figure 21). These transection data confirmed that the source of this additional EMG increase was current escaping the cuff to activate the nearby and shorter SL path innervating both the CA and CT muscles, similar to previous studies (Nicolai et al 2020). Off-target activation of the SL fibers occurred at lower stimulation amplitudes for contact locations closest to the SL, which branches off at a ventral-medial angle to the main trunk (figure 5). When the ED 50 values for these animals were calculated, a truncated DRC, which only included the first plateau, was used for the analysis.

Location-specific differences for evoking changes in HR
As described above, we observed statistically significant differences in EMG responses between contact locations closest to and furthest away from the efferent mode. However, the ratio of ED 50 from the furthest to closest contact was, on average, only approximately 60% (figure 4; animal #2: and animal #4 (right) at the level of the stimulating cuff, rotated such that the efferent cluster is oriented toward the top of the page. Stimulating contacts are superimposed to replicate the in situ condition and are characterized in color by their distance to the centroid of the efferent cluster ('Efferent Distance' colorbar). A heatmap was interpolated around the perimeter of the histology to indicate the current amplitude required to elicit half-maximal Aα responses ('ED50' colorbar). (b) ENG and EMG dose-response curves for animal #2 (left) and animal #4 (right); 95% CIs are plotted but not visible due to its extremely narrow range 15 . The color of each trace corresponds to the colors of the contacts in panel (a). (c) Extrapolated ED50 heatmap for both neural (left; n = 4) and muscle (right; n = 5) responses, averaged across the cohort and superimposed on a generalized circular illustration of the nerve cross section. The smaller ellipsoid (blue) represents the efferent mode. Same colorbar for ED50 as in panel (a). 0.11-0.41 mA, animal #4: 0.14-0.29 mA). In contrast, the HR responses to VNS were more sensitive to stimulation location (figure 6). Across all animals, there was at least one stimulation contact location that evoked profound bradycardia (>10 BPM) (supplementary material: figure 23). The mean stimulation amplitude required to reduce HR by more than 3 BPM from baseline across all contacts and all animals was 1.3 ± 0.5 mA. The mean maximum decrease in HR (HR max↓ ) across all animals was 34 ± 21 BPM (range: 11-71 BPM; supplementary material: figure 23, animal #3 (smallest HR max↓ ), animal #5 (largest HR max↓ )).
Different contacts could yield broad variations in HR responses in the same animal (supplementary material: figure 23). For example, in animal #5, contact 1 yielded a HR max↓ of 17 BPM at 5 mA (21% change from baseline), whereas contact 4 yielded a HR max↓ of 71 BPM at 4 mA (89% change from baseline). In this animal, we did not stimulate at a current above 4 mA through contact 4 given the magnitude of the bradycardic response. Conversely, other contacts did not produce such a dramatic response, even at higher stimulation amplitudes. Across all animals, the difference of HR max↓ between the contact locations with the smallest and largest HR max↓ values was 35 ± 27 BPM, indicating a profound sensitivity in HR response to electrode location.
Animal #1 uniquely had profound tachycardia or marked bradycardia depending on contact location ( figure 6(b)). The histology showed that the two electrodes eliciting a tachycardic response were adjacent to a tertiary small grouping of fascicles, distinct from the established afferent and efferent clusters seen in the rest of the cohort ( figure 6(a)). This tertiary mode was assumed to be either the aortic depressor nerve or hitchhiking fascicles from the main sympathetic trunk, both of which originate from the pseudounipolar cells of the superior cervical ganglion (SCG); the aortic depressor nerve can be embedded in the cervical vagus in as many as half of pigs (Settell et al 2020).
Across the cohort, stimulation location influenced the magnitude of HR responses. Multiple fiber types can drive HR responses. Nerve fibers responsible for cardiac changes are putatively afferent Aδ fibers from mechanoreceptors, parasympathetic efferent B fibers, and visceral afferent C fibers (Kardon et al 1975, Ford and McWilliam 1986, Jones et al 1995. Consistent with studies suggesting that the HR response to VNS is primarily driven by B fibers (Sabbah et al 2011, McAllen et al 2018, Qing et al 2018, our recorded HR responses as a function of electrode contact location show that electrodes closest to the efferent cluster of fascicles containing both motor and parasympathetic efferents elicited the most profound bradycardia. To visualize and quantify this effect, we rotated the histology to place the efferent fascicle grouping oriented toward the top of the page, as was done for the previous Aα and EMG analyses. Average maximal HR decreases were calculated and fit to a circular nerve model ( figure 6(d)). Visual inspection of this aggregate figure suggests that the electrode closest to the efferent cluster also evoked the most prominent HR response. HR changes across the cohort were compared between the three contacts closest to the efferent cluster and the three contacts furthest from the efferent cluster to test this hypothesis. As expected, the maximum HR decrease of the contacts closest to the efferent cluster (21 ± 10 BPM) was significantly larger than the HR decrease of the contacts furthest from the efferent cluster (9 ± 8 BPM) (t(5) = 4.46, p = 0.0005, paired t-test). Additionally, the Euclidian distance from the stimulating contacts to the evoked HR response was plotted (supplementary material, figure 24).  (c)). (b) Heart rate (HR) dose-response curves. Colors match the stimulating contacts in panel a. Note the marked tachycardia (increased HR) in animal #1, particularly for contacts closest to the additional small cluster of fascicles at the bottom right of the histology. Conversely, both animals show bradycardia (decreased HR) for contacts near the efferent cluster at the top of each micrograph. (c) Colorbar legends for the distance of each contact to the efferent cluster and for the maximum change in HR. Note that these legends are not related to each other (e.g. 2.5 mm distance does not correspond to 0 bpm change in HR). (d) Extrapolated heatmap of the change in heart rate, averaged across the cohort (n = 5) and superimposed on a generalized circular illustration of the nerve cross section.

Computational modeling 3.5.1. Quantitative match of modeled and in vivo thresholds
We implemented models based on the specific nerve morphology and electrode locations from each pig and compared the model activation thresholds for specific fiber types to the corresponding in vivo data. We found strong agreement in onset and saturation thresholds for modeled Aα fibers (13 µm MRG) in the fascicles of the efferent mode versus Aα LIFE responses and for modeled B fiber thresholds (3 µm MRG) in the fascicles of the efferent mode versus in vivo bradycardia thresholds; including vagotopy in the models improved the match. Current amplitudes sufficient to activate model B fibers (i.e. ∼1.5 mA) also robustly activated Aα fibers across all contact locations both in the models and in vivo (figures 7(b) and (c); supplementary material: figures 8-18, even-numbered figures, panels (b) and (c)).
Across the cohort, models demonstrate generalized trends in terms of fiber activation in response to location-specific stimulation. Models predicted thresholds for Aα fibers, with stimulation amplitudes accurate to the same order of magnitude as the threshold resolution of the in vivo LIFE recordings. Our models accurately predicted B fiber onset thresholds as measured in vivo by onset of bradycardia. Models underestimated B fiber saturation thresholds (i.e. HR continued to decrease beyond the highest model threshold values) and did not consistently capture the contact-specific responses In the middle row, we only modeled fibers in the black fascicles (i.e. the efferent mode). (b) and (c) Fascicles labeled as suprathreshold (blue) and subthreshold (gray) in response to a 1.5 mA amplitude pulse for Aα and B fibers. At this stimulation amplitude, Aα fibers are activated in all fascicles, while the B fibers are activated concomitantly only in the fascicles closest to the contact. See data for all animals in supplementary material: Computational modeling data for additional animals: figures 7-18. seen in vivo for the contacts furthest from the efferent mode. It should be noted that models did not include activation of Aδ or C fibers arising from baroreceptors in the aortic arch within the vagus or hitchhiking fibers from the sympathetic trunk that may also influence experimental HR changes and may explain differences between in-vivo HR responses and modeled activation of parasympathetic B fibers innervating the heart.

Vagotopic implications of thresholds across electrode contacts
We plotted the model DRCs with and without nerve vagotopy to assess whether the afferent/efferent fascicular organization determined by histology aligned with the electrode-specific ENG, EMG, and HR responses (example data provided in figure 7; data for all animals provided in supplementary material: Computational Modeling Data for Additional Animals, figures 7-18).
Without nerve vagotopy (i.e. each fiber type placed in every fascicle), we observed only small differences in the DRCs across contact locations for Aα and B fibers (figure 7(a): top row). Models incorporating vagotopy, with Aα and B fibers only in the fascicles of the efferent mode, improved the matches between models and experimental results (figure 7(a): middle row) and yielded more distinct contact-specific responses. Together with the experimentally documented location-specific differences in physiological responses, the model data confirm that contact location-specific responses to stimulation can be attributed to the spatially distinct locations of efferent fibers.

Discussion
Multiple cervical VNS clinical studies report intolerable, therapy-limiting side effects resulting from neck muscle activation (The Vagus Nerve Stimulation Study Group 1995, Handforth et al 1998, De Ferrari et al 2017. Biophysical principles dictate that largediameter nerve fibers activate at lower stimulation amplitudes than small-diameter nerve fibers (Tasaki 1953). Thus, neck muscle activation is an expected limitation since smaller diameter fibers putatively cause therapeutic effects while large diameter A fibers cause muscle contractions. However, both therapeutic effects and side effects do not have a binary on/off switch, but manifest over a dynamic range due to biophysical complexities.
The results presented here demonstrate clear location-specific differences in threshold and saturation for EMG and HR responses. Although a multicontact stimulation device such as the ImThera has the potential to shift both these curves in favorable directions, the shift does not overcome the fact that motor efferent nerves and therefore unwanted side effects occur prior to changes in HR. Future strategies such as habituation, desensitization, and alternative stimulation waveforms may be possible to further shift these curves and increase the effectiveness of treatment to side effect ratio, within certain limits. Additionally, other multi-contact cuff designs and stimulation paradigms which take advantage of the increased length of the VN, as opposed to the hypoglossal which the ImThera was designed for, could generate similar or preferred results (Ohsawa and Inui 2009, Dali et al 2018, Shulgach et al 2021. As the ImThera device has previously received an FDA IDE approval for human studies of stimulation of the hypoglossal nerve, the path to evaluate these results in human studies is likely to be dramatically accelerated in comparison to completely new electrode designs.

ImThera cuff design and safety considerations
The ImThera multi-contact cuff was approved for sale in Europe (CE Marked) to treat sleep apnea in 2021 and received an FDA IDE to conduct a U.S. Pivotal Trial in 2014. The electrode configuration was designed to stimulate the motor efferents of the human hypoglossal nerve innervating the muscles of the tongue to treat obstructive sleep apnea. As motor nerve fibers can be activated at low current amplitudes, the electrodes did not need to be designed to support higher current amplitudes needed to activate the Aδ, B, and C fibers putatively associated with the therapeutic effects of VNS. The 2 mm diameter disk electrodes have a smaller surface area than the existing LivaNova bipolar helical cuff leads, and therefore the charge density applied at the same current amplitude is larger, potentially resulting in deleterious electrochemical reactions.
We followed the protocol outlined by Wilks/ Ludwig et al (2017) that was previously used to support FDA IDEs to assess the maximum safe limit for stimulation from these electrodes (supplementary material: figure 1). VT analysis suggested that the average safe limitation was 5.6 ± 0.34 mA. Consequently, the maximum current applied in this acute study was 5 mA. Given that thresholds tend to increase in chronic implantations due to the formation of a scar around the electrodes, it may be necessary to pursue electrode coatings strategies to increase the electrochemical surface area of the electrodes to activate B and C fibers safely in human subjects, especially given subject-to-subject variability in the human anatomy.
It is also worth noting that each of the six ImThera contacts (contact diameter ∼2 mm) spans approximately one fifth of the circumference of the pig and human cervical vagus (Hammer et al 2018, Pelot et al 2020. This large stimulating surface area results in overlapping coverage of the underlying neuroanatomy. Thus, it may be necessary to pursue smaller electrodes to optimize for spatial selectivity based on electrode location. However, as noted above, this would further increase charge density at a given amplitude, potentially necessitating coating strategies to safely activate higher threshold therapeutic fibers.

Implications of pertinent anatomy 4.2.1. Current escape activates the SL branch
In two of six animals, a secondary plateau was noted in the EMG DRCs as stimulation amplitudes were increased. This phenomenon was unlikely influenced by time or repeated stimuli because stimulation amplitudes were randomized. Additionally, previous VNS work in pigs using the LivaNova helical cuff showed that high stimulation currents can indirectly activate the nearby SL branch via current escape from the insulation, leading to additional activation neck muscles in the majority of animals (Nicolai et al 2020). The Nicolai study of the LivaNova lead used bipolar stimulation, which generates a much sharper falloff of the electrode field compared to the monopolar stimulation used in this study (1/r 2 vs. 1/r). Nonetheless, the Nicolai data showed off-target activation of the nearby SL through current leakage from the insulation in 9 of 12 pigs at current amplitudes ranging from 0.5 to 2.5 mA.
In contrast, the ImThera cuff may be better suited for preventing current escape than a helical device in part due to a greater amount of insulation: both in terms of length from stimulating contact to the edge of insulation and thickness of insulation. This may explain the relatively limited number of subjects (2/6 herein vs. 9/12 in Nicolai) in which indirect activation was observed compared to studies using the helical design (Nicolai et al 2020). The relative proximity of the specific contacts that resulted in offtarget activation of the SL played a role within both these subjects; previous studies showed that the main branch point of the SL off the cervical vagus trunk is much closer to the accessible surgical window in pigs than in humans (Settell et al 2020). However, the external branch of the SL still lies relatively close to the most common electrode placement in humans, and therefore somatic nerve activation via field escape from the insulation is likely still an issue in human subjects (Kierner et al 1998).

Nerve vagotopy dictates function
The ideal VNS intervention is one where therapyinducing effects are maximized while concomitant therapy-impeding side effects are minimized. In this study, we observed varied location of the efferent cluster within the VN, as well as a varied relational positioning of stimulating contacts to the efferent cluster (supplementary material: figure 4). Yet, we observed neck muscles and HR responses dependent on their proximity to underlying neuroanatomy. However, it is essential to note that at no contact location in any animal were thresholds for desired cardiac responses lower than the thresholds for unwanted neck muscle activation via the RL. It was, however, possible to elicit cardiac responses without causing activation of the SL through current leakage outside the insulation in four out of six animals. Further, in one of the animals experiencing muscle activation via current leakage, it was still possible to stimulate at high amplitudes in three of six contacts distant from the SL without causing additional muscle activation.
Although EMG saturation thresholds in this study were lower than thresholds for cardiac effects (figure 8), it may be possible to habituate and desensitize the patient to these side effects (Penry andDean 1990, Ardell et al 2017). Clinical guidance recommends a two week healing period, followed by a five to eight week titration period in which patients are habituated to therapeutic doses (Fisher et al 2021, LivaNova 2021. While the exact mechanism of habitation is unclear, the prevalence of sideeffects has been shown to diminish over time (Morris and Mueller 1999, Ben-Menachem 2001, Ardell et al 2017. Additionally, it may be possible to alter stimulation paradigms to favorably target a specific type of nerve fiber (Vuckovic et al 2008, Dali et al 2016, Fitchett et al 2021, and/or to use location in conjunction with more fiber-selective stimulation strategies. Habituation, desensitization, and alternative stimulation paradigms may provide enough of a shift to further increase the effectiveness of VNS due to locationspecific stimulation. The most striking example of contact-specific stimulation yielding varying results was evoking both bradycardia and tachycardia in animal #1. In this animal, the two stimulating contacts closest to a hitchhiking, third grouping of fascicles traced back to the SCG of the sympathetic trunk resulted in a HR increase of over 50 BPM, while three other contacts resulted in an average decrease of 13 BPM (figure 6, animal #1). Cross-connections between the VN and the sympathetic trunk (Settell et al 2020) and sympathetic fibers within the VN (Seki et al 2014) have previously been reported. These hitchhiking fibers may be an important consideration for understanding the underlying neuroanatomy responsible for mediating competing increases and decreases in sympathetic tone.
Prior studies using the LivaNova helical cuff in canines recorded tachycardia at low current amplitudes putatively due to the activation of afferents which result in parasympathetic inhibition, with bradycardia only being induced at higher levels of stimulation (Ardell et al 2017). This change between tachycardia at low levels of stimulation and bradycardia at high levels of stimulation is known as the 'neural fulcrum.' Tachycardia was not observed at low levels of stimulation in most animals in this study. However, isoflurane was used in this study which is known to act on glutamate receptors. Although blocking glutamate-mediated responses would not blunt cholinergic parasympathetic efferent to the heart, it may blunt the afferent pathway driving reflexive parasympathetic inhibition. Another possibility is that the relatively simple cervical vagus fascicular organization in canines consisting of 1-3 fascicles (Onkka et al 2013, Yoo et al 2013 has significantly different functional responses to VNS than the more complex fascicular organization of the pig vagus (Pelot et al 2020, Settell et al 2020.

Computational modeling
Our modeling results support the in vivo findings: the ImThera cuff in monopolar configuration provides no opportunity to activate B fibers in the pig VN before robust concomitant activation of side effect producing Aα fibers, however, evoked threshold and ED 50 are responsive to location-based stimulation. Furthermore, our paired experimental and modeling data highlight the importance of accounting for individual-specific nerve morphology, nerve vagotopy, and distinct fiber types when attempting to target specific physiological responses.
Although previous studies modeled stimulation of individual-specific nerves to target selectively groups of fascicles that innervate distinct muscles (Schiefer et al 2008, Dali et al 2018 and others demonstrated stimulation of functionally distinct fascicles of the VN in vivo (Ordelman et al 2013, we are the first to quantify selective stimulation with individual-specific models that incorporate nerve vagotopy. The accuracy of the model was enabled by incorporating controlled cuff orientation, nerve vagotopy and accurate conduction velocities for the fibers evoking distinct physiological responses, and models that correspond oneto-one with nerves stimulated in vivo achieved the strongest contact-specific match for Aα and B fiber thresholds.
Across animals, contact locations, and fiber types, our models with vagotopy slightly underestimated in vivo thresholds at all response levels, indicating systematic error either in our models or in vivo signal classification. Models incorporating nerve vagotopy accurately excluded fascicles from our calculation of modeled response that did not contain the fibers that contributed to the recorded in vivo signal, and this is reflected by the differences we observed between DRCs for models with and without vagotopy.
Models enable rapid, reproducible, and efficient quantification of therapy design choices on nerve response to stimulation. Specifically for VNS, patients need reduced side effects, which will require avoiding activation of side effect producing fibers both within the nerve trunk proper and in nearby branches (e.g. the SL) (Nicolai et al 2020).

Limitations
Cardiac effects of isoflurane and fentanyl are well characterized, but their competing interaction is complex. Isoflurane can lead to vagal inhibition, causing an increase in HR (Kato et al 1992, Huang et al 1997, Picker et al 2001. In contrast, fentanyl decreases HR due to an increase in vagal tone (Laubie et al 1977, Griffioen et al 2004. All recordings analyzed during this study were obtained during periods when isoflurane concentration was less than 2%. Efforts were made to keep the isoflurane concentrations minimal during the entire procedure. However, due to animal safety and welfare concerns, isoflurane concentrations occasionally rose above 2%. Synaptic blunting caused by isoflurane anesthesia may require activation of a larger number of efferent parasympathetic fibers or afferent fibers from mechanoreceptors to elicit a HR response, and therefore the anesthesia may have affected the activation and saturation thresholds, as well as the total change in HR (Herring et al 2009, Baumgart et al 2015.
Putative nerve fibers responsible for cardiac changes are afferent Aδ fibers from mechanoreceptors, parasympathetic efferent B fibers, and visceral afferent C fibers (Kardon et al 1975, Ford and McWilliam 1986, Jones et al 1995. However, we were unable to obtain LIFE recordings from Aδ, B, or C fibers. Prior studies suggest that C fibers are not activated within clinically tolerated stimulation intensities (Yoo et al 2013) and that Aδ fibers may only play a small role in the cardiac response, especially in the right VN (Ito and Scher 1973). Therefore, B fiber measurements are the best predictors for cardiac activity (Sabbah et al 2011, Qing et al 2018. However, we were unable to record B fiber activity and analyze this link in the pathway between applied stimulus and cardiac response. This is unsurprising, as other studies have either only intermittently captured recordings from these fiber types using LIFE electrodes (Nicolai et al 2020, Verma et al 2022, or demonstrated that measuring changes in potential due to these fibers may require surgically removing the epineurium of the VN in large animal models (Yoo et al 2013). Instead, we recorded changes in HR as therapeutic effect and recorded EMG activity as side effects, analogous to clinical neck pain, throat tightening, cough, dyspnea, or hoarseness (The Vagus Nerve Stimulation Study Group 1995, Handforth et al 1998, De Ferrari et al 2017, although these may be oversimplifications. We may have therefore come to ambitiously broad conclusions, albeit having sound inferential logic, meriting further validation in future studies.
Acute surgical trauma often resulted in fluid and edema at or near the electrode/nerve interface. Such fluid was wicked away via surgical gauze as needed. This fluid buildup may have impacted certain variables such as the spatial distribution of applied currents and thus stimulation thresholds. In addition, the surgery may have impacted additional biophysiological phenomena, such as indirect activation of nearby anatomical structures. Therefore, further studies with chronically implanted electrodes are needed to understand the roles that edema, fluid, scar tissue formation, and biofouling play in the responses described in this paper.
The computational models have limitations to their complexity for maintaining feasibility. Specifically, we assumed a constant nerve cross section along the model length, a radially uniform correction to account for tissue shrinkage during histology, and a physics-based method for reshaping the nerve to a circular cross section in the cuff electrode and repositioning the fascicles accordingly. Additionally, factors intrinsic to the nerve such as conduction differences between fat and epineurium, and extrinsic factors such as edema, may impact the accuracy to which the models simulate the in vivo environment. Our models may have overcorrected for tissue shrinkage, as shown by observed differences in experiment and modelbased determination of contact location after the nerve was deformed to the inner diameter of the cuff. The models maintained the distances between contacts as the cuff was expanded to accommodate the nerve; therefore, accurate contact placement relative to the fascicles inside the nerve required that we have an accurate parameter for correcting tissue shrinkage from histological processing, which has exhibited a large range (Hursh 1939, Friede and Samorajski 1967, Stickland 1975, Boyd and Kalu 1979. Further, we did not model distributions of diameters of each fiber type, and some of the largest B fibers could be recruited before the smallest Aα fibers. There were limitations in the validation of model data against corresponding in vivo measurements. We were limited in the resolution of amplitudes tested in vivo; therefore, for each onset and saturation threshold in the in vivo data, there was a range of range within which the response started or saturated. Further, the lowest and highest thresholds across fascicles are known in the models, but we do not know precisely which fascicles contributed to the associated ENG/EMG signal in vivo. We assumed that the Aα and B fibers were both uniformly distributed across the fascicles of the efferent mode, but the specific fibers that evoked the ENG, EMG, and HR responses each constitute an unknown subset of those fiber locations (Kronsteiner et al 2022).

Conclusion
Therapeutic outcomes of VNS are limited in part by the stimulation-induced side effect of neck muscle activation. In this study, we sought to characterize the fundamental relationship between the underlying anatomical organization of fiber types within the cervical vagus trunk in pigs and the functional outcomes of location-specific electrical stimulation in terms of both effect and side effect. Functional and histological data were used to develop and validate computational models of location-specific activation. Both in vivo and in silico results demonstrate location-specific differences in threshold and saturation for responses related to effect and side effect. While we were unable to evoke intended effect prior to evoking side effect, we were able to demonstrate location-specific stimulation can avoid fibers responsible for evoking tachycardia, can eliminate current escape activating nearby muscle groups, and can activate fibers responsible for intended effect at lower amplitudes. Therefore, a multi-contact stimulation device such as the ImThera cuff may exploit the underlying fascicular bimodality of efferent versus afferent fiber types and while not mitigating side effect, may shift the ratio of therapeutic effect to noxious side effect in a favorable direction. This shifting of effect to side effect ratios in addition to common clinical practices such as desensitization and habituation, may potentially increasing the effectiveness of treatment.
Previous studies determined the importance of certain VNS parameters such as stimulation amplitude and frequency. These data build upon other studies and stress an equally important parameter: stimulation location.
Baumgart for their help in collecting electrophysiological data, Maria LaLuzerne for providing logistical and experimental support, Jennifer Frank and Keri Graff for providing technical veterinary support, and to Daniel Marshall for segmenting the nerve histology used to define nerve cross sections in ASCENT. A sincere thank you to Nishant Verma and Dr Wolf-Ekkehard Blanz for proofreading the article.
The opinions expressed in this article are the author's own and do not reflect the view of the National Institutes of Health, the Department of Health and Human Services, or the United States government.

Author contributions
J W and K L are scientific board members and have stock interests in NeuroOne Medical Inc., a company developing next-generation epilepsy monitoring devices. J W also has an equity interest in NeuroNexus technology Inc., a company that supplies electrophysiology equipment and multichannel probes to the neuroscience research community. K L is also a paid member of the scientific advisory board of Cala Health, Blackfynn, Abbott, and Battelle. K L also is a paid consultant for Galvani and Boston Scientific. K L and A S are consultants to and co-founders of Neuronoff Inc. None of these associations are directly relevant to the work presented in this manuscript.
Additionally, R V and J B are employed by LivaNova USA Inc., a vagus nerve stimulation company, and hold stock or stock options. R V and J B contributed primarily to the initial conceptualization and financial support for the program. The authorship team asserts that their conflict did not influence the collection, analysis, or interpretation of study results. The authorship team asserts that financial contributions of LivaNova USA funded study-related expenses but were not allocated for authorship.

Conflict of interest
The remaining authors declare that the research was conducted in the absence of any additional commercial or financial relationships that could be construed as a potential conflict of interest.