Effect of dispersion in nerve on compound action potential and impedance change: a modelling study

Objective: Electrical impedance tomography (EIT) is capable of imaging fast compound electrical activity (compound action potentials, or CAPs) inside peripheral nerves. The ability of EIT to detect impedance changes (dZ) which arise from the opening of ion channels during the CAP is limited by the dispersion with distance from the site of onset, as fibres have differing conduction velocities. The effect is largest for autonomic nerves mainly formed of slower conducting unmyelinated fibres where signals cannot be recorded more than a few cm away from the stimulation. However, as CAPs are biphasic, monophasic dZ are expected to be detectable further than them; testing this hypothesis was the main objective of this study. Approach: An anatomically accurate FEM model and simplified statistical models of 50-fibre Hodgkin–Huxley and C-nociceptor nerves were developed with normally distributed conduction velocities; the statistical models were extended to realistic nerves. Main results: Fifty-fibre models showed that dZ could persist further than biphasic CAPs, as these then cancelled. For realistic nerves consisting of Aα or Aβ fibres, significant dZ could be detected at 50 cm from the onset site with signal-to-noise ratios (SNR, mean  ±  s.d.) of 2.7  ±  0.2 and 1.8  ±  0.1 respectively; Aδ and rat sciatic nerve—at 20 cm (1.6  ±  0.03 and 1.6  ±  0.06), rat vagus—at 10 cm (1.6  ±  0.05); C fibres—at 1–2 cm (2.4  ±  0.02). Significance: This study provides a basis for determining the distance over which EIT may be used to image fascicular activity in electroceuticals and suggests dZ will persist further than CAPs if biphasic.


Introduction
Electrical impedance tomography (EIT) allows imaging of impedance changes (dZ) in neural tissue which arise from the opening of ion channels over milliseconds (Oh et al 2011, Aristovich et al 2016, Fouchard et al 2016, Faulkner et al 2017. It has the potential to be used in neuroprosthetics (Hope et al 2018) and the new area of electroceuticals (Famm 2013, Waltz 2016. This latter aims to treat disease by electrical stimulation of autonomic nerves. EIT with a cylindrical nerve cuff offers a novel means to image the compound action potential (CAPs) in fascicles within autonomic nerves and so provides a method to avoid offtarget effects by selective stimulation of identified fascicles.
In Electroceuticals in humans, one aim might be to undertake fast neural EIT with a cuff around the cervical vagus nerve and modulate activity in abdominal organs which could be one metre distant. Fast neural EIT relies on imaging of impedance changes over time and so the differences over milliseconds during the CAP. Dispersion has the effect of smoothing out these differences and so poses a challenge for recording EIT images in relation to the CAP in electroceuticals.

Dispersion in nerves
Due to variability in sizes of fibres in nerves and proportionality of conduction velocity to the fibre diameter (Waxman 1980), propagation velocities of individual fibres vary. As a result, the amplitude of the compound action potential being an aggregate sum of all individual action potentials decreases along the nerve following 2 I Tarotin et al its initiation; this is commonly defined as temporal dispersion (Freeman 1972, Dorfman 1984, Olney et al 1987, Taylor 1993, Schulte-Mattler et al 2001 The effect of dispersion is much greater in small diameter unmyelinated nerves. It has mainly been studied in peripheral nerves consisting of mainly large myelinated fibres because they produce signals of higher magnitudes and possess lower stimulation thresholds Torebjörk 1973, Torebjörk andHallin 1974). For example, there was a 36% reduction in the compound action potential when stimulating the human ulnar nerve at aboveelbow and wrist regions and recording at the fifth digit (Olney et al 1987). In Taylor (1993) and Schulte-Mattler et al (2001), large human peripheral nerves including median, ulnar, common peroneal and tibial were stimulated and compound motor action potentials were recorded a few tens of cm away from the stimulus. These studies showed only a small amplitude and area decay which were in the range of 5%-45% per meter of nerve length. In contrast, in the olfactory nerve of the cat (Freeman 1972) which comprises mainly small unmyelinated fibres 0.1-0.5 μm in diameter, CAPs could not be recorded further than 2.5 mm from the stimulation site. In vagal C fibres, there was a >50% decrease at 4 mm from the point of activation (Chang et al 2015).

Fast neural EIT
Similar to traditional EIT, fast neural EIT is able to reconstruct an image of change in conductance of a tissue by injection of electrical currents and recording surface voltages . Fast neural EIT has been shown to be capable of imaging fast changes in the somatosensory cortex of the rat (Aristovich et al 2016) and was successfully applied to rat sciatic nerve in vivo to obtain real-time images in its tibial and peroneal fascicles . For this, the posterior tibial and common peroneal nerves were electrically stimulated and the activity was recorded with the flexible multi-electrode cuff placed on the main sciatic nerve trunk few centimetres away from the site of stimulation. The spatial and temporal resolutions of the images were 100 μm and 0.3 ms respectively which were shown to be better than those obtained using inverse-source analysis. The results presented in that study served as a proof of concept and demonstrated the high potential of fast neural EIT to image inside nerves at a fascicular level.
However, about half of the fibres in the rat sciatic nerve are myelinated (Schmalbruch 1986, Fouchard et al 2016 with larger diameters and conduction velocities and therefore they occupy the majority of the nerve crosssectional area. Only these large fibres were imaged in Aristovich et al (2018) because they had lower activation thresholds and significantly lower levels of dispersion in comparison to the unmyelinated ones (for detailed explanation see section 1.3 below).

Implications of dispersion for fast neural EIT
Using EIT for imaging inside autonomic nerves is more challenging as there is a lower signal-to-noise ratio. This has motivated the development of numeric models of ion channel opening and the consequent impedance changes in peripheral nerve with the aim of suggesting optimal electrical and geometrical parameters for EIT during the CAP. These have so far been published for unmyelinated giant squid axons with Hodgkin-Huxley (HH) ion channels (Hodgkin and Huxley 1952) and mammalian C-nociceptors (Tigerholm et al 2014) with up to eight identical fibres with simultaneously initiated action potentials (AP) (Tarotin et al 2019). However, real autonomic nerves such as the cervical vagus nerve comprise thousands of similar small fibres. More than twothirds of those are unmyelinated in mammals (De Neef et al 1982, Asala and Bower 1986, Prechtl and Powley 1987, Soltanpour and Santer 1996 and humans (Shimizu et al 2011, Verlinden et al 2016; these have slow conduction velocities (CV) of about ~0.5-2 m s −1 (Coleridge and Coleridge 1984). These fibres have differing diameters and so CV. This causes dispersion of action potentials along the nerve during propagation. Consequently, it is challenging to record a signal of sufficient amplitude when moving away from a point of stimulation, and this effect is proportionately greater for unmyelinated fibres than faster conducting fibres. For example, in rat sciatic nerve, significant impedance changes could only be recorded in fast fibres with low dispersion more than 5 cm from the onset site . Impedance measurements in the walking leg nerve of the crab which are unmyelinated (Boone 1995) also demonstrated high dispersion with the possibility of CAP recording only up to 1.6 cm from the stimulus. In addition, only fast fibres were possible to be imaged with fast neural EIT in the rat sciatic nerve experiment  while the activity of C fibres was not visible.
In spite of this issue, there is potential for EIT to image unmyelinated fibres at greater distances than the CAP. This is because extracellular action potentials of individual fibres may be expected to disperse sooner than the related impedance changes. This is because CAPs are usually biphasic or triphasic (Harper and Lawson 1985, Gold et al 2006, Agudelo-Toro and Neef 2013, Ghitani et al 2017 so they cancel when separated by dispersion. This effect may be expected to be much less for impedance changes which are mainly monophasic (Faulkner et al 2017; in principle, therefore, dZ should decrease slower than CAP. We have explored the possibility of recording dZ further than CAPs by developing 3D FEM computational models comprising 50 fibres with HH (Hodgkin and Huxley 1952) or C-nociceptor (Tigerholm et al 2014) ion channels with normally distributed sizes and propagation velocities. These models were bi-directionally coupled with the external space which allowed injection of electric current via external electrodes and simultaneous external recording at various distances down the nerve. Simplified statistical models matching the accurate ones were developed to accelerate computations and allow dZ simulations of complex nerves with >10 k axons. If computationally optimised, the developed accurate 3D FEM models could be extended to up to thousands of fibres with various geometrical and electrical properties and with the addition of connecting tissue. Such full 3D models would be a valuable tool for simulation of nerve behaviour in normal or abnormal conditions and under different external stimuli.

Purpose
The purpose of this study was to investigate the impedance change behaviour in complex compound nerves relative to CAP. A specific interest was to test the hypothesis that impedance changes can be recorded further down the nerve from the site of its stimulation than the CAP.
Specific questions addressed in this study were: 1. Is it possible to record dZ further than CAP? 2. Are there differences between models in the effect of dispersion on dZ? 3. What is the largest distance from the site of stimulation at which dZ can be recorded?
a. For nerves consisting of Aα, Aβ, Aδ and C fibres; b. for realistic nerves such as the right vagus and sciatic nerves of the rat.

Experimental design
The study was accomplished by development of 3D fully coupled models of nerve comprising 50 active fibres with normally distributed conduction velocities: giant axons of the squid (Hodgkin and Huxley 1952) or mammalian C-nociceptors (Tigerholm et al 2014). These models were based on previously developed accurate active models of single fibres coupled with external space (Tarotin et al 2019). They utilized the finite element method (FEM) approach and were built in COMSOL Multiphysics software (COMSOL Inc., USA) in conjunction with MATLAB using LiveLink for MATLAB interface. Due to the different size of the axons involved, diameters of the nerve models and circumferential electrode rings were 2.4 and 0.01 cm for the HH and C fibre models, respectively. The width of the electrode rings was adjusted accordingly to 100 and 5 μm respectively (table 1). Although a nerve diameter of 2.4 cm nerve exceeds the human anatomical range, it needed to be this large to contain 50 giant axons of the squid, each having 1 mm diameter. Such a nerve model is computationally efficient and is a useful basis to answer the questions about dispersion stated above. Realistic EIT experimental conditions were simulated in which a linear electrode array was used for dZ measurements (Oh et al 2011, Aristovich et al 2015. This enabled calculation of action potentials and impedance changes by external electrodes along the nerve in order to explore if the initial hypothesis were supported. The detailed models were complemented with simplified statistical models which allowed dZ of much larger nerves with >10 k axons to be simulated without extensive computation; these were still in good agreement with the 3D FEM models. Statistical models were based on the histological studies which provided morphological properties of fibres and their conduction velocities. Measurements with FEM models were accomplished using direct current; alternating current was added in the statistical models. The reason was that simulations at DC demanded significantly smaller time to compute: 8 h per single simulation versus 20 h for the 225 Hz AC case; this time increased with increasing frequency because of the appearance of transient effects and their interaction with the active membrane. Modelling (Tarotin et al 2019) showed that the dZ reached maximum at DC and no differences in biophysical origin of dZ recorded using AC or DC were observed.

Model setup 2.2.1. Accurate FEM models
The current models included 50 1D active nerve fibres and a 3D extracellular space bi-directionally coupled together. The developed coupling feedback system allowed simultaneous recording of intracellular action potentials and extracellular field created by the membrane currents and the injected current (Tarotin et al 2019). The fibres were simulated as cables with active voltage-dependent experimentally validated ion channels using either Hodgkin-Huxley model of the giant axon of the squid or mammalian C fibre with 10 active ion channels and voltage-dependent ions' concentrations (Tigerholm et al 2014).
The interaction between fibres was omitted as it was shown to not significantly affect the measurements (Tarotin et al 2019). The external space was represented as a cylinder with a uniform electrical conductivity of the extracellular medium of 10 mS cm −1 (Elia and Lamberti 2013). EIT current (table 1) was injected and the result-ant voltage was recorded by means of external electrodes. The mesh of HH and C fibre models comprised 1.6 M and 530 k tetrahedral elements, respectively.
The main equations and electrical parameters of the models can be found in Tarotin et al (2019) and in the appendix. The external field was simulated with volume conduction Poisson's equation; the current was injected as a constant flux through the ring electrodes. Stimulation of the fibres was implemented via the inclusion of an activating function into the main equations of the nerve fibres. Membrane currents flowing out of the fibres to the external space were represented as a normal flux through the boundaries of the modelled axons.
Fibres were uniformly distributed inside the insulated cylinder with the diameter depending on the type of fibres (table 1). The size of the cylinder defined the diameters of the electrodes as well as the amount of the external space (table 1, figure 1). Impedance of the electrodes did not affect the measured dZ because a 4-electrode measurement paradigm (Schwan and Ferris 1968) was used for time-difference single impedance measurements (Holder 2004).
Dispersion was implemented by normally distributing propagation velocities of the fibres by randomizing their intracellular resistivities. Their mean values were taken from the developed validated models: 0.05 kΩ·cm and 0.02 kΩ·cm for the model with HH axons, and 0.0354 and 0.01 kΩ·cm for the one with C fibres (table 1) (Hodgkin andHuxley 1952, Tigerholm et al 2014). Standard deviations were chosen to have a significant visible dispersion along the modelled lengths of the nerves (60 cm for HH and 2 cm for C fibre, table 1) and so that it is similar to experimental recordings (Gasser 1950, Boyd andKalu 1979) (table 2).
For each model, four sets of electrodes were placed along the nerve; each set consisted of two injecting ones through which the constant small direct current not affecting the membrane was applied: 30 μA (40 μA cm −2 ) for HH and 6.3 nA (4 mA cm −2 ) for C-fibre nerves; the recording electrode used for measuring the resultant external voltage in respect to ground (figure 1). Impedance changes recorded at DC have been shown to have the largest values and the same physiological nature as the AC ones (Tarotin et al 2019), so they provided a good starting point for simulations. Electrodes had diameters of 2.4 cm and 0.01 cm and widths of 0.1 cm and 5 μm for models with HH and C axons respectively (table 1). This large electrode diameter for the HH axons was chosen to accommodate all giant squid axons (d = 1 mm) inside it and so that the packing density would be similar to the nerve with C fibres (figure 1).
The location of the recording electrode in each set was at 10, 19, 25, 35 cm from the activation point for the HH nerve and 0.4, 0.8, 1, 1.4 cm for C fibre nerve; distances between recording and injecting electrodes (Δx R ) for HH and C fibre cases were 0.2 and 0.01 cm; distances between injecting electrodes (Δx R ) were 0.002 and 0.5 cm respectively (table 1). The action potentials were initiated simultaneously at the end of the nerve.
The COMSOL model files saved in Matlab format are provided in the EIT-lab GitHub repository online as well as on the Physiological Measurement website 1 .

Signal processing
Two simulations with and without current injection were done for each location of the electrode sets. Each simulation lasted 40 ms for HH and 30 ms for C fibre models. Overall, 4 electrodes locations · 2 states = 8 simulations were done for each type of the nerve. The process of dZ extraction was identical to the ones utilized in recent studies , Tarotin et al 2019. In case of DC injection which was performed in the current study (figure 2), it included (a) subtraction of the signals with and without the injected current to eliminate action potentials but preserve impedance change and (b) subtraction of the mean baseline voltage in order to express the impedance change in terms of absolute voltage change in μV. If percentage values were needed, normalization of the obtained absolute dZ was done. The dZ could be expressed in μV because the constant current was injected and the phase shift between the current and the measured voltage is close to zero (Cole and Curtis 1939). Thus, voltage traces in figures 2, 5 and 6 are referred to as dZ as they signify impedance changes and can be directly compared with dZ expressed in percent (Tarotin et al 2019).
The negative integral areas under the CAPs and dZ curves recorded by each electrode set were calculated to compare the effect of their cancellation due to dispersion. Because the absolute values of these areas had different orders of magnitude and in order to compare their behaviour, they were normalized in respect to the ones at the shortest computed distances (table 1): 10 cm for HH and 0.4 cm for C fibre models.  (table 2). The AP was induced at the end simultaneously for all the fibres; the current was injected via two external electrodes (blue) and the electric field was recorded by an external electrode (green) located before the injecting ones with respect to ground. Four sets of electrodes were used to record impedance change accompanying action potential propagation (table 1). (b) Example of tetrahedral FEM mesh used in the models (section 2.1). (c) Side view of the uniform distribution of the fibres used in the models. Table 2. Parameters of the statistical models (Gasser 1950, Rushton 1951, Birren and Wall 1956, Boyd and Kalu 1979 3). Single fibre dZ signals measured at DC and at AC currents (1 kHz for HH and 2 kHz for C fibre) were used for simulations. To increase the absolute magnitude of impedance changes measured at AC, the EIT current used for obtaining them was 10 times larger than at DC, due to the increase in safe range values (Tarotin et al 2019). Because dZ is linear with respect to the amplitude of the injected EIT current (Holder 2004, the absolute μV amplitudes of the single fibre signals were linearly corrected to the same current level as in the FEM model. Before summation, these signals were scaled in amplitude due to the difference in electrode diameters used for single fibre and 50-fibre simulations and therefore different distances from the electrode to the fibre. The recorded amplitude reduction due to the increase in distance to the fibre was proportional to 1/r 2 , according to Coulomb's law, therefore the used single fibre signals were decreased (d 50 /d 1 ) 2 times where d 50 is the diameter of the electrodes used in 50-fibre models (2.4 and 0.01 cm, table 1) and d 1 are the diameters used for single fibre models which equalled 0.6 cm and 10 μm for HH and C fibre model respectively (Tarotin et al 2019). Thus, the single fibre signals inserted into the developed simplified 50-fibre model were decreased 100 times in the case of C fibres and 16 times for HH axons (table 1). The parameters of the velocities' distribution were chosen so that they correspond to the distribution of axoplasm resistivities in the FEM model (table 1). As the conduction velocity of a nerve fibre equals the ratio of its length constant to the time constant, it can be written as follows (Johnston and Wu 1995): where λ (m) and τ (s) are length and time constants of the nerve, r ax (cm) is the axon radius, ρ m (Ω · cm), ρ i (Ω · cm 2 ) and C m (F cm −2 ) are a membrane and intracellular resistivities and membrane capacitance. Thus, CV is inversely proportional to the square root of intracellular resistivity. Given that mean values are 15 m s −1 for HH model and ~0.6 m s −1 for C fibre model, S.D. obtained from (1) were 2.8 and 0.07 m s −1 for these fibre types respectively (table 1).
The source code for 50-fibre statistical models is available online and on the website of Physiological Measurement 1 .

Statistical modelling of dZ in large complex nerves
As it was impossible to implement computationally-heavy realistic 3D FEM models for realistic nerves containing thousands of fibres, simplified statistical models containing up to 100 k different fibre types were developed. The implemented approach was based on the inverse proportionality of the recorded impedance changes to the cross-sectional area of the measured nerves; the decrease in fast neural impedance change is observed because the intracellular resistivity of fibres is lower than the external one. When ion channels open, the injected current flows to the area with lower resistivity and we observe a small decrease in the dZ of the system 'nerve + external space'. Approximating the nerve fibre as a cylinder with a constant diameter, its resistance can be written as: where ρ is the axoplasm resistivity (Ω · cm), L n (cm), S n (cm 2 ) and d (cm) are length, cross-sectional area and diameter of the nerve. So, the resistance of the nerve is inversely proportional to cross-sectional area meaning that the recorded dZ has the same dependence. Simple statistical models were developed to find the maximal distances at which impedance changes may be theoretically measured. This was accomplished for four types of nerves consisting of a single type of fibres: Aα, Figure 2. Signal processing utilized to extract impedance changes (dZ) from the recorded voltages, similar to Tarotin et al (2019). The AP signal without the injected current was subtracted from the one with the injected current; the mean baseline voltage was subtracted to express the signal as a change of voltage; if percentage values were needed, the resultant signal was normalized.
Aβ, Aδ and C as well as for realistic sciatic and right vagus nerves of the rat (table 2). Compound action potentials were not simulated in these models as the main interest was to study how dZ disperses with the distance from the onset site.
These models were implemented in the same way as the simplified 50-fibre models. In the models, the compound impedance changes at DC and at 2 kHz were computed as a sum of modified dZs of all single fibres at different distances from the stimulation point (figure 3). Conduction velocities in each fibre were assumed to be constant and normally distributed; the values for normal distributions were taken from the exper imental data found in the literature (tables 2 and 3). Modifications of the previously simulated dZ of a single realistic mammalian C fibre depended on the fibre type (figure 3): the dZ latency for fast fibres linearly decreased in accordance to the latency of their APs (table 2); the dZ amplitude increased proportionally to their cross-sectional area (2) (figure 3). Conduction velocity (mean and S.D.) and latency data for Aα, Aβ and Aδ fibres was taken from Boyd . Schematic representation of impedance change dispersion in the developed model. The same is applicable to action potentials. The models in this study were implemented in three steps: (1) modification of the dZs of each fibre in accordance with their cross-sectional area (2) and action potential latencies; (2) summarizing the AP or dZ of single fibres having constant but different normally distributed conduction velocities; (3) scaling of the resultant signals in accordance with the nerve diameter and conductivity of connective tissue (4) and adding the experimental noise.
Step 2 on this picture demonstrates on a simple three-fibre example, how the compound dZ (red) is formed from the dZ of single fibres (blue) with slightly different velocities. In each fibre, the action potential propagates with constant velocity V i ; the compound dZ is equal to the sum of dZ i in all single fibres. It is seen that the amplitude of compound dZ decreases with distance. The shape of dZ used this picture was obtained from simulations for realistic single C fibre (Tarotin et al 2019); it was inserted into the developed statistical multiple-fibre models with velocity distribution taken from the literature (table 2). and Kalu (1979) for motor and cutaneous nerves of the cat. Data for C fibres was obtained from diameters distribution (Gasser 1950) and experimentally found equation (3a) (Gasser 1950, Rushton 1951 where v unmyel and v myel (m s −1 ) are conduction velocities for unmyelinated and myelinated fibres and d fibre (μm 2 ) is the cross-sectional area of the fibre. A number of fibres in each type of nerve was chosen so that they fit into the 0.1 mm 2 rat vagus nerve (Soltanpour and Santer 1996) given that around 60% of nerve cross-sectional area is occupied by fibres (Birren and Wall 1956) (table 2). The numbers were ( figure 6): 300, 500, 4800 and 78 000 fibres for Aα, Aβ, Aδ and C fibres respectively.
All the parameters for realistic nerves (table 3) were based on morphometric data for sciatic (Schmalbruch 1986) and vagus nerves of the rat (Soltanpour and Santer 1996) and the experimentally found CV versus diameter dependencies for unmyelinated C fibre (Gasser 1950, Rushton 1951, as well as myelinated Aβ and small Aδ fibres (3b) (Boyd and Kalu 1979). Each nerve type was simulated with up to 50 cm length and 50 ms duration.
Because some of the parameters in the model of a single C fibre (Tarotin et al 2019) differed from the realistic nerve model, the dZ taken from it had to be scaled before fitting it into the models with multiple axons. There were two main differences: first, extracellular space in a single fibre model was simulated as a high conductive saline, while in reality it consists of more resistive connective tissue; second, the diameter of the recording electrode was small (10 μm) so that it could fit only one fibre. Thus, the scaling coefficients k conn and k el were introduced to account for these discrepancies; the dZ obtained after summation of all single dZ across all fibres (figure 3) was scaled by dividing it by the product of these coefficients (4). The connective tissue was assumed to consist of 3% endoneurium, 3% perineurium and 94% epineurium with resistivities 1.211, 1.211 and 47.8 kΩ · cm respectively (Choi et al 2001). The weighted average of these values equalling 2.6 kΩ · cm was used which was around k conn = 26 times smaller than 0.1 kΩ · cm previously used for saline (Tarotin et al 2019). To account for larger distance from measuring electrode to the fibres, and because the electric field is inversely proportional to the squared distance from the object, the scaling coefficient was chosen to be equal to k el = (r 1 /r 0 ) 2 where r 1 and r 0 are new and previous distances to the fibre respectively equalling to electrode radii. In the single fibre model (Tarotin et al 2019) electrode diameter was 10 μm; the complex nerves had diameters equalling 357 μm for the right rat vagus nerve and a single fibre type nerves (table 2): the number of fibres in them were chosen to match this diameter. The diameter was 700 μm for rat sciatic nerve (table 3) (Schmalbruch 1986). Thus, k el was equal ~1270 and 4900 for these nerves respectively.
The resultant scaling coefficient was (tables 2 and 3): where ρ sal and ρ conn (kΩ · cm) are resistivities of the saline and connective tissue composed of epineurium, endoneurium and perineurium, d el,1 and d nerve (μm) are the diameters of the electrodes in the single fibre simulations and the simulated nerve respectively. For scaling, the final compound dZs found in step 2 of the model (figure 3) were divided by the resultant coefficient k Σ (table 2). Using the sequence of the steps explained above, dZ for all nerves were obtained at 2, 4, 10, 20 and 50 cm from the stimulation site. Experimentally observed noise of 0.5 μV RMS  was added to the resultant compound dZ as the last step ( figure 3). Then, a signal-to-noise ratio (SNR) was computed as a ratio of the amplitude of a pure compound dZ signal to the standard deviation of the added noise. Ten random models with the parameters specified above were computed for statistics. The maximum distance at which dZ can be theor etically measured and possibly imaged was calculated. The impedance change signal was considered measurable if SNR (ratio of mean signal to S.D. of the noise) was more than 1 (⩾0 dB); the dZ was possible to be Table 3. Parameters of the right vagus and sciatic nerves of the rat (Gasser 1950, Rushton 1951, Boyd and Kalu 1979, Soltanpour and Santer 1996, Shimizu et al 2011 imaged if SNR ⩾ 4 (⩾12 dB) ). All the results were compared with the available experimental data using SNR and maximum distances of signals' detection. Source code for statistical models of complex nerves is available online and on the Physiological Measurement website 1 .

Hodgkin-Huxley axons model
The amplitudes of the CAPs simulated with the realistic 3D FEM model (figure 1) decreased with the distance from the AP initiation point (i.p., figure 4): they were ~7 mV at 10 cm, 6 mV at 19 cm, 4.5 mV at 25 cm and 3 mV at 35 cm from the stimulation point. Duration of the negative phase of the compound AP increased from 4 to 5, 6 and 7 ms along the same distances respectively ( figure 4(a)). dZ had the same behaviour. It decreased from −3.5% (−7 μV) at 10 cm, to 2.9% (−5.8 μV) at 19 cm, 2.5% (5 μV) at 25 cm and 1.8% (3.5 μV) at 35 cm from stimulus while duration of its negative phase increased from 9 to 11, 12 and 15 ms at the same distances ( figure 4(a)).
The area under the CAP fell to about 62% of the one at 10 cm from the location of nerve stimulation. In contrast, the area under dZ signal curves was close to constant at these distances ( figure 4(c)).
Necessary statistics (standard deviations) in the calculated areas were obtained with the computation of 100 simplified statistical models with 50 fibres and the same velocity distributions as in the FEM model. The CAP and dZ at DC current in these models were close to the ones obtained with the accurate FEM model (figures 5(b) and 4(a)); they also allowed computation of dZ dispersion recorded at AC (1 kHz in HH case) current. The areas under the CAP computed with these models were (in respect to the ones at 10 cm, figure 4(d)): 83.1% ± 2.9% at 19 cm, 75.6% ± 3.8% at 25 cm and 67.9% ± 4.3% at 35 cm. Impedance changes measured at AC at the same distances were also decreasing with distance: 84.4% ± 2.7%, 77.5% ± 3.3% and 71.5% ± 3.3% with respect to the dZ at 10 cm from the stimulus. The area under the dZ curves measured at DC was constant and independent on the distance from AP initiation point (i.p.) (red line in figure 4(d), left). These values are in close agreement with the results obtained using the accurate FEM model.
However, the dependence of the areas under the compound AP and dZ on the distance from the point of nerve activation were different from the HH nerve (figure 4(c)). In contrast to the HH model, the area under CAP stayed constant; the area under dZ decreased but fluctuated approximately 80% from the maximum value due to the noise associated with weak convergence of the C fibre FEM model requiring to solve 22 PDEs in parallel (see section 4).
As for the nerve model with HH axons, a 50-fibre statistical model with C-nociceptors was implemented to compute statistics for the obtained values. In the developed model, the CAPs and dZ recorded at DC were in a good agreement to the ones simulated with the FEM model (figures 5(b) and 4(b)). This made it possible to implement 100 of these models with the same CV distributions as well as to add impedance change measured with an AC current (2 kHz) which was impossible to do with the computationally-heavy FEM model. The same as in the FEM model, areas under the CAP in respect to the area measured at 0.4 cm were practically constant (blue dashed line at figures 4(d)) at 99.5% ± 0.4%, 99.5% ± 0.5% and 99.3% ± 1.4% at 0.8, 1 and 1.4 cm from the stimulus. In contrast to the HH model, dZ measured at DC decreased with distance from 86% ± 5.3%, to 79.2% ± 5.1% and 68.8% ± 4.0% at the same locations. Compound dZ measured with AC current decreased slower than the one at DC: 97% ± 0.4%, 95.8% ± 0.6%, 93.7% ± 0.8% at the same sites in respect to the dZ measured at 0.4 cm.

Models of mixed diameter fibre nerves
Due to the computational heaviness of the developed 50-fibre FEM models, the simplified statistical models of realistic nerves consisting of nerve fibres of a single type or mixed types were implemented. The parameters utilized for the development of these models were found experimentally in various studies (tables 2 and 3).
Signal-to-noise ratio decreased with distance from stimulation for all nerves, but the maximum distances at which dZ could be measured were larger for nerves with larger and faster fibres. dZ of 0.1 mm nerves consisting of Aα and Aβ fibres were visible at up to at 50 cm from the stimulus (figure 6) with SNR at DC equalling to 2.7 ± 0.2 and 1.8 ± 0.1 respectively at this location (figure 7). Aδ fibres had sufficient SNR = 1.6 ± 0.03 at 20 cm from initiation point but it fell below the noise to 0.6 ± 0.02 at 50 cm. These values were similar for rat sciatic nerve with the SNR falling from 1.6 ± 0.06 at 20 cm to 0.6 ± 0.02 at 50 cm making the signal at this distance undetectable. dZ in the right branch of the vagus nerve of the rat consisting of small Aδ and C fibres was visible at 10 cm (SNR = 1.6 ± 0.05) but was indistinguishable at 20 cm where SNR fell below 1. C fibres dZ was only distinguishable at 1 cm from stimulation point with signal-to-noise equalling to 2.4 ± 0.02 there. SNR obtained with the dZ computed at 2 kHz (figure 5(a)) were in close agreement with the ones at DC (figure 7, table 4).
All these results were in fair agreement with experimental data (see section 3.4). Thus, the theoretically maximum distances of dZ measurement where SNR fell below 1 were (figure 7): >50 cm for Aα and Aβ fibres, ~40 cm for Aδ fibres and the sciatic nerve of the rat, ~15 cm for the rat vagus and about 3 cm for C fibres (table 4). These distances are at the limit of visibility because, at low values of SNR, a lot of averaging will be necessary to distinguish the sought-for signal from the unwanted noise.  (1)).
As was found in , dZ could be reliably imaged when SNR > 4 which significantly decreases the possible distances to 35 cm for Aα and Aβ fibres, ~8 cm for Aδ and rat sciatic, 4 cm for rat vagus and <1 cm for C fibres.

Comparison with experimental data
The maximum distances and SNR simulated for different types of fibres and realistic rat vagus and sciatic nerves closely matched experimental data in the literature. Low level of CAP dispersion in Aα and Aβ fibres where they could be seen at up to a meter from stimulation was shown in Olney et al (1987), Taylor (1993) and Schulte-Mattler et al (2001); it corresponds to the found low dispersion in dZ for these types of fibres.
High dispersion of C fibres was also shown in Freeman (1972), Boone (1995) and Chang et al (2015): in Freeman (1972), CAP of mainly unmyelinated olfactory nerve of the cat could be recorded at up to 2.5 mm from the stimulus; CAP of the vagal C fibres of the mouse was shown to fall >50% at 4 mm (Chang et al 2015). This is similar to the results obtained for 50 C fibres where CAP decreased ~2-fold at 1.4 cm from the site of stimulation Figure 5. (a) Examples of single HH and C fibre action potentials (AP) and impedance changes (dZ) which were used for implementing the simplified 50-fibre models (section 2.1.2) as well as for extension to nerves with realistic morphology (C fibre only, section 2.3). The EIT AC current was 10 times larger than DC to increase the absolute measured dZ. The signals were taken from (Tarotin et al 2019). (b) Compound action potentials and impedance changes at the same distances as were used in the FEM model (table 1). Lines of different colours represent AP and dZ simulated using different statistical models with the same mean and S.D. parameters (tables 2 and 3) to obtain necessary statistics ( figure 4(d)). Velocity distribution of the AP and dZ signals was based on the values of resistivities used in the FEM model (table 1) and their proportionality (1). The signals were summarized and scaled due to the increased diameter of the electrode from a single fibre model to 50-fibre model (table 1, details are in the text).
(figures 4(a) and 5(b)); the 50-fibre dispersion is lower due to smaller number of fibres. CAP dispersion for complex nerves was not simulated in this study, however, its shape for C fibre closely matches the one of dZ at 2 kHz (figure 5(a)) which makes them comparable. Also, in the experiment on unmyelinated crab nerve (Boone 1995), the dZ was shown to be visible up to 16 mm from the onsite. Thus, as the modelled dZ in C-fibre nerve could not be seen further than 1-2 cm (figure 7), these results closely agree with the experiments.
In recently accomplished EIT imaging of rat sciatic nerve , the highest obtained SNR was equal to 8 after averaging which is also in fair agreement with the obtained values at 4 cm for this nerve (6.2 ± 0.1, figure 7) In the same study, the C fibre response was not visible at ~4 cm (length of the rat sciatic nerve) which also agrees with the current study. The literature on the dispersion for right vagus nerve of the rat was not available, however, the results obtained for this nerve can be a reference of what to expect in this nerve as well as in human vagus nerve which is the main aim for neuromodulation of the internal organs it supplies.

Summary of results
(1) Accurate coupled FEM models of the nerves comprising 50 HH axons or 50 C fibres and their simplified statistical equivalents showed different behaviour of compound AP and dZ with distance  (table 2) and with realistic rat sciatic and right vagus nerves (table 3). Compound dZ were formed of single dZ of a C fibre (figure 3) recorded at DC and scaled in time and amplitude depending on the fibre's AP latency and cross-sectional area. The number of fibres in single type models were chosen to fit into a 0.1 mm 2 nerve which approximately equals the size of the rat vagus nerve (table 2). For detailed explanations, see section 2.3. dZ were computed at different distances from stimulation point: 1, 4, 10, 20 and 50 cm (1 cm was omitted for simplicity in A fibres and sciatic nerve); the colour legend is embedded into the graphs. The average experimental white noise with RMS of 0.5 μV  was added to the recordings after the simulations (figure 3); it was omitted in the figure to improve its readability.
from the site of stimulation. In the HH case, the action potentials cancel out with distance due to their biphasicity while the impedance changes measured at DC do not. The inverse was true for the nerve consisting of C fibres due to differences in the shapes of the studied signals. The dZ recorded at AC current showed the decrease with distance in both models (figure 4(d)). Thus, it was revealed that for compound impedance changes to be seen further from the stimulus than compound action potentials, the shape of these signals must satisfy certain condition: dZ needs to be more monophasic than the CAP, as in the case of HH axons at DC (figures 4(a) and 5(a)).
(2) By extension of the study to real nerves consisting of Aα, Aβ, Aδ or C fibres, the theoretical maximal distances at which dZ could be recorded were obtained. These distances were more than 40 cm for a 0.1 mm nerve consisting of A fibres and only up to ~3 cm for the one with C fibres. These findings agree with experimental data on impedance changes and action potentials dispersion (Olney et al 1987, Taylor 1993, Boone 1995, Schulte-Mattler et al 2001. SNR for rat sciatic nerve at the distance of its length (4 cm) obtained in this study was close to the experimental one . In the same study, C fibres could not be measured, which also agrees with the results obtained in the current work (figures 6 and 7). Values for the rat vagus nerve were obtained to serve as an expectation guideline for further studies which need to be carried out with a purpose of its subsequent imaging and selective stimulation.

Limitations and technical difficulties
One of the limitations of the current study was that only one accurate 3D FEM nerve model was undertaken for each nerve type consisting of only 50 fibres. This was due to the lack of computational resources: the C fibre model demanded a system of 20 nonlinear equations (Tigerholm et al 2014, Tarotin et al 2019 to be solved at the 10 6 FEM elements at each time step which demanded around 100 Gb of RAM and a week of computations on a 2-CPU machine. However, the developed statistical simplifications were in a close agreement with the accurate models in terms of CAPs and dZ amplitudes (figures 4(b) and 5(b)) and areas (figures 4(c) and (d)). Therefore, they were used for performing necessary statistical analysis. The extension of the models to different types of nerve fibres and realistic nerves relied on experimentally found distributions of conduction velocities and fibre diameters as well as their relation to each other. These values were based on a limited number of studied nerves and therefore were approximate. Also, the choice of scaling coefficients for the transition from the modelled to experimental conditions was based on limited literature and simple assumptions providing qualitative results which were in a fair agreement with the literature. To obtain such results as well as to predict the behaviour of studied signals on various nerves and experimental conditions was the original purpose of our study.

Answers to the stated questions 4.3.1. Is it possible to record dZ further than CAP?
The level of dispersion of CAP and dZ signals highly depended on the shape of these signals for each particular case. In general, if multiple phases are significantly expressed in the compound action potential of a nerve, the compound dZ may be measured further from the stimulus than the CAP. Such nerves with largely multiphasic extracellularly recorded compound action potentials include the ones containing many fast fibres like rat sciatic nerve ; recordings in unmyelinated nerves , Oh et al 2011 confirmed that their CAPs are monophasic which is in accordance with this study.
Multiphasicity in the action potential does not have the same significance in the impedance change due to the difference in the origin of these signals. Action potentials represent changes in the flow of the current through the membrane while dZ-changes in its impedance. In particular, a negative phase in the AP signal (or positive phase in the extracellular AP, figure 5(a)) appears mainly due to opening of the potassium ion channels; however, this induces decrease in the impedance, the same as when Na channels open during the positive AP phase (negative EAP). A positive phase in the dZ is induced either by deactivation of Na channels in the end of AP in the DC case or phase change in the injected AC current (Tarotin et al 2019); this does not change the shape of the action potential.
These findings are relevant for assessing the feasibility of EIT imaging inside long autonomic nerves for development of neuromodulation techniques (Famm 2013); the vagus nerve is a good target for them as it has access to various internal organs (Berthoud and Neuhuber 2000). For this, SNR of spontaneous activity from internal organs in this nerve should be sufficient at the cervical level.

Are there differences between models in the effect of dispersion on dZ?
In terms of dZ and CAP dispersion, HH giant axons of the squid and mammalian C fibres differed significantly. Because the CAPs of the HH model had larger positive phase than C fibres (figure 5), the APs cancellation in the HH case was quicker than in the C fibre case (figure 4); the converse was true for impedance changes. Consequently, it is possible to measure dZ further than CAP in the HH case but not in the C fibre case. This result can be extended to nerves made of any fibre types: the more multiphasic the CAP is, the more the increment in distance for dZ measurement can be reached in respect to CAP measurement.

What is the largest distance from the site of stimulation at which dZ can be recorded 4.3.3.1.For nerves consisting of Aα, Aβ, Aδ and C fibres
As may be expected from the values of the standard deviations of the conduction velocities of the simulated fibre types as well as their sizes proportional to the amplitude of the dZ (table 2, figure 3), the amplitude of the dZ decreases with distance significantly slower for the nerve consisting of large fast Aα or Aβ fibres. These nerves produce signals significantly larger than the noise level even at half a meter from the stimulation point. dZ in Aδ nerve was found to be detectable at up to 35 cm; the effect of dispersion becomes much stronger for C fibre nerve where the reliable signal can be obtained up to 2-3 cm from the stimulus (figures 6 and 7).
The maximum distance of dZ recording which is where SNR approaches 1 is not equal to the SNR required for dZ imaging. It was shown that for reproducible imaging of fast impedance changes the SNR of 4 is required which substantially decreases the above distances.

For the realistic vagus and sciatic nerves of the rat
The SNR of the dZ in the rat sciatic nerve was shown to decrease to 1 at ~35 cm, in the vagus nerve-at 15 cm; these values defined the maximum distances of dZ measurement for these nerves. The SNR values for the sciatic nerve of the rat show that only A fibres can be imaged at its length (~4 cm) (figure 7) which is in agreement with experiment . The same is valid for the vagus nerve: its SNR is higher than the one for C fibres (~6 at 1 cm and 1 at 15 cm, figure 7) because it also contains A fibres. However, the majority of its fast-myelinated fibres direct into the motor recurrent laryngeal nerve (Gacek et al 1977), they cannot be used as a channel for neuromodulation. Conversely, the autonomic part of the vagus leading to the internal organs is mainly unmyelinated (Agostoni et al 1957), therefore, it is expected to be hardly measured and therefore imaged further than 2-3 cm from the location of stimulation (section 3.3). Thus, an issue arises for imaging its spontaneous activity originating from different organs at a cervical level located tens of cm away. It may be possible to overcome it, by recording close to the stimulus, or concentrating on fast fibres, or perhaps recording for longer periods during changes in state.

Conclusion
Due to variability in the conduction velocities of fibres composing nerves, it is challenging to record compound activity externally at a distance from a point of stimulation. The effect of dispersion is especially strong in unmyelinated fibres whose CAPs cannot be reliably recorded starting from a few centimetres from initiation. The accurate 50-fibre 3D FEM and statistical multi-fibre models developed in this study demonstrated that, for the nerves containing fibres with non-monophasic action potentials, like HH axons, the evoked impedance changes could be measured and possibly imaged with EIT, at greater distances than CAPs could be recorded. The reason is mainly that the bi-phasic action potentials of these fibres cancel out when desynchronised while the impedance changes do not. If taken together with the proportionality of the dZ to nerve cross-sectional area, this enables estimation of the maximal distances at which impedance changes could be measured as well as the SNR expected at these distances. The model predictions are in agreement with the available experimental data.
Systems of differential equations for ionic currents for both types of fibres as well as their gating variables can be found in Hodgkin and Huxley (1952), Tigerholm et al (2014) and Tarotin et al (2019).
The membrane current from the fibres flowing into external space was simulated as I m | Γ = σ e ∇V e · n = C m dV m dt + I ion (V m ), on Γ. (A.5) Here, V m and C m are transmembrane potential and membrane capacitance in (mV) and (μF cm −2 ); I ion (V m ) sum of the ionic currents for HH axon or for C fibre. The current injection via external electrodes was simulated using activating function (Rattay 1989) so that the equation (A.2) becomes where r ax is the radius of the fibre in (cm); ρ i is the resistivity of the axoplasm, (kΩ · cm).