Overcoming temporal dispersion for measurement of activity-related impedance changes in unmyelinated nerves

Objective. Fast neural electrical impedance tomography is an imaging technique that has been successful in visualising electrically evoked activity of myelinated fibres in peripheral nerves by measurement of the impedance changes (dZ) accompanying excitation. However, imaging of unmyelinated fibres is challenging due to temporal dispersion (TP) which occurs due to variability in conduction velocities of the fibres and leads to a decrease of the signal below the noise with distance from the stimulus. To overcome TP and allow electrical impedance tomography imaging in unmyelinated nerves, a new experimental and signal processing paradigm is required allowing dZ measurement further from the site of stimulation than compound neural activity is visible. The development of such a paradigm was the main objective of this study. Approach. A finite element-based statistical model of TP in porcine subdiaphragmatic nerve was developed and experimentally validated ex-vivo. Two paradigms for nerve stimulation and processing of the resulting data—continuous stimulation and trains of stimuli, were implemented; the optimal paradigm for recording dispersed dZ in unmyelinated nerves was determined. Main results. While continuous stimulation and coherent spikes averaging led to higher signal-to-noise ratios (SNRs) at close distances from the stimulus, stimulation by trains was more consistent across distances and allowed dZ measurement at up to 15 cm from the stimulus (SNR = 1.8 ± 0.8) if averaged for 30 min. Significance. The study develops a method that for the first time allows measurement of dZ in unmyelinated nerves in simulation and experiment, at the distances where compound action potentials are fully dispersed.


Introduction
Electroceuticals or bioelectronic medicines [1] are a novel emerging set of techniques aimed at treating diseases by selective stimulation of nervous tissue and neuromodulation of the internal organs innervated by it. The main nerve of interest is the vagus nerve (VNS) which is the longest autonomic nerve in the body serving as an interface between the central nervous system and major internal organs. Electrical stimulation of the VNS and neuromodulation of the internal organs or areas of the brain supplied by it is a clinically approved technique for treatment of drug-resistant epilepsy [2,3] and depression [4], and has a great potential for treatment of heart failure [5], rheumatoid arthritis [6] and a variety of inflammatory diseases via modulation of the cholinergic antiinflammatory pathway [7][8][9]. However, despite the good clinical efficacy and great potential of VNS, it is prone to adverse effects with an incidence rate of up to 50% [10,11] originating from non-targeted electrical stimulation of the VNS that induces modulation of all organs supplied by it leading to undesired physiological effects. To reduce side effects, selective stimulation of the specific fascicle within the nerve leading to a particular organ in the single direction of this organ is required.
To allow selective stimulation of fascicles in the nerve, it is essential to know their precise location within this nerve. Moreover, organisation of a closed feedback loop neuromodulation of the desired organ relies on the measurement of functional activity of the fascicle supplying this organ. Both capabilities can be achieved with fast neural electrical impedance tomography (EIT), a novel method capable of imaging electrical activity in nerves in their cross-section.

Temporal dispersion
Although fast neural EIT is reliable for imaging neural activity in the brain and peripheral nerves, imaging in autonomic nerves is more challenging for two reasons. First, autonomic nerves mainly consist of small unmyelinated C fibres [23,24] producing lower impedance changes; second, conduction velocities (CVs) of C fibres are significantly slower and more variable [25,26]. As a result, the amplitude of the compound action potentials (CAPs) being an aggregate sum of action potentials (APs) of individual fibres rapidly decreases along the nerve from the site of its activation. This effect is referred to as temporal dispersion (TD) [26][27][28][29] which leads to a fall of the CAP below the noise threshold beyond a few centimetres or in some cases even a few millimetres from the stimulus [22,26,30].
The last point is especially important for the achievement of the goals stated for bioelectronic medicines. First, mammalian, including human, VNSs are largely unmyelinated [23,31]. Second, for selective stimulation of the specific fascicles and closed-loop neuromodulation of the internal organ which they supply, the neural activity propagating from the organ must be recorded and imaged at the cervical level. The length of the nerve from the internal organs to the neck in humans is around half a metre [32] that is significantly larger than the theoretical limit of approximately 4 cm allowed by TD in unmyelinated fibres [22]. To overcome this limitation, a method capable of measuring dZ further from the onset than allowed by TD is required.
The feasibility of this objective can partly be justified by the fact that CAPs in nerves have more pronounced phasic nature than dZs which are mainly monophasic, as was measured experimentally in crabs [33], rats [12,19] and large animals [34] as well as confirmed in the current study (figure 6 in section 3). Therefore, the expectation is that CAPs will decrease in amplitude much faster than dZ so that dZ will be visible further from the stimulation point than CAPs are.
In addition, study [35] showed that when the stimulation paradigm had been changed from a traditional continuous nerve stimulation to application of a high-frequency series of stimuli separated by resting intervals, it was theoretically possible to record dZ at 20 cm from the stimulus with SNR of 4 if averaging for 30 min, while the traditional approach cannot be used at >5 cm [22]. However, this approach was based on a simple statistical model of an arbitrary C fibre nerve which significantly differs from real mammalian autonomic nerves. Also, the experimental parameters such as the noise, geometry of the electrodes and electrical parameters were arbitrarily chosen based on the previously performed experiments in peripheral nerves [12].
Therefore, in order to provide accurate predictions and determine whether it was realistic to measure dZ in autonomic nerves at distances from the stimulation where CAPs were cancelled out, this approach must have been optimised. It should have accounted for realistic nerve histology and fibre composition as well as included a variety of experimental parameters utilised in impedance measurement experiments performed with mammalian nerves. Then, following the development and optimisation of the method, the predictions must have been verified experimentally.

Purpose
The main purpose of this study was to develop and optimise a method to overcome TD and allow EIT recordings of phasic activity at far distances from the onset along the nerve where compound activity is dispersed. This brings up the opportunity to image fascicles of the autonomic nerves with fast neural EIT at the distances where the CAPs fall below the noise and facilitate the development of bioelectronic medicines for selective stimulation of the VNS, neuromodulation of internal organs and treatment of associated drug-resistant disorders. Specific questions addressed in this study were as follows: (a) What are the optimal stimulation and signal processing strategies and parameters producing the largest impedance changes at different distances (15, 20 and 50 cm) from the onset? (b) How much averaging is required (1) to obtain a measurable signal (SNR > 1) and (2) to image neural activity with EIT (which requires SNR > 4 [18]) at 15, 20 and 50 cm from the site of stimulation? (c) Are simulated results confirmed with the experimental data?
The first part of the study included the development of the accurate model of dispersion in the porcine subdiaphragmatic nerve (SN) which was followed by the development and optimisation of the method for overcoming dispersion and recording dZ at far distances from the stimulation site.

Experimental design
The study was divided into the following steps: (a) Development of an experimentally driven statistical model of the SN of the pig.
The model combined (a) previously developed accurate finite element (FEM) model of a C fibre [36] and (b) statistical model for simulation of TD of CAPs and dZs [22] in a complex nerve, consisting of a composition of unmyelinated and myelinated fibres [37]. The parameters of the model were chosen on the basis of ex-vivo experimental recordings (CAP and dZ) obtained using a SN of the pig subjected to repetitive continuous stimulation. As a result, the model could accurately simulate CAPs and impedance changes at ∼3 cm from the stimulus and could be therefore utilized for the development and optimisation of a new method for overcoming TD and measurement dZ at further distances, where CAPs were cancelled out due to dispersion. (b) Development of a method for overcoming dispersion. Using the developed model, the optimal stimulation and signal processing paradigms to record impedance changes at the distances where compound APs were dispersed were obtained. For this, the modelled nerve fibre was subjected to repetitive stimuli at various frequencies. For dZ extraction, the resultant signals were processed in two ways-(1) averaging and band-pass filtering as single spikes and (2) averaging as trains of spikes allowing band-pass filtering around the whole train thus significantly reducing the bandwidth of the filter. The optimal parameters for dZ measurement at various distances from stimulation were determined, the SNR at these distances was obtained. The final predictions of the model have then been subjected to experimental verification using a preparation of the porcine SN ex vivo (N = 18).

Versatile statistical model of dispersion in nerve 2.2.1. FEM model of a single fibre
Impedance changes accompanying neural activity were obtained using the FEM models of mammalian C fibre (d = 1 µm) bi-directionally coupled with the extracellular space (figure 1) [36]. Two variants of the FEM model were designed in the study. In the first variant ( figure 1(a)), the C fibre (d = 1 µm) was surrounded by a cylinder of extracellular space with the diameter equal to the one of the ring electrodes (D = 0.01 mm). The second variant ( figure 1(b)) was designed to represent the ex-vivo experiments performed using multielectrode cuffs. When the silicone cuff is placed tightly around the nerve, as in the ex-vivo experiment performed in the current study (figure 2), there is a significantly smaller amount of saline solution inside the cuff (in the electrodenerve interface) than along other parts of the nerve. The shape of the recorded CAP, in this case, differs from the one recorded using hook electrodes [33], where the saline solution occupies uniform volume along the nerve (as in figure 1(a)). Therefore, to simulate the activity of the C fibre in the conditions Figure 1. Axisymmetric representation of the model of a mammalian C fibre, based on [36,38]. The axon is depicted by a blue line, the axis of symmetry is shown by the red dash-dotted line. The AP was induced from the end of the fibre; DC or AC were applied through two external electrodes (blue) placed 7.8 and 7.82 mm from the axon's end (distance between injecting electrodes ∆xI = 0.02 mm); the electric field was recorded by an external electrode (green) placed before the injecting ones, 7.7 mm from the proximal end of the fibre (distance between recording and injecting electrodes ∆xR = 0.1 mm). (a) Initial model [36]; (b) the same model with the change in the structure of extracellular space to simulate the condition of ex-vivo experiments with the nerve cuff. For this, the width of the extracellular space was increased 50-fold outside the region along the nerve where electrodes were located (L cuff = 0.3 mm). Shapes of extracellular action potential (EAP) and impedance change (dZ) produced with each model as well as their triangular FEM meshes are depicted on the right. For dZ measurement, constant sinusoidal current (1-6 kHz, 200-300 µA) was injected through two last electrodes on each cuff, and voltage was measured on the remaining electrodes in respect to the last electrode on the last cuff. The 2nd cuff which was closest to the stimulation site was used for control and model development as C fibre CAPs were clearly measurable there, cuffs 3 and 4 were far enough so that the APs corresponding to C fibres were dispersed, only myelinated fibres' CAPs were visible. (b) Multi-electrode cuff design. Six identical electrodes per cuff with a surface area of 0.46 mm 2 each were used. similar to those when the cuff is present, the width of the extracellular space cylinder was increased 50fold outside the 0.3 mm region along the nerve where electrodes were located ( figure 1(b)). The electric field in the extracellular space was simulated using volume conduction Poisson's equation, and the fibre was modelled using the conductance-based Hodgkin-Huxley type Tigerholm model of a mammalian C nociceptor [38]. The Tigerholm model was chosen as it could accurately represent the active behaviour of mammalian C fibre (porcine SNs mainly consist of C fibres [31]) and was validated experimentally [38]. In addition, the model has previously shown itself capable of optimising experimental parameters for obtaining the largest impedance changes in unmyelinated nerves as well as explaining the biophysical origin of experimental recordings [36].
The Tigerholm model contained ten ion channels and variable concentrations inside and outside of the fibre; these concentrations, together with other parameters, affected the membrane conductance which, in turn, led to changes in electric potentials inside and outside the membrane (equations (27) and (28) in the supplementary material of [36]). The impedance change (dZ) signal measured in the model is related to a change in the membrane conductance which is, among other parameters, influenced by variable extra-and intracellular ionic concentrations.
For measurement of impedance change, constant direct current of small amplitude (1.25 µA, or 4 mA cm −2 ) was continuously applied to the fibre during stimulation through two external ring electrodes. The reason was that using DC was less computationally intensive than AC, especially at high frequencies, and it was previously shown that dZ obtained with small amplitude DC and AC had the same biophysical origin [36].
Voltage was measured simultaneously using the recording electrode placed before the injecting ones with respect to ground (figure 1), the same as in [36]. Since the constant current was applied, and because the phase shift ∆φ between the injected current and measured voltage is close to zero [34,39], the measured impedance change dZ may be expressed as follows: In the equation, impedance change dZ is equal to the relative change of the impedance Z(t AP ) when AP passes under the electrodes with respect to the baseline impedance of the system Z = Z(t). Using (1), the complex dZ and absolute |dZ| can be expressed in terms of the measured voltages V = V(t) and V(t AP ). To extract dZ, DC was applied twice in different polarities with positive and negative electrodes switched [36]. These voltage signals were subtracted from each other so that the identical APs were cancelled out, while the dZ which modulates the recorded voltage, doubled, and could be easily extracted.
In the experiments performed in the study, this subtraction step was not required as high AC frequencies (⩾2 kHz) were used and CAPs were removed by bandpass filtering around the carrier frequency (as the characteristic frequency of the CAP is <1 kHz [40]) ( figure 4).
The parameters and geometrical design of the model as well as all the equations describing it can be found in [36].

Statistical model of a complex nerve
The APs and impedance changes (dZs) obtained with the FEM model were incorporated into the statistical model for simulation of the dispersed compound signals of a multi-fibre porcine SN. The modelled nerve included 40 000 unmyelinated C fibres and 4000 fast myelinated fibres [37] fibres uniformly distributed in the cross-section of the nerve, with normally distributed CVs (table 2). Normal distribution of the CV was considered based on the histological studies in various nerves [23,41,42] where the distributions of fibres' diameters (which are directly proportional to CVs [22]) close to normal were obtained; the same assumption was also made in the recent modelling study [22].
The compound dZ was obtained as the sum of single impedance changes of each fibre at the required distance from the site of stimulation: where CAP and dZ are APs and impedance changes of the compound nerve, EAP i C , EAP i M , dZ i C and dZ i M represent single extracellular AP and dZ corresponding to C fibres and myelinated fibres respectively (figure 1).
Due to the difference in geometric and electrical parameters utilised in the single-fibre FEM and the experimentally-based statistical multiplefibre model, the single APs and dZs simulated with the FEM model (subsection A) were scaled in amplitude before summation, as in [22]. However, compared to [22] where electrode diameter and connective tissue resistivity were only considered, multiple parameters differing between the FEM and statistical models were introduced, so that the predictions of the final model closely agree with experimental data.
For scaling, the obtained extracellular AP and dZ were multiplied by the coefficients k AP = ∏ i k APi and where k APi and k dZi accounted for various parameters that differed in the FEM and statistical models and affected AP and dZ respectively: The considered parameters affecting the dZ were [36]: the electrodes' diameter (d), extracellular resistivity (ρ), the current density of the applied dZ measuring current (J), the distance between recording and injecting electrodes (∆x R , figure 1), the distance between injecting electrodes (∆x I , figure 1). The latter three parameters were found to strongly influence the dZ measurements in the previously performed study [36]. For scaling APs, the electrodes' diameter (distance between the electrode and the fibre) and conductivity of the extracellular space were considered, as other parameters did not have any effect on its value. The final equations for scaling APs and dZs are presented below: In the equations, index stat corresponds to the value used in the experiments and in the statistical model, and FEM-in the FEM model. The detailed strategy for evaluation of the scaling coefficients as well as their values are presented in supplementary material (available online at stacks.iop.org/JNE/19/026054/mmedia) and in table 1(A). In addition to the introduced scaling coefficients, CVs of the fibres in the model were chosen so that simulated CAP and dZ closely match the experimental ones recorded using the ex vivo preparation of the porcine SN. The porcine SN mainly consists of unmyelinated C fibres and is a good representation of the subdiaphragmatic branches of human VNSs [31].
For the preparation, nerves of 20-25 cm length were sourced from the terminally anaesthetized pigs used in other experimental studies. The nerves were held in an organ bath perfusion chamber filled with oxygenated Krebs-Ringer solution kept at ∼30 • C. Three silicone rubber cuffs with six radially arranged electrodes each were placed around the nerve 3, 15 & 20 cm from a cuff for electrical stimulation so that the TD of CAPs can be observed (figures 2(a) and (b)). The stainless-steel electrodes (0.2 × 2.3 mm 2 , figure 2(b)) embedded into a medical-grade silicone rubber base were fabricated using a laser cutter and coated with PEDOT:pTS providing the lowest contact impedance and phase shift (∼300 Ω and 1.5 • at 1 kHz) among the popular coating electrode materials [34].
The choice of six-electrode design was done (1) for verification purposes-dZ must be equal on electrodes 1-4 if the measuring current is applied through electrodes 5 and 6, and (2) to account for possible failures of one or more electrodes on the cuff. These may have happened due to multiple reasons including un-soldered connection, broken wire, increased contact impedance due to initially low-quality or detached PEDOT:pTS coating as well as bad contact of the particular electrode and the nerve.
With the designed setup, dZ and CAPs were recorded using the 2nd cuff placed at approximately 3 cm from the site of the onset with respect to the electrode on the last cuff using continuous stimulation with frequency f stim = 2 Hz, current I stim = 20-40 mA depending on the thickness of the nerve, pulse width PW = 50 µs, frequency and amplitude of the applied impedance measuring current f AC = 1-6 kHz, I AC = 200-300 µA (28 nerves in total).
To obtain satisfactory SNR, averaging was required; for this purpose, CAPs were recorded for 20 s, and the dZ-for 10 min. Based on the recorded CAP and dZ, the CV (mean and S.D.) of the fibres constituting the statistical model were determined. This was done using times of the negative peaks of the CAPs related to fast (Ab) and slow (C) fibres as well as their widths. Knowing the distances of the recording cuff from the stimulation, the CV could be calculated.
The average level of the Gaussian noise present in the ex-vivo experiments (3.5 µV RMS before averaging) was added to the resultant modelled signals as the last step. Addition of noise allowed to determine the optimal parameters for maximisation of the SNR, defined as the amplitude of the signal divided by the standard deviation of the noise. The optimal parameters to be determined include, among others, the bandwidth of filtering which will be influenced by type and levels of noise present in the recordings.
The resulting finalised model could utilise EAP and dZ produced with the FEM model ( figure 1(b)) subjected to an arbitrary stimulation paradigm and, on this basis, it could accurately predict values of CAP and dZ produced by the porcine SN at any distance from the onset of the stimulus. Therefore, the model is versatile and could be used to determine the optimal stimulation paradigm for overcoming dispersion and recording dZ further than CAPs are visible.
The COMSOL and MATLAB model files used for FEM and statistical modelling are provided online in the EIT-lab GitHub repository.

Model setup
The developed model was utilised to determine the optimal paradigm for dZ measurement at various distances from the stimulation site where CAPs are  [22]). For this, the FEM model of a C fibre ( figure 1(b)) was subjected to series of repetitive bipolar monophasic stimuli (50 µs pulse width, 5 nA) which were applied intracellularly across 0.1 mm at the end of the fibre. Two stimulation paradigms were used (figure 3(a)): (1) continuous stimulation with the frequency of 1 and 2 Hz (at these frequencies, the fibre never loses the ability to excitation); (2) stimulation with trains of stimuli of 5, 10, 20 and 50 Hz separated by resting intervals to allow the nerve recovery between the consecutive trains (table 1).
During repetitive stimulation at 5-50 Hz, the amplitudes of the consecutive APs were decreasing until the fibre lost the ability to excitation (figure 3(b)) that also happens experimentally due to the accumulation of potassium ions in the periaxonal space adjacent to the membrane [44,45]. Therefore, in case stimulating trains would last longer than the nerve is capable to be activated, the ratio of nerve firing time to the duration of the stimulation (duty cycle) would be reduced.
Thus, durations of the trains were chosen so that (1) the time when the nerve is in the active state (and hence the dZ) is maximised, that includes maximisation of the number of APs per train and the duty cycle, and (2) the nerve survives for the long term (>∼3 h in the saline bath) in the ex-vivo experiment (figure 2(a)), as increasing SNR to satisfying values may require prolonged averaging. The maximal durations allowing reaching the highest duty cycle were 8, 2.5, 0.75 and 0.25 s at 5, 10, 20 and 50 Hz respectively that was equal to 40, 24, 14 and 10 spikes per train at these frequencies (table 1). The resting time T rest between the 5-50 Hz trains was chosen to be 3 s (figure 3(a) and table 1).
In order to achieve long-term survival of the stimulated nerves so that they are stable during the experimental day period (CAPs are not changing for >6 h), these parameters were adjusted following testing performed in three nerves (table 3 in section 3). For this, CAP amplitudes in the nerves were extracted following their initial stimulation with the maximal train durations provided above and in table 1. Then, if the CAPs' amplitudes did not sustain for the duration of the experiment, the number of pulses and associated train durations were halved-this was repeated up to two times until the CAPs stable over several hours were achieved. To further improve stability, the resting time was also iteratively increased and its effect on the CAPs were observed.
In contrast to the APs, the amplitudes of the single-fibre impedance changes were increasing during stimulation ( figure 3(c)). The reasons for this behaviour are supposedly similar to the ones for the APs ( figure 3(b)): accumulation of potassium ions in the periaxonal space and sodium ions inside the fibre modifies reversal potentials of these ions and sensitises the associated ion channels thus leading to increased total conductance of the fibre during excitation. This effect is expected to improve the experimental dZ response during repetitive stimulation. However, it is hard to evaluate it experimentally since at least 1200 averages are required to reliably detect single spike dZ (10 min averaging at 2 Hz stimulation, section 2.2.2), and this would average out the progressively increasing dZ amplitudes. In addition, the number of spikes in the train and the resting time between trains was found to significantly affect nerve survival (6 pulses/train, 5 s between trains, table 3 in section 3)-it was vital for the purposes of the study, and the gradually increasing behaviour of dZ had therefore not been investigated further.
The obtained dZ trains (figure 3(c)) were incorporated into the statistical model for simulation of the dispersed compound dZ of a multi-fibre nerve described in the previous subsection. The SNR (ratio of mean signal to S.D. of the noise) at 3, 15, 20 and 50 cm from the onset of the stimulus was determined using two signal processing paradigms described in the next subsection. For averaging and noise reduction, the total duration of each simulation was chosen to be 30 min (table 1). 50 models were simulated in total to obtain statistics.

Signal processing
To extract compound impedance changes from the dispersed signals obtained in the statistical model, the following signal processing paradigms were performed; SNRs using these paradigms were obtained.  1(b)). Only a single train is shown at each frequency, although multiple trains were present at >10 Hz during the 10 s interval. In contrast to the APs, the amplitudes of consecutive dZ in the train are increasing. The shape of a single dZ is embedded into the last subplot.
(i) Averaging of single spikes (coherent spikes averaging). This approach has traditionally been used in previous studies involving dZ measurement [12,15,19,33]. In those studies, the recordings were cut into single spikes segments with the time window corresponding to the stimulation frequency around each spike ( figure 4(a)). Those segments were then (1) band-pass filtered using the bandwidths of 100-3000 Hz depending on the characteristic frequency of these signals (characteristic frequency of A fibres ≫ C fibres), (2) demodulated using the absolute of Hilbert transform (since the phase shift induced by the membrane is insignificant [34,36], and the dZ is approximately equal to |dZ|, equation (1) and (3) averaged together. In the current study, BW was chosen to be 200 Hz to account for the characteristic frequency of the simulated non-dispersed dZ ( figure 1). The resulting signals were averaged across all 50 computed models. SNR was calculated using the formula: where A signal is the maximal amplitude of the measured signal, σ noise is the standard deviation of the noise after filtering. (ii) Processing trains of spikes as a whole. Instead of cutting the recordings into single spikes, processing was conducted on the entire trains of spikes ( figure 4(b)). This allowed band-pass filtering with bandwidths around whole trains lasting up to a few seconds. As a result, the theoretical bandwidths for band-pass filtering could be significantly lower at 0.2-8 Hz, determined as BW 2 = 2 / T train where T train is the duration of the specific train (table 1). This could allow a significant reduction of the noise without longcontinued averaging. However, to differentiate the dZ signal of the C fibres from the signal of myelinated fibres, stimulation artefacts and low-frequency noise present in the recordings, these theoretical filtering bandwidths had to be increased (table 3 in section 3). SNR at 3, 15, 20 and 50 cm from the onset of the stimulus was computed according to (6), where A signal and σ noise corresponded to a new approach for signal processing. Due to the presence of stimulation artefacts (see the detailed description below) as well as fast myelinated fibres in the experimental recordings (figures 7 and 9), only the ending portions of the dispersed dZ corresponding to the slow C fibres could be recovered. Therefore, the last 100-500 ms in the processed dispersed dZ trains were considered for dZ and SNR calculation (figure 8 in the section 3). Although only the ending portions of the signals were used to measure dZ, the entire spike train signals were filtered: leaving the dZ in the central part of the signal allowed avoiding the filtering edge artefacts. Another artefact appearing due to presence of spikes at the end of the train when filtering the whole signal was of much shorter duration (in the order of tens of milliseconds) than the expected and observed dispersed dZ signal lasting up to 500 ms (figure 8), so they could be easily differentiated. Based on the obtained results, the optimal stimulation paradigm for recording dZ at far distances from the stimulus was determined.
Stimulation artefacts are inevitable when neural activity is evoked and recorded with the use of external electrodes [12,19,33]. Since the membrane of nerve fibres is more resistive than the surrounding connective tissue and physiological solution interface, part of the current applied during stimulation will flow through the conductive pathways along the nerve and an increase in the potential will be therefore measured by the recording electrodes ( figure 5). In addition, the voltage generated by the stimulation current is usually orders of magnitude higher than the recorded physiological voltages, so, the insufficient input range of the amplifier or small charge imbalance between stimulation electrodes may lead to saturation of the recording circuits which can last significantly longer than stimulation pulse itself [46,47]. To avoid saturation in the current study, charge-balanced biphasic stimulation was used together with an actiCHamp amplifier (Brainproducts GmbH, Gilching, Germany) having a wide input range of ±400 mV [16].
The influence of stimulation artefact can generally be minimised by placing the recording electrodes as far from the stimulating ones as possible. This is challenging in ex vivo conditions: surgically, it was not possible to extract intact undamaged SNs longer than ∼20 cm from the anesthetised pig in the current study. However, in vivo approach is much more promising as stimulation and recording can be done in different parts of the VNS (cervical and subdiaphragmatic); this is an essential part of the future work (see Discussion).
Stimulation artefacts obtained in the current study were nonsaturating [46] so that they could be minimised, and genuine CAP and dZ signals could be extracted through filtering. However, these artefacts could not be completely eliminated that severely affected signal processing in the case when the nerve was stimulated by trains of spikes ( figure 4(b)): only ending parts of the dispersed dZ signals had to be studied and the SNR could have therefore been reduced (see sections 3 and 4).

Experimental evaluation of the method using porcine SN ex vivo
All experimental procedures complied with regulations in the UK Animal (Scientific Procedures) Act, 1986 and were reviewed and approved by the Animal Welfare and Ethical Review Board. For experimental evaluation of the developed method, an ex-vivo setup with the SNs of the pig was used. In addition to porcine SNs being similar in fibre composition to the subdiaphragmatic branches of the VNSs in humans [31], they are very endurant-the stability of porcine SNs allowed to frequently manipulate them in the saline bath and conduct experiments lasting for more than 8 h, in accordance with survival times of other mammalian nerves [48,49]. This is a significant advantage over other unmyelinated nerves commonly used in ex-vivo setting such as a walking leg nerve of the crab [33,50] (figure A2 in supplementary material).
The experimental design was the same as in section 2, (b) describing the experimental adjustment of the model. Shortly, porcine SN were held in an organ bath perfusion chamber filled with continuously oxygenated saline solution. Three silicone rubber cuffs each having six radially arranged electrodes made from stainless steel and coated with PEDOT:pTS were placed around the nerve 3, 15 & 20 cm from the same cuff used for electrical stimulation (f stim = 2 Hz, I stim = 20-40 mA, PW = 50 µs, figure 2(a)). Impedance changes were measured by sequential application of the sinusoidal current through two last electrodes on each cuff, and the voltage was recorded on the remaining electrodes on the same cuff in respect to the last electrode on the last cuff ( figure 2(a)). Then, dZ was obtained by demodulation of the recorded voltage using the absolute of the Hilbert transform (1). Parameters of the applied sinusoidal current were: f AC = 1-6 kHz, I AC = 200-300 µA.
The optimal stimulation paradigm for recording dispersed dZ determined with the model was applied to N = 28 nerves to sequentially record dZ using cuffs 3 and 4 at 15 and 20 cm from the onset, where CAPs were dispersed and not measurable ( figure 2(a)). SNRs at these distances were obtained; implications for imaging unmyelinated nerves with fast neural EIT were investigated.
Statistical significance of the recorded dZ was verified using a two-sample t-test algorithm by comparison of the measurement amplitude straight after the stimulation artefacts (600-1000 ms, figure 8) with the dZ at all other points following this period (1000-3000 ms). The overall design of the study, including the workflow from development of the models to experimental verification, is presented in figure 5.
The MATLAB code written for data processing is provided in the EIT-lab GitHub repository; all the recorded unprocessed data will be available online in the SPARC portal (https://sparc.science/) after approval.

Experimental adjustment of the statistical model of the porcine SN
Initial recordings of CAPs and impedance changes at 3 cm from the stimulation site (cuff 2, N = 28 nerves, figure 2(a)) allowed adjusting the statistical model to correspond to experimental data (figure 6 and table 2). The CAPs recorded at cuff 2 were equal to 34 ± 17 µV for myelinated fibres and 69 ± 44 µV for unmyelinated fibres (mean ± s.d., figure 6(a)). The peak of myelinated fibres CAP was observed on average at 5 ms from the stimulus, peak of C fibres CAP-at 30 ms from the stimulus that allowed to determine their CVs: v fast = 8 ± 3 m s −1 , v C = 0.8 ± 0.3 m s −1 (mean ± s.d., table 2). The amplitude of the CAP of fast myelinated fibres was lower in the experiment than in the model because part of it was covered by stimulation artefact thus decreasing its amplitude ( figure 6(a)). The dZ were equal to (1.42 ± 1.11) × 10 −4 % (0.16 ± 0.19 µV) and (6.96 ± 4.61) × 10 −4 % (0.74 ± 0.58 µV) in fast myelinated fibres and unmyelinated C fibres respectively.

Method development and optimisation
Using the designed statistical model of dispersion in the porcine SN, the method for recording dZ further from the site of the onset than CAPs are measurable was developed and optimised.
Images of the dispersed dZ with and without fast myelinated fibres (the real porcine SN case and artificial unmyelinated case) show that the C fibre activity is more dispersed and constitute a low-frequency component in the resultant signal, while myelinated fibres correspond to a high-frequency component (figure 7, blue and red lines). Therefore, to measure dZ of the C fibres in the mixed porcine SN nerve, the ending parts of the signals were considered.
By stimulation of three nerves with trains of pulses of 5, 10, 20 and 50 Hz with 3 s intervals between trains, and following the approach described in section 2.3.1, it was found that to achieve long-term (>6 h) survival of the nerves, the number of pulses in each train must not exceed six and the interval between trains must be increased to 5 s, independent of the applied train frequency. The model was modified accordingly (table 3).
The developed model has shown that the SNR obtained using the coherent spike averaging approach for signal processing (table 1 and figure 4(a)) were larger at short distances from the onset while falling exponentially at longer distances (solid lines in figure 8 and table 4). The second approach involving processing of whole trains (table 1 and figure 4(b)) led to smaller SNR at 3 cm from the stimulus, but increased to >1 at 15 and 20 cm from the site of the stimulus after 30 min of averaging (figure 8 and  table 4). Therefore, the optimal parameters for measurement of dZ at far distances from the site stimulation predicted with the model were: stimulation with 5 or 10 Hz trains, six pulses per train, 5 s interval between trains (table 4). Filtering bandwidth was increased to 10 Hz to differentiate the dZ signal of the C fibres from the signal corresponding to myelinated fibres, stimulation artefacts and low-frequency noise present in the recordings (figures 7 and 9).
Even with the optimal parameters, averaging for 30 min only produces SNR marginally higher than the limit of detectability at 15 and 20 cm from the stimulus (1.8 ± 0.8 at 10 Hz, 15 cm, table 4). Therefore, to measure dZ at longer distances as well as to obtain an SNR of 4 which is minimally required for reproducible imaging of fast impedance changes with EIT [18], longer averaging will be required (see Discussion for details).

Experimental evaluation of the developed approach
The optimal stimulation and signal processing paradigm for recording dispersed dZ in unmyelinated fibres determined in the modelling study (10 Hz trains, 6 pulses/train, 5 s interval between trains, filtering bandwidth 10 Hz) were applied to 28 porcine SNs ex-vivo, out of which N = 18 recordings had a satisfactory level of noise smaller than 4 µV RMS before averaging, in agreement with the developed modified model.
For each nerve, the dZ measurement was sequentially performed using cuff 3 and cuff 4, located at 15 and 20 cm from the site of stimulation respectively (figure 2). Due to the presence of fast myelinated fibres and stimulation artefacts in the recordings (figures 7 and 9), measurements were done at the ending stages of the dispersed dZ, from 0.6 to 1.1 s and from 0.65 to 1.15 s for the used 10 Hz trains at cuff 3 and cuff 4 respectively ( figure 7(b)).
The resulting dZ were equal to (1.11 ± 1.03) × 10 −4 % (0.11 ± 0.10 µV) at 15 cm (cuff 3, figure 2(a)), and (1.17 ± 1.21) × 10 −4 % (0.12 ± 0.10 µV) at 20 cm (cuff 4, figure 2(a)) from the site of stimulation. The SNR at 15 cm after 30 min of averaging was 1.8 ± 0.7, decreasing to 1.7 ± 0.6 at 20 cm that is in agreement with the predictions of the developed model ( figure 8(b) and table 4). The mean absolute value of the determined dZ across all nerves was found to be significantly larger than the mean at every other point in the recording (P < 0.01, N = 18).     The stimulation frequencies where SNR ⩾ 1 and the optimal train frequencies are highlighted in blue (b) The designed model was used for the development and optimisation of a novel stimulation and signal processing paradigm to record impedance changes in unmyelinated nerves at the distances from the stimulus where CAPs are dispersed ( figure 5, step 4). This determined optimal paradigm was the following: 5 or 10 Hz trains with 6 pulses/train and 5 s interval between trains, and subsequent signal processing using 10 Hz band-pass filter to account for the presence of myelinated fibres, low-frequency noise and stimulation artefacts (  8 and table 4). Thus, more than 30 min of averaging would be required to record dZ further than 20 cm from the stimulus. (c) The model's predictions were evaluated experimentally by stimulation of pig SNs by trains of stimuli ex vivo (figures 9, 5 and step 5). The predictions were in good agreement with the simulations: although the levels of noise were significantly reduced compared to the previous dZ measurement experiments [12], SNR obtained after 30 min of averaging (section 3) was at the edge of detectability at 15 and 20 cm from the onset of the stimulus. This was partly due to low-frequency noise and stimulation artefacts observed in the recordings that required increasing the bandwidth of filtering to a minimum of 10 Hz.

Answers to the stated questions
(a) What are the optimal stimulation and signal processing strategies producing the largest impedance changes at different distances (15, 20 and 50 cm) from the onset?
The optimal stimulation strategies producing the largest impedance change signal are stimulation with 10 Hz trains, with 6 pulses per train and 5 s between trains, for the long-term survival of the nerve necessary for the desired chronic implantation. Other parameters were: AC current-200-300 µA amplitude, 1-1.5 kHz frequency, stimulation current-20-30 mA amplitude, 50 µs pulse width.
(b) How much averaging is required (1) to obtain a measurable signal (SNR > 1) and (2) to image neural activity with EIT (which requires SNR > 4 [18]) at 15, 20 and 50 cm from the site of stimulation?
Given the very low experimentally achieved noise of 3.5 µV RMS before averaging, the dZ in the porcine SN at 15-20 cm from the stimulus was at the edge of detectability (SNR < 2) if averaged for 30 min. To obtain larger SNR, longer averaging would be required. At 50 cm from the stimulus, to obtain SNR ⩾ 2 one would need to average for approximately 1.5 h, given the noise decreases with the square root of the number of averages (table 4).
In order to image fast neural activity in fascicles of unmyelinated nerves in their cross-section, a minimal SNR of 4 obtained by sequential switching between 14 electrode pairs is required [18,19]. Therefore, to achieve imaging with the obtained optimal paradigm, and considering that the noise decreases with the square root of the duration of averaging, the minimal averaging time would equal to 30 min · (4/1.8) 2 · 14 electrode pairs ≈ 34 h at 15 cm from the site of stimulation. This time would increase significantly for imaging at further distances in accordance with the table 4. However, even such a long duration of averaging may be clinically feasible, for example, in implantable neuromodulation devices [2][3][4] that may run the specified stimulation paradigm for days or even weeks without any adverse effects to patients. In addition, if the internal organisation of the nerve being stimulated is known, for example, using the feedback system based on the physiological responses (such as respiratory breath rate or heart rate) following stimulation with multi-electrode cuff arrays [51], there is no need in the execution of a full imaging paradigm with 14 electrode pairs. Instead, one could concentrate on specific electrode pairs corresponding to the required fascicles that may help to decrease the duration of averaging by two-or three-fold.
In addition, future plans include transitioning to an in vivo experimental paradigm which has several potential advantages over the ex vivo one used in the  figure 4) of the selected experimental recordings performed with the optimal stimulation paradigm obtained in the modelling study (10 Hz trains, BW = 10 Hz). Noisy spikes in the beginnings of the recordings are stimulation artefacts. Locations where dZ was expected and where it was measured are highlighted in red. current study. Although the level of noise is expected to increase in vivo (to 8-10 µV before averaging [12,19]), the SNR may be improved (given enough averaging) for the following reasons: • The contact between the cuff and the nerve should improve due to the absence of saline solution in and around the cuffs, and this will lead to stronger signals. • The whole length of the VNS can be utilised in vivo that will significantly reduce stimulation artefacts strongly affecting the recorded dZ signals.
In the anaesthetized animal, stimulation can be done in the cervical part of the nerve while the recording cuff can be placed in the subdiaphragmatic part. This was not possible to perform ex vivo due to the limited size of the nerve bath and complex surgical procedures required to remove longer parts of the nerve from the body. Removal of stimulation artefacts will enable measurement of the whole duration of the compound dZ signals, and not just the ending parts which follow stimulation artefacts (figure 9), and this is expected to strongly improve the SNR.
Even with such a long averaging expected to be required for EIT imaging, the new experimental paradigm developed in the study demonstrates the feasibility of recording impedance changes far from the site of stimulation, which cannot be currently achieved using other existing methods.
(c) Are simulated results confirmed with experimental data?
The SNR measured in the performed ex-vivo experiments were in agreement with the values obtained with the developed model of dispersion (figures 8, 9 and section 3)). With the average level of noise achieved in the developed experimental setup (3.5 µV RMS), the maximal distance from the onset where significant reproducible dZ could be measured in 30 min was 15 cm.

Limitations
The current study had several limitations. First, it was the statistical model and not the FEM model which was used for simulation of dispersion in unmyelinated nerves and for obtaining the optimal parameters for dZ measurement. Since computational time and the amount of required resources rise exponentially with the number of fibres in the FEM model, FEM simulation of 40 thousand-fibre nerve would be very computationally intensive and will demand unrealistic time and resources. For instance, it takes about an hour for a 40 ms simulation of the single C fibre FEM model described by a system of 20 differential equations, while a model with 50 fibres of the same type required a whole week on a double-CPU machine with 128 Gb RAM [22]. However, the statistical model was based on the signals simulated with the realistic mammalian C fibre FEM model that brought them into accordance with real experimental data.
Second, the FEM model was not solved for each fibre location in the cross-section of the cylindrical external space, so that the simulations were performed only for the fibre placed in the centre. This was done for two main reasons. First, the system was symmetric (1 µm fibre in the centre of 10 µm external space), so that placement of the fibre in other positions inside the cylinder could be assumed to not significantly alter the shape of the measured signal but to mainly affect its amplitude. Also, during the transition to the multi-fibre statistical model, the diameter of the cuff was re-adjusted so that the uniformly distributed fibres were on average equidistant from the electrodes (supplementary material, equation (3A)), and the error brought by small variations in the shapes of single APs is not expected to significantly contribute into the error of the resultant CAP. In addition, this simplification allowed to significantly reduce computational time: in the current study, the fibre was stimulated with trains of pulses lasting up to 10 s (table 1) that took up to 40 h per single simulation.
Third, the porcine SNs used in the ex-vivo experimental study were short, of up to 20 cm long, that led to large stimulation artefacts overlapping with the dZ measured at these distances. In addition, since around 10% of the porcine SN consists of fast and large myelinated fibres producing larger dZ, they produce the artefactual signal of the same kind as stimulation appearing during the initial phase of the dispersed signal ( figure 7). Therefore, only the ending part of the dispersed dZ not covered under these artefacts had to be studied (figure 9), and the SNR may have thus been reduced. One way to significantly reduce these artefacts is to use longer nerves, and this can be achieved in the in vivo experiments which are planned as the next step of the presented study.

Future work
In order to remove or significantly reduce the stimulation artefacts and the artefacts caused by the presence of the myelinated fibres, the future work will be to perform an experimental study in paralysed pigs in vivo with stimulation on the cervical part of the VNS and recording on the subdiaphragmatic part which is around one metre apart. This will allow measurement of the pure dZ signal not contaminated by this type of artefact. In addition, in vivo setup may provide better contact between the cuff electrodes and the nerve due to the absence of a highly conductive saline interface present in the ex vivo bath; this will potentially lead to larger signals and the SNR.

Conclusion
It is challenging to measure CAPs and impedance changes in unmyelinated nerves starting from a few centimetres from the site of stimulation. The developed experimentally adjusted computational model of TD in nerve allowed simulation of compound APs and dZs at various distances from the variable stimulus. With the model, optimal stimulation and signal processing parameters for dZ measurement were determined. It was shown that stimulation of the nerve by trains of stimuli allows recording impedance changes further from the onset than it is possible with a traditional continuous stimulation and averaging of consecutive spikes. The findings were evaluated experimentally using the porcine SN ex vivo. This work enables a new way for measurement of impedance changes accompanying excitation at distances from the stimulation where standard approaches are not feasible.
The models of dispersion in complex-fibre nerves designed for optimisation of the experimental parameters in the current study can be used by research community for studying the dispersion-related properties in any types of nerves and the development of novel techniques for sensing neural activity in unmyelinated nerves. The latter can be particularly useful for facilitation of the novel field of bioelectronic medicines aimed at neuromodulation of internal organs via stimulation of the VNS: development of novel closed-loop solutions involving stimulation of the nerve in response to the recorded physiological signals will enhance the treatment outcomes for a variate of drug-resistant disorders including, among others, epilepsy, depression, and cardiovascular diseases.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.
The COMSOL and MATLAB model files used for FEM and statistical modelling, as well as the MATLAB code written for data processing are provided online in the EIT-lab GitHub repository [43]. All the recorded unprocessed data will be available online in the SPARC portal (https://sparc.science/) after approval.