Field models: pig and human vagus nerve stimulationWe implemented and simulated finite element models of pig and human vagus nerves with cuff electrodes using ASCENT v1.2.1 (DOI: 10.5281/zenodo.7627427)55, with COMSOL Multiphysics® v5.6 (COMSOL, Inc.; Burlington, MA, USA).We modeled realistic morphologies of pig (n = 7) and human (n = 7) cervical vagus nerves using segmented histology of nerve cross sections56,57,58 that we extruded 50 mm longitudinally (Supplementary Note 1). We made each nerve circular while preserving its cross-sectional area using ASCENT’s deformation feature; fascicles were repositioned during deformation to maintain at least 10 µm between fascicles and between fascicles and the nerve boundary.We modeled two cuff electrodes (Fig. 1). First, we modeled the six-contact ImThera cuff electrode (LivaNova PLC, London, UK). The cuff was 9 mm in length, with an inner diameter of 3 mm and six circular contacts (2 mm diameter each) arranged along two diagonals (Fig. 1a). The ImThera cuff geometry and validated models of pig vagus nerve stimulation (VNS) with the ImThera cuff are published in prior papers35,59. We also modeled the bipolar helical cuff used to treat epilepsy and depression clinically (LivaNova PLC, London, UK) (Fig. 1b), which is available in two sizes. For nerves with a diameter under 3 mm (n = 11), we used the 2 mm diameter cuff; for larger nerves, we used the 3 mm diameter cuff (n = 3). In its unexpanded state, the helical metal contacts of the cuffs spanned 338.5 degrees for the 2-mm cuff and 346.5 degrees for the 3-mm cuff. The contacts were positioned 8 mm apart, measured from their centers. Further details on the helical cuff’s design and validated models of human VNS using this cuff are previously published60. We used the ImThera cuff to generate training data (see “Surrogate model: Training”), and subsequent tests were performed using both the ImThera and helical cuffs.For both cuff geometries, we positioned the cuff at the midpoint of the nerve’s length. Spaces between the cuff and nerve were filled with saline. We modeled a 10 µm saline gap between the ImThera cuff and the nerve and a 100 µm saline gap between the helical cuff and the nerve. We placed the nerve and cuff in a cylinder of muscle (5 mm diameter, 50 mm length). We rotated and shifted each cuff based on the size and position of each nerve’s fascicles. Specifically, we determined the centroid of all fascicles (“fascicle centroid”) by applying a weighted average (using each fascicle’s area) to all individual fascicle centroids, and then found the closest point on the nerve boundary to that fascicle centroid. We then rotated the center of the cuff to align with the closest point. The cuff was then shifted towards the nerve to leave a 10 µm (ImThera) or 100 µm (helical) gap between the cuff and nerve.We assigned materials to the geometry as shown in Fig. 1; see Supplementary Note 2 for material and tissue conductivities. Some human nerves feature “peanut fascicles”, where multiple endoneurial bundles (inner perineurium boundaries) are within a single outer perineurium boundary, separated by a perineurium septum (e.g., Supplementary Note 1 H3). We meshed the perineurium for all peanut fascicles, otherwise we modeled the perineurium using a thin layer approximation:$${\rho }_{{{\rm{surface}}}}=\frac{{{{\rm{thk}}}}_{{{\rm{peri}}}}}{{\sigma }_{{{\rm{peri}}}}}$$
(1)
where ρsurface is the contact impedance of the perineurium boundary (Ωm2), σperi is the conductivity of the perineurium (S m−1), and thkperi is the thickness of the perineurium (m). For human nerve morphologies, we calculated the perineurium thickness for each fascicle based on the original segmentations:$${{{\rm{thk}}}}_{{{\rm{peri}}}}={r}_{{{\rm{outer}}}}-{r}_{{{\rm{inner}}}}$$
(2)
where router, inner are the effective circular radii of the outer and inner perineurium boundaries, respectively. For pig morphologies, we applied a published relationship between fascicle diameter, dfasc (μm), and perineurium thickness for the pig vagus nerve56:$${{{\rm{thk}}}}_{{{\rm{peri}}}}=3.44\,{{\rm{\mu }}}{{\rm{m}}}+0.02547\times {d}_{{{\rm{fasc}}}}$$
(3)
We solved Laplace’s equation:$$\nabla \cdot \left(\sigma \cdot \nabla V\right)=0$$
(4)
using quadratic geometry and solution shape functions once for each electrode contact in each model. When solving for each contact, we applied a 1 mA point source of current in the center of the electrode contact, a floating potential boundary condition for all other electrode contacts, and grounded outer boundaries of the model. Pig finite element models (FEMs) had 35 to 95 million free tetrahedral elements, and human FEMs had 2 to 35 million free tetrahedral elements.Surrogate modelDesign and implementationThe surrogate model (S-MF, Fig. 1g) is a simplified version of the McIntyre-Richardson-Grill (MRG) model of a mammalian myelinated fiber12 (Fig. 1f) that calculates the transmembrane potential (Vm) and all gating variables (m, h, p, s) at every node of Ranvier and time point. S-MF retains only the nodes of Ranvier and the intracellular space connecting nodes, thus assuming perfectly insulating myelin and omitting additional extracellular field layers. S-MF, like the original MRG model, assumes that the fiber is in a large, highly conductive medium (i.e., equivalent circuit nodes representing extracellular voltage are grounded and have nonzero potential only when an extracellular stimulus is applied). The inputs to S-MF are a 2D array defining the extracellular electric field in space (along the length of the fiber, at each node of Ranvier) and in time (according to the chosen stimulation waveform) (Fig. 1d), the diameter of the fiber, and optionally, any intracellular current sources.We solved the system of differential equations of the surrogate model using forward Euler discretization and a staggered timestep. To advance the simulation by one timestep, we first update the gating parameter values from time t − 0.5dt to t + 0.5dt using the analytic solution:$${x}^{n,t+0.5{{\rm{dt}}}}={x}_{{{\infty }}}\left({V}_{{{\rm{m}}}}^{n,t}\right)-\left({x}_{{{\infty }}}\left({V}_{{{\rm{m}}}}^{n,t}\right)-{x}^{n,t-0.5{{\rm{dt}}}}\right){e}^{-{{\rm{dt}}}/{\tau }_{x}({V}_{{{\rm{m}}}}^{n,t})}$$
(5)
$${{x}_{{{\infty }}}}(V_{{{\rm{m}}}})=\frac{{\alpha }_{x}({V}_{{{\rm{m}}}})}{{\alpha }_{x}\left({V}_{{{\rm{m}}}}\right)+{\beta }_{x}({V}_{{{\rm{m}}}})}$$
(6)
$${{\tau }_{x}}(V_{{{\rm{m}}}})=\frac{1}{{q}_{10}^{x}\times \left({\alpha }_{x}\left({V}_{{{\rm{m}}}}\right)+{\beta }_{x}({V}_{{{\rm{m}}}})\right)}$$
(7)
$${q}_{10}^{x}={{{\rm{a}}}{{{\rm{q}}}}_{10}^{x}}^{\frac{{T}_{c}-{{{\rm{bq}}}}_{10}^{x}}{{{{\rm{cq}}}}_{10}^{x}}}$$
(8)
where n is the index of the node of Ranvier, \({V}_{{{\rm{m}}}}^{n,t}\) is the transmembrane potential for node n at time t, and αx and βx are independent rate constant functions, as defined in the MRG paper12, for gating variable x (where x is m, h, p, or s). \({q}_{10}^{x}\) is the q10 value for gating variable x, parameterized by scalars \({{{\rm{aq}}}}_{10}^{x},\,{{{\rm{bq}}}}_{10}^{x},\,{{{\rm{cq}}}}_{10}^{x}\), where Tc = 37 °C is the temperature. Using m to indicate \({m}^{n,t+0.5{dt}}\) (likewise for all other gating variables h, p, and s), we then calculate the ionic current and advance the transmembrane potential from time t to t + dt:$${I}_{{{\rm{ion}}}}^{n,t}={\bar{g}}_{{{\rm{Naf}}}}{m}^{3}h\left({V}_{{{\rm{m}}}}^{n,t}-{E}_{{{\rm{Na}}}}\right)+{\bar{g}}_{{{\rm{Nap}}}}{p}^{3}\left({V}_{{{\rm{m}}}}^{n,t}-{E}_{{{\rm{Na}}}}\right)+{\bar{g}}_{{{\rm{K}}}}s\left({V}_{{{\rm{m}}}}^{n,t}-{E}_{{{\rm{K}}}}\right)+{g}_{{{\rm{L}}}}({V}_{{{\rm{m}}}}^{n,t}-{E}_{{{\rm{L}}}})$$
(9)
$${V}_{{{\rm{m}}}}^{n,t+{{\rm{dt}}}}={V}_{{{\rm{m}}}}^{n,t}+{dt}\times \frac{1}{{C}_{{{\rm{n}}}}}\times \left(\frac{1}{{R}_{{{\rm{a}}}}}\times {{\rm{CNN}}}\left(\left[{V}_{{{\rm{m}}}}^{n,t},{V}_{{{\rm{e}}}}^{n,t}\right]\right)-{I}_{{{\rm{ion}}}}^{n,t}+{I}_{{{\rm{stim}}}}^{n,t}\right)$$
(10)
$${C}_{{{\rm{n}}}}={c}_{{{\rm{m}}}}\times \pi \times {d}_{{{\rm{node}}}}$$
(11)
$${d}_{{{\rm{node}}}}={\beta }_{1}^{{{\rm{n}}}}\times {D}^{2}+{\beta }_{2}^{{{\rm{n}}}}\times D+{\beta }_{3}^{{{\rm{n}}}}$$
(12)
$${R}_{{{\rm{a}}}}=\frac{{\rho }_{{{\rm{a}}}}{\partial }_{{{\rm{x}}}}}{\pi \times {\left(\frac{{d}_{{{\rm{axon}}}}}{2}\right)}^{2}}$$
(13)
$${\partial }_{{{\rm{x}}}}={\beta }_{1}^{{{\rm{d}}}}\times {D}^{2}+{\beta }_{2}^{{{\rm{d}}}}\times D+{\beta }_{3}^{{{\rm{d}}}}$$
(14)
$${d}_{{{\rm{axon}}}}={\beta }_{1}^{{{\rm{a}}}}\times {D}^{2}+{\beta }_{2}^{{{\rm{a}}}}\times D+{\beta }_{3}^{{{\rm{a}}}}$$
(15)
where \({{{\mathbf{\beta }}}}^{{{\bf{n}}}}{{\boldsymbol{,}}}\,{{{\mathbf{\beta }}}}^{{{\bf{d}}}}{{\boldsymbol{,}}}\, {{{\mathbf{\beta }}}}^{{{\bf{a}}}}\in {{\mathbb{R}}}^{3}\) are coefficients of interpolants relating fiber diameter (D) to node diameter (dnode), internodal distance (\({\partial }_{{{\rm{x}}}}\)), and axon internodal diameter (daxon), respectively. We used second-order polynomials for these interpolants to remain consistent with previously developed interpolants for the NEURON MRG implementation55 (Supplementary Note 3.5). \({{\rm{CNN}}}\left(\left[{V}_{{{\rm{m}}}}^{n,t},{V}_{{{\rm{e}}}}^{n,t}\right]\right)\) is a symmetric linear convolution over concatenated Vm and Ve vectors with a three-node-wide kernel centered on node n at time t. \({I}_{{{\rm{stim}}}}^{n,t}\) is the intracellular current injected at node n at time t.Parameters were initialized to the values used in the MRG NEURON model (Supplementary Note 3), with CNN initialized to the second spatial difference. We used backpropagation and gradient descent to optimize 26 parameters of the surrogate model, as detailed in “Surrogate model: training”, to reproduce the responses of the NEURON model.We implemented S-MF in PyTorch v2.061 and deployed it on an NVIDIA RTX A5000 24GB GPU. All arithmetic operations calculating updates to state variables, determining the presence of action potentials, and performing backpropagation were executed on the GPU.TrainingAll training data were generated from NEURON simulations of the MRG myelinated fiber12. All simulations were conducted in NEURON v7.8 and solved using the implicit Euler method with a fixed timestep (dt) of 0.005 ms and tstop of 5 ms.The ultrastructural geometric parameters for the fiber models were interpolated over the originally published fiber geometries (Supplementary Note 3). The training dataset consisted of 65,536 field-response pairs, where the field was a 2D array representing the extracellular potential at every nodal compartment of the model neuron at each timepoint, and the response was a 3D array representing the states of the neuron (Vm, m, h, p, and s) at every nodal compartment at each timepoint. The training data were generated using two vagus nerve morphologies, one pig (Supplementary Note 1 P1) and one human (Supplementary Note 1 H1), instrumented with the ImThera cuff (Fig. 1a). We followed the procedure outlined in Algorithm 1 (Box 1), where S = 65,536 is the number of sampled field-response pairs, C = 6 is the number of contacts in the cuff used to generate the data, N = 53 is the number of nodes of Ranvier per fiber, and T = 1000 is the number of simulated timesteps.All modeled fibers had their central node of Ranvier positioned halfway along the modeled nerve with a random longitudinal offset \(\sim U[-\frac{{{\rm{INL}}}}{2},\frac{{{\rm{INL}}}}{2})\), where INL is the internodal length. NEURON simulations were distributed over 400 CPU cores on the Duke Compute Cluster, using OpenMPI62 for load balancing and parallel HDF5 for efficient data handling. The full training set was ~100 GB and took ~2 h to generate. We used 80% of the generated data for training, and 20% for validation.S-MF was trained by gradient descent on minibatches of data that were randomly sampled from the field-response pairs of the training dataset. Training was performed on an NVIDIA RTX A5000 24GB GPU and took ~5 h. We used double-precision floating-point arithmetic for training to mitigate underflow error during the initial training period; subsequent tests using the trained surrogate were performed using single-precision floating-point arithmetic for improved speed.Twenty-six parameters of S-MF were optimized during training (i.e., θ, “optimizable surrogate model parameters”): maximum ionic conductance values (\({\bar{g}}_{{{\rm{Naf}}}}\), \({\bar{g}}_{{{\rm{Nap}}}}\), \({\bar{g}}_{{{\rm{K}}}}\), and gL), axial intracellular resistivity (ρa), membrane capacitance (cm), coefficients of the dnode interpolant (\({\beta }_{1}^{{{\rm{n}}}},{\beta }_{2}^{{{\rm{n}}}},{\beta }_{3}^{{{\rm{n}}}}\)), coefficients of the daxon interpolant (\({\beta }_{1}^{{{\rm{a}}}},{\beta }_{2}^{{{\rm{a}}}},{\beta }_{3}^{{{\rm{a}}}}\)), q10 parameters (\({{{\rm{aq}}}}_{10}^{m},{{{\rm{aq}}}}_{10}^{h},{{{\rm{aq}}}}_{10}^{p},{{{\rm{aq}}}}_{10}^{s}\)), parameters governing K+ dynamics (Supplementary Note 3.3, denoted sK below), and the values of the symmetric convolutional kernel of CNN (4 values). The objective function for training was to minimize the mean squared error between S-MF and NEURON predictions:$${{{\rm{argmin}}}}_{{\bar{g}}_{{{\rm{Naf}}}},{\bar{g}}_{{{\rm{Nap}}}},{\bar{g}}_{{{\rm{K}}}},{g}_{{{\rm{L}}}},{\rho }_{{{\rm{a}}}},{c}_{{{\rm{m}}}},{{{\rm{aq}}}}_{10}^{m,h,p,s},{{{\mathbf{\beta }}}}^{{{\bf{n}}}},{{{\mathbf{\beta }}}}^{{{\bf{a}}}},{{{\bf{s}}}}^{{{\rm{K}}}}{{\boldsymbol{,}}}{{\bf{CNN}}}}{||}{{{\bf{y}}}{{\boldsymbol{-}}}\hat{{{\bf{y}}}}{{\boldsymbol{||}}}}_{{{\rm{MSE}}}}$$
(16)
where y is the NEURON model prediction, and \(\hat{{{\bf{y}}}}\) is the S-MF prediction. See Supplementary Note 4 for discussion of the significance of reparametrizing the membrane potential ODE with a trained convolutional network.For a given training epoch, let \(\left\langle {{{\bf{x}}}}_{i}\in {{\mathbb{R}}}^{B\times N\times T\times 1},{{{\bf{y}}}}_{i}\in {{\mathbb{R}}}^{B\times N\times T\times 5}\right\rangle\) be a randomly sampled (without replacement) minibatch of field-response pairs from the full training dataset, where \(i\in \{{1,2},\ldots \lceil\frac{{{\rm{training}}}\; {{\rm{set}}}\; {{\rm{size}}}={65,536}}{B}\rceil\}\) is the minibatch index, B is the size of the sampled minibatch (i.e., number of simulated fibers in the sample; B = 64), N is the number of modeled nodes of Ranvier (53 nodes), and T is the number of simulated time steps (5.0 ms/0.005 ms = 1000 time steps); the fourth dimension is 1 for the extracellular potentials of \({{{\bf{x}}}}_{i}\) and 5 for the number of state variables of \({{{\bf{y}}}}_{i}\) (\({V}_{{{\rm{m}}}}\), m, h, p, and s).Let \({{\rm{S}}}:{{\bf{x}}}\in {{\mathbb{R}}}^{B\times N\times T\times 1},{{\bf{s}}}\in {{\mathbb{R}}}^{B\times N\times 5}\to {\hat{{{\bf{y}}}}{{\boldsymbol{\in }}}{\mathbb{R}}}^{B\times N\times T\times 5}\) represent the mapping performed by our surrogate model S from input x (field) to predicted output \(\hat{{{\bf{y}}}}\) (fiber response) with initial condition (IC) s (using the numerical methods described in “Surrogate model: Design and implementation”). Let \({{\bf{SSMRG}}}\) (“steady-state MRG”) represent the state of MRG nodal compartments at rest (Vm = −80 mV and all four gating variables \(x\) at x∞(−80 mV)).Let \({{\rm{Adam}}}({{\boldsymbol{\theta }}},{\nabla }_{{{\boldsymbol{\theta }}}}L,\lambda )\) represent a single gradient descent update to optimizable surrogate model parameters \({{\boldsymbol{\theta }}}\) using the Adam update rule with default hyperparameters45 and learning rate λ, where gradients of a scalar valued loss L (mean squared error) with respect to the surrogate model parameters θ are represented by \({\nabla }_{{{\boldsymbol{\theta }}}}L\).We trained S-MF over 5 epochs using the procedure outlined in Algorithm 2 (Box 2), a variation on truncated backpropagation through time63 with the Adam update rule. Training was performed on temporally sequential disjoint chunks of minibatch data, with chunk length (in terms of number of timesteps) denoted ∆. We used minibatch size B = 64 fibers, ∆ = 50, and \(\lambda=1e-5\).Testing—activation thresholdsWe simulated activation thresholds using both S-MF and NEURON models to quantify the accuracy of S-MF across fiber diameters, stimulus waveforms, nerve morphologies, and cuff geometries. We modeled 6 human and 6 pig vagus nerve morphologies (Supplementary Note 1 P2-7, H2-7) each instrumented with the six-contact ImThera cuff and the helical cuff for a total of 24 finite element models for testing; these nerves were constructed from different histological cross sections than those used to train S-MF.For the ImThera cuff, for each nerve, for each of the six electrode contacts, we randomly (with replacement) selected an associated fascicle. We modeled 8 fiber diameters (5.7 to 14 µm) at the centroid of each selected fascicle and delivered monopolar stimulation from the associated electrode contact. For the helical cuff, for each nerve, we randomly (without replacement) selected up to 10 fascicles (or the number of fascicles in the nerve if the nerve had fewer than 10 fascicles). We modeled 8 fiber diameters (5.7 to 14 µm) at the centroid of each selected fascicle and delivered bipolar stimulation from the cuff. Every modeled fiber consisted of 101 nodes of Ranvier (increased from 53 to reduce the likelihood of end-excitation for optimization and threshold tests), and all modeled fibers had their central node of Ranvier positioned halfway along the modeled nerve.We simulated activation thresholds for S-MF and NEURON models using a bisection search algorithm to determine a threshold within 1% tolerance, where activation was defined as Vm crossing −20 mV with a rising edge at a node of Ranvier 5 nodes from either end of the fiber. We used a time step of 0.005 ms and simulated a total time (tstop) of 5 ms.We evaluated six waveforms: monophasic rectangular, symmetric biphasic rectangular, sawtooth, exponential, sinusoid, and Gaussian (Fig. 6). For each waveform, we tested five pulse widths (\(\delta\)): 0.1, 0.2, 0.5, 0.75, and 1 ms. We compared the activation thresholds predicted by S-MF (TSURR) and by NEURON (TNRN) using the relative percent threshold error (\(100\%\times \frac{{T}_{{{\rm{SURR}}}}-{T}_{{{\rm{NRN}}}}}{{T}_{{{\rm{NRN}}}}}\)). In total, for nerves instrumented with the ImThera cuff, we tested 17,280 thresholds (12 nerves × 6 contact-fascicle pairs × 6 waveforms × 5 pulse widths × 8 fiber diameters), and for nerves instrumented with the helical cuff, we tested 25,200 thresholds (12 nerves × up to 10 fascicles × 6 waveforms × 5 pulse widths × 8 fiber diameters).Fig. 6: Waveform definitions.a Monophasic rectangular. b Symmetric biphasic rectangular. c Sawtooth. d Exponential. e Sinusoid. f Gaussian.We compared the performance of S-MF with the published Peterson threshold estimator14 designed to predict activation thresholds for monophasic rectangular waveforms. Their approach involves calculating a driving force estimator termed MDF2 that uses a weighted sum of second spatial differences over 19 adjacent nodes of Ranvier:$${{{\rm{MDF}}}}_{2}={\sum}_{j}{W}_{\left|n-j\right|}({{\rm{PW}}},d)\cdot {\Delta }^{2}{V}_{{{\rm{e}}}}^{j}$$
(17)
n is the node index for which MDF2 is being calculated. Weights \({W}_{\left|n-j\right|}\) are the ratio of depolarization caused by the current injected at node \(j\) to the depolarization caused by current injected when \(j=n\) for a given pulse width (PW) and fiber diameter \(d\) in a linearized MRG model (generated by replacing all nonlinear ion channels with a fixed specific conductance = 0.007 mS cm−2). Rather than using their published MDF2 values, we reimplemented their procedures to calculate threshold MDF2 for the specific fiber diameters and pulse widths that we tested for our MRG implementation, to avoid any interpolations of fiber diameter or pulse width. We validated our implementation of the Peterson surrogate model by reproducing the published threshold errors across different electrode placements (Supplementary Note 5).Testing—other non-linear responses to stimulationWe evaluated S-MF’s accuracy in reproducing non-linear responses to stimulation other than activation thresholds: unidirectional action potential (AP) propagation, bidirectional direct current (DC) block of AP propagation, AP annihilation by collision, excitation and block responses to kilohertz frequency signals, and state-dependent modulation of activity at conventional stimulation frequencies.For unidirectional propagation and bidirectional DC block, we randomly selected a fascicle and an electrode contact in a pig vagus nerve model instrumented with the ImThera cuff (Supplementary Note 1 P3). We placed a 12 µm diameter fiber in the centroid of the selected fascicle and delivered stimulation using a 0.75 ms cathodic monophasic rectangular pulse from the selected electrode contact. We delivered stimulation at progressively increasingly amplitudes to observe transitions in the neuronal response from activation to unidirectional propagation to block to re-excitation (Fig. 3a) in both the NEURON and S-MF models. We used a timestep (dt) of 0.005 ms and a tstop of 7.5 ms.For kilohertz frequency stimulation, we randomly selected 7 fascicles and 1 electrode contact in a pig vagus nerve model (Supplementary Note 1 P2) instrumented with the ImThera cuff. We placed 5.7, 8.7, and 14 µm diameter fibers at the centroids of the selected fascicles and delivered stimulation from the selected electrode contact. We delivered extracellular sinusoidal stimulation at 1, 2, 5, and 10 kHz starting at t = 0.5 ms across a range of amplitudes (0 to 10 mA for the 5.7 µm fibers, 0 to 5 mA for the 8.7 µm fibers, and 0 to 2 mA for the 14 µm fibers). We initiated action potentials at a rate of 100 Hz at one end of each fiber using a 0.1 ms intracellular anodic current pulse (2 nA) starting at t = 50 ms and recorded the number of action potentials (defined as \({V}_{{{\rm{m}}}}\) crossing −20 mV with a rising edge) arriving at the other end of each fiber after t = 50 ms. In all cases, we simulated a total of 100 ms with a dt of 0.001 ms. We compared the population response (mean number of action potentials and 95% confidence interval across the 7 sampled fiber locations) for each stimulation frequency and amplitude for each fiber diameter between NEURON and S-MF (Fig. 3b).For AP collision, we initiated action potentials at both ends of 5.7, 8.7, 10, and 14 µm fibers using a 2 nA, 0.1 ms intracellular rectangular current pulse. We recorded snapshots of Vm along the length of the fibers at various timepoints to observe how the propagated APs interacted with each other (Fig. 3c).For state-dependent modulation of activity at conventional frequencies, we modeled the effects of extracellular stimulation across a range of suprathreshold and subthreshold amplitudes on fibers firing intrinsically at different mean frequencies due to random intracellular current pulses. We selected a human nerve model (Supplementary Note 1 H2) instrumented with the helical cuff and randomly selected 7 of its fascicles (without replacement). We modeled 5.7, 8.7, and 14 µm diameter fibers (101 nodes each) at the centroid of each fascicle. We first calculated activation thresholds in NEURON for every fiber for a single 0.1 ms symmetric biphasic pulse using a bisection search algorithm, without any intrinsic activity. For both the NEURON and S-MF models, we then stimulated extracellularly with trains of symmetric biphasic pulses at {30, 50, and 100 Hz}, pulse widths of 0.1 ms, and amplitudes from −50% to +50% of the NEURON-calculated mean threshold (across the seven sampled locations for the given fiber diameter). We simultaneously delivered intracellular current pulses to one end of the fiber following five Poisson random generated sequences at each of 3 mean rates µIFR ∈ {10, 50, 100 Hz}, totaling 15 random intrinsic firing patterns for each fiber location and diameter for each extracellular stimulus parameter combination (frequency × amplitude). We simulated for t = 100 * (100/µIFR) ms using a dt of 0.005 ms with and without extracellular stimulation. We recorded \({V}_{{{\rm{m}}}}\) at the fifth node of Ranvier from the end of the fiber distal to the location of intracellular stimulation. We calculated the spike times for APs arriving at the \({V}_{{{\rm{m}}}}\) recording node (the timepoints at which \({V}_{{{\rm{m}}}}\) crossed −20 mV in the positive direction). We compared spike times with and without extracellular stimulation using the SPIKE-synchronization metric64, \({S}_{C}\), which is a mean over aggregated coincidence indicators \(C({t}_{k})\) between two spike trains:$${S}_{C}=\frac{1}{M}{\sum }_{k=1}^{M}C({t}_{k})$$
(18)
where M is the total number of spikes (across both spike trains) and \({t}_{k}\) is the \({k}^{{{\rm{th}}}}\) spike time in the pooled spike train (in the case of an exact match, \(k\) counts both spikes). \(C(t)\) is calculated per spike in the original pair of spike trains and subsequently pooled:$$C\left({{t}_{i}}^{\left(1\right)}\right)=\left\{\begin{array}{c}1\, {{\rm{if}}}\, {\min }_{j}\left(\left|{{t}_{i}}^{\left(1\right)}-{{t}_{j}}^{\left(2\right)}\right|\right)\, < \,{{\tau }_{{ij}}}^{(1,2)}\\ 0\, {{\rm{otherwise}}}\end{array}\right.$$
(19)
$${{\tau }_{{ij}}}^{(1,2)}=\min \left\{{{t}_{i+1}}^{\left(1\right)}-{{t}_{i}}^{\left(1\right)},{{t}_{i}}^{\left(1\right)}-{{t}_{i-1}}^{\left(1\right)},{{t}_{j+1}}^{\left(2\right)}-{{t}_{j}}^{\left(2\right)},{{t}_{j}}^{\left(2\right)}-{{t}_{j-1}}^{\left(2\right)}\right\}/2$$
(20)
where \({{t}_{i}}^{\left(1\right)}\) is the \({i}^{{{\rm{th}}}}\) spike time in train 1 (and likewise for \({{t}_{j}}^{\left(2\right)}\)). \({S}_{C}\) is 0 if and only if the two spike trains do not contain any coincidences, and \({S}_{C}\) is one if and only if each spike in each spike train has exactly one matching spike in the other spike train. We used this metric to measure the extent to which the extracellular signal transformed the intrinsic firing signal (Fig. 3d).Box 1 Algorithm 1 – Training Data GenerationSolve 1 pig and 1 human vagus nerve model instrumented with the ImThera cuff (FEM)ns ← number of fascicles in modeled pig nerveif pig morphology thensample field longitudinally at all fascicle centroids (one fiber location per fascicle)elsesample field longitudinally with uniform density at ns fiber locations within fasciclesdonefor i in 1 to S dofor j in 1 to C do \({{\bf{v}}}\in {{\mathbb{R}}}^{N}\leftarrow\) randomly selected voltage distribution from pig or human nerve \({{{\bf{s}}}}_{{{\rm{t}}}}\in {{\mathbb{R}}}^{T}\leftarrow\) randomly generated monophasic rectangular pulse, A ~ U[−0.2,0.2) mA, ▷ see Fig. 6a δ ~ U[0, 2) ms, ∇ ~ U[0, 2) ms \({{{\bf{v}}}}_{e}^{j}{{\in }}{{\mathbb{R}}}^{T\times N}\leftarrow {{{\bf{s}}}}_{{{\rm{t}}}}\otimes {{\bf{v}}}\)end for\(D \sim U[{\mathrm{5.7,\, 14}}) \,{{\rm{\mu }}}{{\rm{m}}}\)\({{{\bf{b}}}}_{c}{{\in }}{{\mathbb{R}}}^{T\times N}\leftarrow {\sum }_{j=1}^{6}{{{\bf{v}}}}_{e}^{j}\)simulate MRG fiber (N nodes, fiber diameter D) with extracellular potential boundary condition bcsave field-response pair (all N nodes, all T timesteps, all state variables (\({V}_{{{\rm{m}}}},\, m,\, h,\, p,\, s\)))end forBox 2 Algorithm 2 – Disjoint Truncated Backpropagation Through TimeIn the following, we use NumPy array indexing syntax, where a colon ‘:’ indicates all indices along that dimension and ‘−1’ indicates the final indexfor epoch in 1 to 5 dofor i in 1 to \(\left\lceil\frac{{{\rm{training}}}\, {{\rm{set}}}\, {{\rm{size}}}}{B}\right\rceil\) do get \({{{\bf{x}}}}_{i}\), \({{{\bf{y}}}}_{i}\) ▷ sample minibatch i without replacement from training set \(\widehat{{{\bf{yc}}}}\) \(\leftarrow {{\rm{None}}}\) for c in 0 to \(\left\lfloor T/\Delta \right\rfloor\) do ▷ chunk index \({{\bf{xc}}}\leftarrow\) \({{{\bf{x}}}}_{i}[:,:,c\times \Delta :c\times \Delta+\Delta,:]\) ▷ get disjoint chunk of \({{{\bf{x}}}}_{i}\), length ∆ \({{\bf{yc}}}\leftarrow\) \({{{\bf{y}}}}_{i}[:,:,c\times \Delta :c\times \Delta+\Delta,:]\) if \(c\) \(=\) 0 then ▷ get disjoint chunk of \({{{\bf{y}}}}_{i}\), length ∆ \({{\bf{s}}}\leftarrow {{\bf{SSMRG}}}\) ▷ if first chunk, set \({V}_{m}=-80{{\rm{mV}}}\) and gating variables at ∞ else \({{\bf{s}}}\leftarrow \widehat{{{\bf{yc}}}}\left[:,:,-{{\boldsymbol{1}}},:\right]\) done ▷ otherwise, set IC as previous last state predicted by surrogate \(\widehat{{{\bf{yc}}}}\leftarrow {{\rm{S}}}\left({{\bf{xc}}},\, {{\bf{s}}}\right)\) ▷ calculate surrogate model prediction for input xc and IC s \(L\leftarrow {{||}{{\bf{yc}}}-\widehat{{{\bf{yc}}}}{||}}_{{{\rm{MSE}}}}\) ▷ calculate MSE loss between surrogate prediction and NEURON \({{\rm{Adam}}}(\theta,\, {\nabla }_{\theta }L,\, \lambda )\) ▷ perform gradient descent stepend forend forend forOptimization for spatial selectivityWe designed and implemented two optimization methods (gradient-free and gradient-based). We applied both optimization methods to define rectangular pulses and arbitrary waveforms to achieve spatially selective fiber activation in 12 nerve models (Supplementary Note 1 P2-7 and H2-7) instrumented with the six-contact ImThera cuff. We used NEURON to validate all parameters returned by optimizations that used S-MF.The resting inner diameter of the original ImThera cuff model (which was used to generate training data and activation threshold data for surrogate model validation) was up to 1.9× larger than the diameter our human nerve specimens (without shrinkage correction). Therefore, for human nerve optimization tasks, we inflated all nerve morphologies by 25% to correct for 20% shrinkage in the tissue preparation process26 and instrumented the resulting nerves with a shrunk version of the ImThera cuff. We reduced the cuff’s resting inner diameter from 3 mm to 2.5 mm, but retained the cuff length, insulation thickness, number of contacts, contact size, and lengthwise contact spacing. Edge-to-edge contact spacing of the shrunk ImThera cuff was reduced from 1.05 mm to 0.96 mm. Circumferential inter-contact spacing was reduced from 1.34 mm to 1.11 mm; this spacing was preserved as the cuff expanded to fit around nerves with an inner diameter larger than 2.5 mm.Optimization for spatial selectivity: tasksThe overall task in all cases was to maximize activation of target fibers and avoid activation of off-target fibers, with different objective functions implemented based on the optimization method (see “Optimization for spatial selectivity: algorithms”). We defined the functional organization of the pig vagus nerves using immunohistochemistry which showed that afferents and efferents are segregated on separate halves of the nerve41, as confirmed with in vivo imaging65 and physiological recordings35. Target and off-target fiber populations were assigned based on this histological designation of sensory/motor (afferent/efferent) fascicles, with sensory fibers targeted for selective activation (Supplementary Note 1 P2-7). Analogous immunohistochemical data are not yet available for human vagus nerves, but this functional organization reflects the gross anatomical structure of the nodose ganglion—where the nodose contains the cell bodies of visceral vagal afferents, next to which efferent fibers pass—and therefore we assumed a similar vagotopy between the pig and human nerves. In human nerves (Supplementary Note 1 H2-7), we randomly partitioned the nerve approximately into two semicircles and randomly assigned all fascicles in each semicircle as either target or off-target. We placed a single 5.7 µm diameter fiber at the centroid of each fascicle, each with 101 nodes of Ranvier and a node placed at the longitudinal midpoint of the nerve.Stimulation with rectangular pulsesWe optimized the stimulation amplitudes delivered through all six contacts of the ImThera cuff for symmetric biphasic rectangular waveforms (PW = 0.4 ms, delay = 0.4 ms). For each nerve, this yielded an optimization problem with 6 degrees of freedom. We compared optimized performance using the ImThera cuff with the best selectivity predicted for bipolar stimulation with the helical cuff. Specifically, we determined the optimum amplitude for stimulation with the helical cuff (using the same biphasic rectangular waveform) by calculating activation thresholds (in NEURON) for every fiber in the nerve and subsequently determining the threshold amplitude that minimized the weighted binary cross-entropy loss (Eq. (21)).Stimulation with arbitrary waveformsWe optimized the current delivered at each timestep (dt = 0.005 ms) between 0.2 and 1.2 ms (the ‘nonzero period’) from all six contacts of the ImThera cuff. For each nerve, this yielded an optimization problem with 1200 degrees of freedom (200 current amplitudes corresponding to 200 timesteps for each of the six contact). We enforced charge balance for each contact by subtracting the mean of the nonzero period from the nonzero period for each contact before evaluating the neural response.Optimization for spatial selectivity: algorithmsWe developed two optimization methods: gradient-free (using an adapted differential evolution algorithm) and gradient-based (using gradient descent). For the gradient-free method, we compared the performance between the NEURON and S-MF models. Conversely, the gradient-based method leveraged the differentiability of the surrogate; therefore, we only used the surrogate model.Gradient-free optimization (differential evolution)We implemented an adapted differential evolution (DE) algorithm66 to optimize stimulation parameters (Fig. 7). We defined an initial population of candidate solutions (i.e., sets of stimulation parameters) and iteratively updated the candidate solutions based on their performance on the target optimization task (i.e., selective activation of target neurons, as defined in “Optimization for spatial selectivity: tasks”). Thus, the DE approach is stochastic and does not rely on the differentiability of the loss function with respect to the parameters being optimized, i.e., it is gradient-free. We selected DE for its relative ease of implementation, applicability to real-valued loss functions, widespread use, and continued development67.Fig. 7: Algorithmic flowchart for DE/best/1/bin.G = generation number; Gmax = maximum number of generations permitted; Xi,G = population member/candidate solution i in generation G; \({{{\bf{X}}}}_{{{\rm{r}}}0,G}\), \({{{\bf{X}}}}_{{{\rm{r}}}1,G}\) = mutually exclusive random members selected from population in generation G; \({{{\bf{X}}}}_{{{\rm{best}}},G}\) = best population member in generation G; rand(0,1) = random number from uniform distribution in interval [0, 1); j = index for each element in \({{{\bf{X}}}}_{i,G}/{{{\bf{U}}}}_{i,G}/{{{\bf{V}}}}_{i,G}\); CR = crossover rate; F = mutation rate; f = loss criterion; NP = population size. a Mutation operation. b Crossover operation. c Selection operation.We initialized the population of candidate solutions using random Latin hypercube sampling. For all DE optimizations, we used bounds of (−0.3, 0.3) mA. For a population of size \(n\), a random Latin hypercube was constructed by splitting each of the d dimensions to be sampled into n equivalently sized intervals over the respective bounds. Each dimension corresponded to a parameter to be optimized in our model of extracellular peripheral nerve stimulation: for the rectangular pulses, this corresponded to one amplitude for each contact, yielding a 6-dimensional optimization problem, and for arbitrary waveforms, this corresponded to 200 current amplitudes for each contact, yielding a 1200-dimensional optimization problem. We then randomly selected a value from a uniform distribution in each of the n intervals, and we generated the initial population of candidate solutions by randomly shuffling these values to produce n sets of parameters. We used populations of 100 and 300 candidate solutions for the NEURON and S-MF optimizations, respectively.For all optimizations using the DE algorithm, we calculated the loss function as the weighted binary cross-entropy (WBCE):$$f({{\bf{X}}})={WBCE}={\sum }_{n=1}^{N}\frac{{f}_{n}}{F}\left((1-{{A}}_{n})\times {{\mathrm{ln}}}(1-{\hat{{A}}}_{n})+{{A}}_{n}\times {{\mathrm{ln}}}{\hat{{A}}}_{n}\right)$$
(21)
where N is the total number of fascicles in the nerve, fn is the cross-sectional area of fascicle n, F is the summed cross-sectional area of all fascicles in the nerve, \({{A}}_{n}\in \left\{{\mathrm{1,0}}\right\}\) is the target activation (activated or not) of the fiber in fascicle \(n\), and \({\hat{{A}}}_{n}\in \left\{{\mathrm{1,0}}\right\}\) is the predicted activation (using NEURON or S-MF) of the fiber in fascicle \(n\). Fiber activation was determined by the presence of at least one AP, defined as a rising \({V}_{m}\) crossing 0 mV, in either node of Ranvier 5 nodes from each end of the fiber.Each generation, the population of candidate solutions was updated through a process of mutation (Fig. 7a), crossover (Fig. 7b), and selection (Fig. 7c). A “trial vector” \({{{\bf{U}}}}_{i,G}\) was constructed for each \({i}^{{{\rm{th}}}}\) candidate solution \({{{\bf{X}}}}_{i,G}\) that was evaluated in the current generation G using the “best/1/bin”68 mutation rule and crossover operation:$${{{\bf{V}}}}_{i,G}={{{\bf{X}}}}_{{{\rm{best}}},G}+F\cdot ({{{\bf{X}}}}_{{{\rm{r}}}0,G}-{{{\bf{X}}}}_{{{\rm{r}}}1,G})$$
(22)
$${{{\bf{U}}}}\, _{i,G}^{j}=\left\{\begin{array}{c}{{{\bf{V}}}}\, _{i,G}^{j}{\mbox{if}}{{{\rm{rand}}}}_{j}\left(0,\, 1\right)\, < \,{{\rm{CR}}}\\ {{{\bf{X}}}}_{i,G}^{j} \, {\mbox{otherwise}}\end{array}\right.\forall j\in 0,\,1\ldots d$$
(23)
where \({{{\bf{X}}}}_{{{\rm{best}}},G}\) is the best candidate solution in the current generation, F is the mutation rate, \({{{\bf{X}}}}_{{{\rm{r}}}0,G}\) and \({{{\bf{X}}}}_{{{\rm{r}}}1,G}\) are mutually exclusive randomly selected candidate solutions from the current generation, rand(0, 1) is a random sample from a uniform distribution in the interval (0, 1) (resampled for every j), CR is the crossover rate, and d is the dimensionality of the optimization problem. If the resulting trial vector \({{{\bf{U}}}}_{i,G}\) exceeded the problem bounds in any dimension, the mutation and crossover operation was repeated up to 10 times, after which those dimensions which remained outside the bounds were randomly reinitialized to values within the bounds. This is a feasibility-preserving strategy for boundary constraint handling that has shown advantages over Darwinian (i.e., where boundary violations are discouraged by adding a penalty term to the loss function) and repair (i.e., where dimensions in which the boundary has been violated are placed back within the bounds following some rule, such as reflection against the boundary) methods in recent empirical studies using DE69. \({{{\bf{U}}}}_{i,G}\) then replaced \({{{\bf{X}}}}_{i,G}\) in generation G + 1 if \(f({{{\bf{U}}}}_{i,G})\, < \,f({{{\bf{X}}}}_{i,G})\), i.e., if the loss of \({{{\bf{U}}}}_{i,G}\) was less than the loss of \({{{\bf{X}}}}_{i,G}\), or equivalently, \({{{\bf{U}}}}_{i,G}\) had a better performance.We used parameter adaptation to improve convergence70. At each generation, a crossover rate \({{\rm{C}}}{{{\rm{R}}}}_{{{\rm{i}}}}\) was generated independently for each member of the population, sampling from a normal distribution with mean CR and standard deviation of 0.1, and subsequently clipped to [0, 1]. Let \({{{\bf{S}}}}_{{{\rm{CR}}}}\) be the set of all such crossover rates in the current generation for which the resulting trial vector \({{{\bf{U}}}}_{i,G}\) successfully replaced the corresponding target vector \({{{\bf{X}}}}_{i,G}\) (i.e., had a lower loss). In the next generation, \({{\rm{CR}}}\) was updated such that \({{\rm{CR}}}=\left(1-c\right)\times {{\rm{CR}}}+c\times {{\rm{mea}}}{{{\rm{n}}}}_{{{\rm{A}}}}\left({{{\bf{S}}}}_{{{\rm{CR}}}}\right)\), where meanA is the arithmetic mean and c is a positive number between 0 and 1. Likewise, a mutation rate \({F}_{i}\) was generated independently for every member of the population by random sampling from a Cauchy distribution with location parameter F and scale parameter 0.1, truncated to 1 if \({F}_{i}\, > \,1\) and regenerated if \({F}_{i}\le 0\). Using the same terminology as before, in the subsequent generation, F was updated such that \(F=\left(1-c\right)\times F+c\times {{\rm{mea}}}{{{\rm{n}}}}_{{{\rm{L}}}}\left({{{\bf{S}}}}_{F}\right)\), where \({{\rm{mea}}}{{{\rm{n}}}}_{{{\rm{L}}}}\left({{{\bf{S}}}}_{F}\right)\) is the Lehmer mean = \(\frac{{\sum }_{F\in {S_F}}{F}^{2}}{{\sum }_{F\in {S_F}}F}\). We initialized F to 0.8, CR to 1.0, and we used a value of c = 0.1 as in the original paper.We terminated the DE operation after 200 or 500 generations (Gmax) for the NEURON and S-MF optimizations, respectively, or if a solution with a theoretic minimum WBCE loss (0) was found.Gradient-based optimization (gradient descent)Gradient descent operates on scalar-valued differentiable loss functions. Surrogate models implemented using AxonML are differentiable and implemented in PyTorch; therefore, neuromodulation parameters can be designed using gradient descent for certain optimization tasks.We optimized stimulation parameters using the procedure outlined in Algorithm 3 (Box 3). Let \({{\rm{S}}}:{{\bf{x}}}\in {{\mathbb{R}}}^{B\times N\times T\times 1},{{\bf{s}}}\in {{\mathbb{R}}}^{B\times N\times 5}\to {\hat{{{\bf{y}}}}{{\boldsymbol{\in }}}{\mathbb{R}}}^{B\times N\times T\times 5}\) represent the mapping performed by S-MF S from input x (field), to predicted output \(\hat{{{\bf{y}}}}\) (the fiber response, i.e., membrane potential and gating variables for all nodes at all simulated timepoints) with initial condition s, where B is the number of fibers in the nerve model, N is the number of modeled nodes of Ranvier per fiber, and T is the number of simulated timepoints. Let SSMRG (steady-state MRG) represent the state of MRG nodal compartments at rest (Vm = −80 mV and all four gating variables \(x\) at \({x}_{\infty }(-80\,{{\rm{mV}}})\)).Let \({{\rm{MF}}}:{{{\bf{p}}}\in{\mathbb{R}}}^{d}\to {{\bf{x}}}\in {{\mathbb{R}}}^{B\times N\times T\times 1}\) represent the mapping from the set of optimizable stimulus parameters \({{\bf{p}}}\) (where d is the degrees of freedom; \(d\) = 6 in the waveform constrained case and \(d\) = 1200 in the arbitrary waveform case) to the S-MF input \({{\bf{x}}}\) (tensor representing the extracellular field at every node of ever fiber at every timestep). Let \({{\rm{Ranger}}}({{\bf{p}}},{\nabla }_{{{\bf{p}}}}L,\lambda )\) represent a single gradient descent update to optimizable stimulus parameters \({{\bf{p}}}\) using the Ranger update rule71 (Rectified Adam72 with LookAhead73 and Gradient Centralization74) and learning rate \(\lambda\), where gradients with respect to the stimulus parameters \({{\bf{p}}}\) of a scalar valued loss \(L\) (weighted quotient (WQ) loss, Eq. (26)) are represented by \({\nabla }_{{{\bf{p}}}}L\). Let \({{\rm{WBCE}}}:{{{\bf{A}}}}\in {\left\{{0,\,1}\right\}}^{N},\hat{{{\bf {A}}}}\in {\left\{{0,\,1}\right\}}^{N}{\mathbb{\to }}{\mathbb{R}}\) represent the weighted binary cross-entropy (Eq. (21)) between target activations \({{\bf{{A}}}}\) and predicted activations \(\,\!\hat{{{\bf{{A}}}}}\) (where activation, in contrast to \(\hat{{{\bf{y}}}}\), is a derived Boolean quantity indicating whether the fiber propagated an action potential, defined, as before, by \({V}_{m}\) crossing −20 mV in the positive direction at 5 nodes from either end of the fiber). Finally, let \({{\rm{WQ}}}:{\hat{{{\bf{y}}}}\in{\mathbb{R}}}^{B\times N\times T\times 5}{\mathbb{\to }}{\mathbb{R}}\) be the weighted quotient loss on S-MF prediction \(\hat{{{\bf{y}}}}\):$${m}_{{{\rm{on}}}}(\hat{{{\bf{y}}}})={\sum}_{n\in {{{\bf{n}}}}_{{{\rm{on}}}}}{f}_{n}\times {m}_{n}(\hat{{{\bf{y}}}})$$
(24)
$${m}_{{{\rm{off}}}}(\hat{{{\bf{y}}}})={\sum}_{n\in {{{\bf{n}}}}_{{{\rm{off}}}}}{f}_{n}\times {m}_{n}(\hat{{{\bf{y}}}})$$
(25)
$${{\rm{WQ}}}(\hat{{{\bf{y}}}})=\sqrt{\frac{d}{{n}_{{{\rm{c}}}}}}\frac{{m}_{{{\rm{off}}}}(\hat{{{\bf{y}}}})}{{m}_{{{\rm{on}}}}(\hat{{{\bf{y}}}})}$$
(26)
where \({{{\bf{n}}}}_{{{\rm{on}}}}\) are the indices of all target fascicles, \({{{\bf{n}}}}_{{{\rm{off}}}}\) are the indices of all off-target fascicles, \({f}_{n}\) is the cross-sectional area of fascicle \(n\), \({m}_{n}(\hat{{{\bf{y}}}})\) is the sum over the \(m\) gating parameter in \(\hat{{{\bf{y}}}}\) corresponding to the 20 end nodes of Ranvier (10 at each end, including the terminal nodes) of the fiber in fascicle \(n\), and \({n}_{{{\rm{c}}}}\) (n = 6) is the number of electrode contacts in the model for which there were optimizable parameters.We leveraged the concurrency of our S-MF GPU implementation to evaluate the neural response and calculate gradients across all 12 nerve models in parallel. This allowed us to solve all 12 optimization problems, with rectangular pulses or arbitrary waveforms, simultaneously with modest overhead.Box 3 Algorithm 3 – Gradient Descent for Stimulus Parameter OptimizationWe used wd = 0.01*N (where N is the number of fascicles in the nerve), γmax = 200/N, η = 0.6, and ns = 200Requires wd, η, ns, \({\gamma }_{\max }\), \({{{\bf{ A}}}}\) ▷ weight decay, learning rate decay, # gradient descent steps, maximum gradient norm, target fiber activations\({{\bf{p}}}\, {\leftarrow }\,\{0\}^{d}\) ▷ initialize stimulus parameters to 0\(\lambda \leftarrow 2\) ▷ initialize learning rate to 2\({{{\bf{p}}}}_{{{\rm{best}}}},\, {\omega }_{{{\rm{best}}}},\, {{{\mathscr{L}}}}_{{{\rm{best}}}}\leftarrow\) Nonefor i in 1 to ns do ▷ optimization loop \({{\bf{x}}}{{\boldsymbol{\leftarrow }}}\, {{\rm{MF}}}({{\boldsymbol{p}}})\) ▷ construct field from stimulus parameters \(\hat{{{\bf{y}}}}{{\boldsymbol{\leftarrow }}}\, {{\rm{S}}}({{\bf{x}}}{{\boldsymbol{,}}}\, {{\bf{SSMRG}}})\) ▷ use surrogate to predict neural response \(\omega \leftarrow {{\rm{WBCE}}}({{\bf{{A}}}}{{,}}\, \hat{{{\bf{ A}}}})\) ▷ calculate weighted binary cross entropy loss \({{\mathscr{L}}}{{\mathscr{\leftarrow }}}{{\rm{WQ}}}(\hat{{{\boldsymbol{y}}}})\) ▷ calculate weighted quotient if i = 1 then \({{{\bf{p}}}}_{{{\rm{best}}}},\, {\omega }_{{{\rm{best}}}},\, {{{\mathscr{L}}}}_{{{\rm{best}}}}\leftarrow {{\bf{p}}},\,\omega {{\mathscr{,}}}\,{{\mathscr{L}}}\) ▷ initialize best params and losses to initial values done if \(\omega \, < \,{\omega }_{{{\rm{best}}}}\) then ▷ if improvement in WBCE \({{{\bf{p}}}}_{{{\rm{best}}}},\, {\omega }_{{{\rm{best}}}}\leftarrow {{\bf{p}}},\, \omega\) ▷ update best params and WBCE done if \(\omega={\omega }_{{{\rm{best}}}}\) and \({{\mathscr{L}}}\, {{ < }}\, {{{\mathscr{L}}}}_{{{\rm{best}}}}\) then ▷ if improvement in WQ and WBCE still best \({{{\bf{p}}}}_{{{\rm{best}}}},\, {\omega }_{{{\rm{best}}}},\, {{{\mathscr{L}}}}_{{{\rm{best}}}}\leftarrow {{\bf{p}}},\, \omega {{\mathscr{,}}}\, {{\mathscr{L}}}\) ▷ update best params, WBCE, and WQ if \({\omega }_{{{\rm{best}}}} \, < \, 1\) then \(\lambda \, {\leftarrow }\, \lambda*\eta\) ▷ reduce learning rate if near minimum done done \({{\boldsymbol{\gamma }}}_{{\bf{p}}} \leftarrow \nabla_{{\boldsymbol{\!p}}}{{\mathscr{L}}}\) ▷ calculate gradients \({{{\boldsymbol{\gamma }}}}_{{{\bf{p}}}}\leftarrow \min \left(1,\frac{{\gamma }_{\max }}{{||}{{{\boldsymbol{\gamma }}}}_{{{\bf{p}}}}{||}}\right){{{\boldsymbol{\gamma }}}}_{{{\bf{p}}}}\) ▷ clip gradient norm \({{\bf{p}}}\leftarrow {{\bf{p}}}\odot(1\, {{\boldsymbol{-}}}\, {{\rm{wd}}})\) ▷ apply weight decay on stimulus parameters \({{\rm{Ranger}}}({{\bf{p}}},\, {{{\boldsymbol{\gamma }}}}_{{{\bf{p}}}},\, \lambda )\) ▷ perform gradient descent step on weighted quotient lossend forreturn \({{{\boldsymbol{p}}}}_{{{\rm{best}}}}\) ▷ return best params (corresponds to min. WQ for min. WBCE)Statistical analysisWe compared performances of the various optimization methods using a non-parametric paired test (two-sided Wilcoxon signed-rank test) between the WBCE values of generated solutions. We used a non-parametric test because most distributions of WBCE values were non-Gaussian (as verified by the Shapiro–Wilk test). We adjusted significance values (with an uncorrected critical \(p\) value = 0.05) using the Holm-Bonferroni method75 to account for multiple comparisons; we reported both \(p\) and corrected significance values. All Boxplots are displayed using the Tukey method (center line, median; box limits, upper and lower quartiles; whiskers, last point within 1.5× interquartile range of nearest hinge).Reporting summaryFurther information on research design is available in the Nature Portfolio Reporting Summary linked to this article.
Highly efficient modeling and optimization of neural fiber responses to electrical stimulation
