A Huygens’ surface approach to rapid characterization of peripheral nerve stimulation

Purpose: Peripheral nerve stimulation (PNS) modeling has a potential role in designing and operating MRI gradient coils but requires computationally demanding simulations of electromagnetic fields and neural responses. We demonstrate compression of an electromagnetic and neurodynamic model into a single versatile PNS matrix (P-matrix) defined on an intermediary Huygens’ surface to allow fast PNS characterization of arbitrary coil geometries and body positions. Methods: The Huygens’ surface approach divides PNS prediction into an extensive precomputation phase of the electromagnetic and neurodynamic responses, which is independent of coil geometry and patient position, and a fast coil-specific linear projection step connecting this information to a specific coil geometry. We validate the Huygens’ approach by performing PNS characterizations for 21 body and head gradients and comparing them with full electromagnetic-neurodynamic modeling. We demonstrate the value of Huygens’ surface-based PNS modeling by characterizing PNS-optimized coil windings for a wide range of patient positions and poses in two body models. Results: The PNS prediction using the Huygens’ P-matrix takes less than a minute (instead of hours to days) without compromising numerical accuracy (error ≤ 0.1%) compared to the full simulation. Using this tool, we demonstrate that coils optimized for PNS at the brain landmark using a male model can also improve PNS for other imaging applications (cardiac, abdominal, pelvic, and knee imaging) in both male and female models. Conclusion: Representing PNS information on a Huygens’ surface extended the approach’s ability to assess PNS across body positions and models and test the robustness of PNS optimization in gradient design.

Eq. 1 The magnetic vector potential is computed from a given current source using the Biot-Savart formulation in integral form: Eq. 2 We then obtain the electric scalar potential by solving the differential form where is the spatially varying tissue conductivity. The electric scalar potential is modeled by a hexahedral finite element mesh (1 mm 3 cell size) yielding 76M mesh elements in our male model and 48.8M in the female model. Each cell is labeled by the tissue-specific conductivity obtained from the IT'IS LF database. We dramatically reduced computation time by taking advantage of the fact that the body model mesh is fixed for all Huygens' basis functions. This allowed us to reuse the FE system matrix on the left-hand side of Eq. 6 after the initial assembly (the most computationally intensive step). In our implementation, the E-field for each Huygens' basis can be computed in 2-3 minutes. The fast solver was validated by comparing the fields to the low-frequency solver in a commercial package (Sim4Life, Zurich MedTech).

PNS oracle calculation
The E-field of each basis function was projected onto the fibers of the body's nerve atlas. Integration along the nerve generated the electric potentials ( ) along each of the ∼1900 nerves in the body. Using the pre-assigned axon diameter for each nerve, we computed the PNS oracle metric. The PNS oracle corresponds to the reciprocal of the PNS threshold coil current in units of Ampere, i.e., the PNS oracle has units of 1/A but can be easily related to the expected gradient strength at that current. The oracle estimates a nerve's PNS threshold based on the spatial characteristics of the electric potential ( ) induced along the nerve path and the axon diameter of the nerve. The PNS oracle is defined by where ( ) is the node of Ranvier distance (which is a function of the axon diameter ), ( ) is a spatial kernel, and ( ) is the myelination calibration factor. The PNS oracle metric is based on the second spatial difference of the electric potential ( ) across consecutive nodes of Ranvier (similar to previous definitions of the neural activating function and modified driving function). The result is convolved with a kernel ( ) to describe current redistribution across nodes of Ranvier. The factor ( ) accounts for the effect of different levels of myelination on the nerve's excitability. Both the spatial kernel ( ) and the myelination factor ( ) are calibrated for a specific coil current waveform: in this work we use a trapezoidal waveforms of varying rise times (500 μs flattop duration, 32 bipolar pulses). This calibration was done by comparing to the numerical titration of the double-cable nerve model in 64,000 nerve segments.
The oracle computation only uses operations linear in the applied current. This ensures that the PNS oracle for a coil configuration can be represented as a linear combination of the oracles of the individual basis functions using the same weights that would be used to formulate the coil configuration's B-field from its basis function B-fields. The linearity also allows the basis function PNS responses to be assembled into a matrix form (the P-matrix). Given a vector listing the weights of the basis functions for a specific winding pattern, the total PNS response and resulting PNS threshold can be quickly computed using total = and PNS threshold = 1 max{| total |} Eq. 5 This P-matrix product can be incorporated as a constraint in the gradient winding optimization or potentially used as a real-time safety watchdog during scanner operation.

Computational runtime
The Huygens' P-matrix computations fall into two categories, computationally intensive steps intended to be performed once and computationally efficient steps intended to be performed for each coil geometry and patient position. The computationally intensive steps included the pre-calculation step of the Huygens' P-matrix which consists of EM calculations and compilation of the PNS oracles into the Pmatrix. These steps were performed on a high-end Linux server (using 300 GB of memory and 20 cores), taking 3-5 days per body model. The simpler calculations are expected to be done multiple times (e.g. for every body position) and were done a high-end laptop (using 8 GB of memory and 6 cores). These simpler steps included computing the mapping matrix M for a given body position, application of this matrix to the Huygens' P-matrix and evaluation of PNS thresholds by multiplying the resulting P-matrix by the current weight vector describing the coil configuration. Evaluating PNS for a given coil winding pattern is fast (∼5 seconds) since in this case the Huygens' basis set is mapped to a single "excitation source", i.e., a single gradient coil. The mapping matrix calculation required for coil winding optimization is more expensive, as it requires a mapping of the Huygens' basis set to the thousands of "excitation sources", i.e., the current density elements on the coil former (7595 in our example), requiring ∼1 minute.
Evaluation of PNS thresholds for a given trial current basis set is a very fast matrix-multiply of the coil geometry specific P-matrix PC with the candidate basis vector (requiring milliseconds). The FEM electromagnetic simulations constitute by far the most expensive step of the precomputation phase. As previously described, these field simulations require relatively high spatial resolution (∼1 mm3) to correctly capture local E-field hotspots frequently occurring, e.g., in thin layers of fatty tissue through which most nerves run in the body. This leads to large FEM simulation domains with up to 76M mesh cells. Additionally, the body model simulation domain cannot be truncated to reduce computation time. We found, that truncation of the simulation domain (e.g., excluding portions of the model with negligible B-field amplitudes) led to unpredictable constructive or destructive error propagation and incorrect PNS thresholds. We dealt with the additional computational demand incurred in the whole-body simulations by developing our own in-house EM field solver based on the lightweight, heavily parallelized MFEM C++ library. This allowed us to reduce computation time from ∼45 minutes to 2-3 minutes per Huygens' basis. Supplementary Table S1: Summary of PNS threshold curve parameters ΔGmin(0-p) (minimum stimulating waveform zero-to-peak gradient amplitude), tchron (chronaxie time), as well as ΔGmin(p-p) (minimum stimulating waveform peak-to-peak gradient excursion) and SRmin (minimum slew-rate needed to achieve stimulation). The PNS threshold curves were obtained by simulation of trapezoidal waveforms with 500 μs plateau time and 16 bipolar cycles.