Computational and Mathematical Methods in Medicine
Volume 2020 (2020), Article ID 5436807, 11 pages
https://doi.org/10.1155/2020/5436807

Research Article

A Lead Field Two-Domain Model for Longitudinal Neural Tracts—Analytical Framework and Implications for Signal Bandwidth

G. Fischer,¹ M. Kofler,² M. Handler,¹ and D. Baumgarten1,3

¹Institute of Electrical and Biomedical Engineering, UMIT-Private University for Health Sciences, Medical Informatics and Technology, Eduard-Wallnoefer-Zentrum 1, 6060 Hall in Tirol, Austria
²Department of Neurology, Hochzirl Hospital, Zirl, Austria
³Institute of Biomedical Engineering and Informatics, Technische Universitaet Ilmenau, Ilmenau, Germany

Correspondence should be addressed to G. Fischer; gerald.fischer@umit.at

Received 30 November 2019; Revised 7 April 2020; Accepted 21 April 2020; Published 29 May 2020

Academic Editor: Valeri Makarov

Copyright © 2020 G. Fischer et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Somatosensory evoked potentials are a well-established tool for assessing volley conduction in afferent neural pathways. However, from a clinical perspective, recording of spinal signals is still a demanding task due to the low amplitudes compared to relevant noise sources. Computer modeling is a powerful tool for gaining insight into signal genesis and, thus, for promoting future innovations in signal extraction. However, due to the complex structure of neural pathways, modeling is computationally demanding. We present a theoretical framework which allows computing the electric potential generated by a single axon in a body surface lead by the convolution of the neural lead field function with a propagating action potential term. The signal generated by a large cohort of axons was obtained by convoluting a single axonal signal with the statistical distribution of temporal dispersion of individual axonal signals. For establishing the framework, analysis was based on an analytical model. Our approach was further adopted for a numerical computation of body surface neuropotentials employing the lead field theory. Double convolution allowed straightforward analysis in the frequency domain. The highest frequency components occurred at the cellular membrane. A bandpass type spectral shape and a peak frequency of 1800 Hz was observed. The volume conductor transmitting the signal to the recording lead acted as an additional bandpass reducing the axonal peak frequency from 200 Hz to 500 Hz. The superposition of temporally dispersed axonal signals acted as an additional low-pass filter further reducing the compound action potential peak frequency from 90 Hz to 170 Hz. Our results suggest that the bandwidth of spinal evoked potentials might be narrower than the bandwidth requested by current clinical guidelines. The present findings will allow the optimization of noise suppression. Furthermore, our theoretical framework allows the adaptation in numerical methods and application in anatomically realistic geometries in future studies.

1. Introduction

Somatosensory evoked potentials (SEPs) are a well-established tool for studying the conduction of stimulated peripheral activity in sensory pathways to the somatosensory cortex in clinical routine [1–3]. They allow the evaluation of volley conduction along the entire sensory pathway. For example, signals evoked by a stimulus applied to the tibial nerve at the ankle may be obtained at the popliteal fossa, over the sacral, lumbar, thoracic, and cervical spinal cord and over cortical regions. An early standard for performing such recordings is defined in [4] and provides the basis for routine recording up to now.

While in [4], bandwidth recommendations go up to 3 kHz, more recent experimental studies on high-frequency components of neural signals focus on a more narrow range of only several hundreds of hertz [5–7]. Most of these studies investigated intracortical data. As high-frequency components appear damped in a lead which is distant from the source, only a few studies exist reliably detecting components in a high-frequency band on the body surface [8].

Besides experimental data, modeling can be used for predicting spectral properties of bioelectric signals. In particular, in cardiac electrophysiology, model-based spectral analysis is a well-established tool [9, 10], and recent studies show impressive correlation of simulated and experimental data even for relatively small atrial signal components [11].

Clark and Plonsey [12] developed an early analytical model for computing the extracellular potential of a single axon. Their model was based on a partial differential equation (PDE). However, due to the complex microanatomical structure of nerves or neural tracts consisting of thousands of axons organized in fascicles, PDE-based models are computationally demanding. They require modeling of each single fiber. Some recent studies have applied PDE models for studying stimulation of peripheral nerves by electrodes [13, 14]. However, for studying stimulation, only relatively short time intervals (significantly less than 1) need to be considered. For modeling lead field evoked neural activity (i.e., activity recorded by a lead located distantly from the sources of the electric field), much longer time intervals (significantly more than 10) are required rendering PDE approaches computationally unattractive.

In this study, we present a novel theoretical framework using distributed dipoles along a central axis of a nerve, fascile, or neural tract for efficiently computing the lead field generated by such a structure. Thus, a tool is provided, which allows for multiscale modeling by linking cell models with body surface potentials. Our approach allows using statistical distribution functions for computing the potential field generated by a large number of axons. Furthermore, the use of Fourier transform (FT) allows a straightforward investigation of spectral signal properties.

We aimed at developing our theoretical framework for analytic expressions and, therefore, chose a simplified half-space geometry which approximates longitudinal conduction of a volley within a spinal tract. An early study revealed low-amplitude () traveling waves generated by longitudinal spinal volley conduction, which can be observed essentially along the thoracic compartment of the spine [15]. It must not be confused with stationary synaptic spinal signals in the lumbosacral region (generated within the spinal cord gray matter) and in the cervical region (generated by subcortical structures such as the brain stem). We will discuss the generalization of our approach to other longitudinal neural structures at the end of the manuscript.

2. Methods

In this section, we present the key features of our analytical framework. All detailed mathematical derivations are listed in the Appendices (available here) which are provided as an online supplement.

2.1. Lead Field Two-Domain Model—Concept

The concept underlying the lead field two-domain model is depicted in Figure 1. A neural tract contains several hundreds or thousands of axons. For a volley, action potentials traveling along the tract activate individual axons with a significant temporal dispersion. Here, significant dispersion means that the variation in the onset of activation may be larger than the duration of a membrane action potential pulse (de- and repolarization) and individual axons may fire even more than one pulse.

Figure 1: Concept of the lead field two-domain model. (a) A neural tract of diameter contains a large number of individual fibers (light red) surrounded by myelin layers (gray). (b) For computation of a lead field potential in the distant field point , the individual fibers can be considered on a central axis of the neural tract (red, active domain). The current return in the volume conductor of average bulk conductivity generates the extracellular potential field .

We aimed at computing the electric potential in the lead field of a neural tract. Here, the distance of a field point (i.e., a lead) from the tract is significantly larger than a typical diameter of the tract. At a sufficiently large distance , the individual location of each axon inside the tract becomes negligible. For computing a lead field potential, we assumed that all axons are located on a central axis of the tract. A membrane action potential traveling along individual axons drives intracellular currents. Typically, action potentials are conducted with a speed in the order of 40 ms⁻¹ to 80 ms⁻¹ and have a duration of 0.5 ms to 1.0 ms. Thus, the activated segment of each axon has a length of some centimeters. Note that this is much larger than the diameter of an individual axon, larger than the distance between adjacent nodes of Ranvier, and larger than the diameter of the entire tract. Therefore, we neglected anatomical microstructures such as myelin layers and assumed an average bulk conductivity for the extracellular current return pathway.

Thus, we considered two domains: an active intracellular domain and a passive extracellular or bulk domain. In the first step, we computed the signal generated by a single axon. In the second step, we applied superposition of dispersed axonal signals for obtaining the compound action potential (CAP) of a neural tract.

2.2. Axonal Dipoles Model

We considered the signal generated by a single axon located at the central axis of a neural tract. Figure 2 gives an overview on the axonal dipole model. Geometrically, the axon is located within a conducting half-space parallel to the body surface at the depth . Along the axon ( direction), the membrane action potential propagates with a velocity . This generates spatial gradients driving impressed intracellular currents. We interpreted these source currents as a distribution of moving dipoles creating an electric potential in the surrounding tissue and at the body surface. The body surface potential created by a unit dipole at the observation electrode (located at the origin) was obtained by

Figure 2: (a) Geometry of the model. The axon is located at a distance from the body surface. The observation electrode is placed at the origin. (b) Membrane potential along the axon. At the peak of is located at . The membrane potential propagates with a velocity in the direction. Membrane potential gradients generate distributed moving dipole sources . They are shown in (b) and (a). Due to the opposite orientation of de- and repolarization, a quadrupole source pattern is obtained. (c) The lead field is an odd function with a single oscillation. The distance between the two extrema is .

Dipole strength is proportional to the first temporal derivative of the membrane potential (see Appendix A.1). We translated the movement of the membrane action potential along the axon into a convolution term and obtained the body surface potential created by the axon

Here, is the conductivity of the intracellular space inside the axon and is the axon diameter. For integrating all dipoles along the axon dimension , we introduced the auxiliary time shift and replaced by .

Applying convolution theory and time continuous Fourier transform (FT), the frequency spectrum of the axonal signal was obtained by the product of the FT of two terms. The first (source) term is the temporal derivative of the membrane potential. The second (volume conductor) term reflects the lead field of a moving dipole source

2.3. Neural Pathway CAP

A neural pathway is formed by a large number of axons. For a volley traveling along this pathway, activation displays a temporal dispersion in individual axons. This relates to a time shift of action potentials between axons. For enabling an analytical treatment, we first assumed equal diameter and equal conduction velocity in all axons. The remaining variation in axonal activity is modeled by a statistical distribution of activation dispersion .

The spectrum of the CAP was obtained by multiplying the spectrum of a single axon with the spectrum of the statistical distribution of the dispersion (see Appendix B)

In the time domain, the CAP was obtained by a convolution of the signal of a single axon with the distribution of time shift . Assuming a normal distribution of time shifts (standard deviation ), we obtained in the frequency domain (see Appendix B)

The spectrum of the distribution function (5) has a maximum at (DC) and continuously converges to zero with increasing frequency, acting as a low-pass component. This holds true also for any other statistical distribution of the dispersion. According to basic theory, the Fourier transform of a probability density function is called the characteristic function of the distribution. In Appendix B.1, we summarize some basic theorems on characteristic functions [16], proving the low-pass characteristic for any distribution. Appendix B.2 depicts examples of characteristic distribution functions in order to exemplify the described low-pass filter effect.

In the time domain, the CAP signal was obtained by inverse Fourier transform of (4).

2.4. Model Analysis

We observed that our analytical model represents the frequency spectrum of the lead field potential by the product of three different terms:

(i)The membrane potential derivative (source term, obtained from a cell model)(ii)The lead field (moving dipoles)(iii)The dispersion term

Figure 3 displays this observation schematically. We analyzed the properties of each of the three components and their interaction.

Figure 3: Signal generation model in the frequency domain (see text).

Source term. Figure 4 exemplarily depicts the action potential of a sensory axon [14] together with its first temporal derivative and the frequency spectrum of the derivative. The direct current content in the spectrum is zero as the membrane potential always starts at resting potential before activation, and it returns again to resting potential after activation.

(a)

(b)

Figure 4: (a) Membrane action potential (dashed) of a sensory axon [14] and its temporal derivative (solid). For each trace an individual scale was applied on the ordinate. Time basis was chosen such that the potential maximum is at . (b) Normalized absolute amplitude of the frequency spectrum of the action potential derivative in a double logarithmic plot. The maximum was observed at 1.83 kHz. The amplitude converged to zero for low and high frequencies.

Furthermore, the derivative of an action potential traveling along the axon is a piecewise continuous function. Thus, we concluded from basic theory of Fourier transform that the spectrum approaches zero also with increasing frequency (see Figure 4 and Appendix B). Thus, the frequency spectrum of the source term is of “bandpass type.” We obtained a coarse estimate of the peak frequency in the spectrum of the derivative by taking the inverse of the time interval of membrane activation (i.e., the duration of an activation cycle from the beginning of depolarization to the end of repolarization)

Lead field. The lead field transmits the source term to the potential at the observation electrode . The lead field is an odd function (see Figure 2). Thus, the integral under the function yields zero, and again, there is no DC content in the spectrum. Furthermore, the volume conductor term is continuous. Again, the amplitude spectrum approaches zero with increasing frequency. Thus, also the frequency spectrum of the lead field term is of “bandpass type.”

The lead field term yields a function with a single oscillation and point-symmetric shape (see Figure 2). For coarsely estimating a center frequency of this peak, we considered the distance between the two extrema of the lead field (). Within this interval, half oscillation occurs. Thus, the center frequency was estimated by taking the inverse of twice the time a dipole needs for passing the distance between the two extrema

Here, was a property of a model component, not of the signal. For coarsely estimating a center frequency of the axonal signal (i.e., the signal at the output of the moving dipole component in Figure 3), we considered also the duration of an activation cycle :

Dispersion. For obtaining a coarse estimate of the center frequency of a lead field compound action potential of a neural pathway, we considered that temporal dispersion increases the duration of the CAP signal with respect to the duration of an axonal signal. The width of the dispersion pulse was estimated by as about 95% of the individual activation pulses are within the interval of ±2 from the center of the normal distribution [17].

Only four biophysically relevant parameters enter this estimation: duration of an activation cycle , conduction velocity , depth , and standard deviation of dispersion .

2.5. Simulations

For numerical evaluation of axonal potentials (2), we replaced integration by discrete summation. In Appendix A.2, we describe this in detail and generalize our concept to multiple observational electrodes or leads located at the coordinates on the body surface. The discrete model is depicted in Figure 5. As outlined in Appendix A.2 a linear equation was established for computing the discrete axonal potentials

Figure 5: Schematic drawing of the discrete model underlying the numerical approximation. At the body surface, the potential is calculated in electrodes (black circles, coordinates ). The axon is subdivided into source points at spacing (red circles).

Here is the lead field matrix relating the source points to the leads and is the dipole matrix holding the primary source terms. The propagation of the membrane action potential at constant speed (schematically depicted in Figure 2) causes a systematic spatial displacement of the source dipoles in each row (i.e., at each time step) of . This is equivalent to the convolution in the temporal domain.

We simulated the action potential of a sensory axon applying the ionic current model described in [14] by using the software NEURON (http://www.neuron.yale.edu/neuron/). For minimizing numerical effects in our analysis, we selected fine grids for discretization. According to the analysis in Section 2.4, we expected highest frequencies at the cell membrane. Therefore, we chose the smallest time step 2 µs for investigating the properties of the membrane action potential. For simulating body surface potentials, we increased the time step to 5 µs. This yielded a Nyquist-frequency of 100 kHZ which is sufficient for accurately investigating spectral features of the body surface potential up to several kHz. We linked spatial discretization to the chosen time step by selecting . We simulated all signals in a time interval of -100 ms to 100 ms. We applied time discrete FT for computing all spectra obtaining a frequency resolution of 5 Hz.

In the first step, we varied the axon diameter in three steps (7 µm, 10 µm, and 13 µm) in the NEURON model for investigating the effect on the membrane action potential and conduction velocity. For studying the dependency of results from physiologic variation, we assigned three different values for the depth of the neural pathway (35 mm, 50 mm, and 65 mm) [18, 19] and for the conduction velocity (45 ms⁻¹, 60 ms⁻¹, and 75 ms⁻¹) [14, 15].

According to equation (2), the amplitude of the axonal potential depends on the ratio of the intracellular conductivity to the bulk conductivity . The intracellular conductivity is the inverse of the cytoplasmic resistivity. We used the same value as it was used for the NEURON model of sensory axons [14, 20] at body temperature. Bulk conductivity is significantly smaller and varies considerably between tissue types [14, 21]. Thus, we assumed a mean ratio . In order to account for the uncertainty in this ratio, we doubled and halved this value in our computations.

For all combinations of parameters, we computed the axonal potential . For illustrating the spatial field distribution, we simulated the potential for 10 µm axon diameter at 50 mm depth within a square of 200 mm edge length centered at the origin (2.5 mm grid spacing).

For simulating neural potentials, we assigned three values to the standard deviation of temporal dispersion (1.0 ms, 1.25 ms, and 1.5 ms) and assumed the largest dispersion at the slowest conduction velocity. We chose an axon diameter of 10 µm [14, 22] and assumed that axons define a neural pathway.

3. Results

Source term. For the three axon diameters (7 µm, 10 µm, and 13 µm), the NEURON model delivered a conduction velocity of 43.6 ms⁻¹, 60.0 ms⁻¹, and 75.1 ms⁻¹, respectively. Variation in diameter had only a minor effect on the action potential. For all three simulations, peak-to-peak amplitude (maximum at activation vs. minimum during hyperpolarization) was in the range of 107.9 mV to 108.8 mV. The time span between the two potential extrema was in the range of 480 µs to 482 µs. Time for activation (i.e., from crossing of the -70 mV level to the maximum) was in the range of 92 µs to 94 µs. Thus, we obtained a pulse duration (de- and repolarization) of yielding  Hz. We arbitrarily assigned to the maximum of the membrane action potential. As predicted by our model analysis, the highest frequency in the signal generation chain was observed at the cell membrane (derivative of the action potential). Here, we obtained a bandpass type frequency spectrum with significant signal components in the kHz range (see Table 1 and Figure 4).

Table 1: Results for parameter setting  mm ms⁻¹, and  ms.

Peak frequency
Bandwidth
Frequency estimates
Amplitude

(Hz)
(Hz)
(Hz)
(peak-to-peak)

AP derivative,
1830
547 to 4545
1750
3.1 

Axon potential, φA
300
100 to 661
338
1.24 

CAP, φN
125
50 to 227
126
0.32 

Lead field. For the three axon diameters (7 µm, 10 µm, and 13 µm), we obtained peak-to-peak amplitudes (global maximum vs. minimum) of 0.44 nV, 1.24 nV, and 2.61 nV, respectively, at and a depth of 50 mm. Taking into account also the uncertainty in conductivity and variations in depth, body surface peak-to-peak amplitude generated by single spinal cord fibers may vary in the range of approximately 0.1 nV to 10 nV. The peak frequencies in the frequency spectrum of the three axonal potentials were 215 Hz, 300 Hz, and 315 Hz, respectively. Peak frequency increased with increasing diameter due to the increase in conduction velocity.

The potentials computed for an axon diameter of 10 µm at a rectangular grid are depicted in Figure 6. The field pattern reflects the underlying quadrupole source. When displacing an observation point along the axon, the signal displayed a temporal offset, while its morphology was preserved. In the rectangular direction, amplitude decreased with increasing distance from the source, as did the peak frequency. In the spatial and temporal domain, activation of the axon underneath an observation point was reflected by a potential minimum. Due to the asymmetric shape of the action potential, the minimum was slightly shifted away from . However, this effect amounted to only 0.06 ms. Two small positive peaks were observed, before and after the activation minimum, respectively. This signal morphology was caused by the quadrupole source pattern shown in Figure 2. The morphology was slightly asymmetric with the second maximum being slightly larger than the first. The difference in the amplitude of the positive peaks was mainly caused by the small hyperpolarization wave (about 2 mV).

(a)

(b)

Figure 6: (a) Axonal potential generated by a single axon of 10 µm diameter at a depth of 50. The conductivity ratio was set to . Isopotential lines are shown in steps of 0.1 nV. Circles indicate the location of the four observation leads which are depicted in (b). (b) The green solid trace depicts the axonal potential at the origin. Shifting the observation lead by 50 mm along the axon (green dashed line) introduced a time shift of 0.83 ms while preserving signal morphology. Lead displacements in a direction perpendicular to the axon are indicated by a magenta trace () and an orange trace (). With increasing distance from the axon amplitude decreased (peak-to-peak amplitude 1.24 nV, 0.89 nV, and 0.44 nV). Furthermore, the shape of the negative peak broadened which was also reflected by a decreasing peak frequency (300 Hz, 265 Hz, and 210 Hz).

Figure 7 depicts the axonal potential at a fixed conduction velocity of 60 ms⁻¹ for variable depths. As expected, axonal body surface amplitude was inversely proportional to approximately the square of the depth.

(a)

(b)

Figure 7: (a) Axonal signals for a conduction velocity of 60 ms⁻¹. Color coding blue, green, and red correspond to a depth of 35 mm, 50 mm, and 65 mm, respectively. (b) Absolute amplitude obtained for the frequency spectrum of the three axonal signals. Color coding is identical as in the time domain. The spectral amplitude was normalized by the peak value obtained at .

In our simulations, the volume conductor (lead field term) reduced the peak frequency of axons signals on the body surface. For the chosen parameter range, axonal peak frequency was between 200 Hz and 500 Hz. All axonal signals displayed spectra with a single peak as predicted by our model analysis. For this reason, both at low near DC frequencies and at high frequencies above several kHz, the amplitudes in the spectrum were very low. With increasing depth, the peak frequency decreased, in accordance with the prediction obtained from (8) for the estimated center frequency .

Dispersion. The low-pass filter effect of temporal dispersion described in Section 2.3 further reduced the frequency content. For the CAP, peak frequency was in the range of 90 Hz to 170 Hz. The estimated center frequency was slightly above peak frequency and well within the computed frequency band. Thus, the frequency band was essentially defined by the inverse of the time span in which the majority of activation pulses pass the segment of length. The peak-to-peak amplitude of the signal was in the order of magnitude of , which is a typical value for short latency somatosensory evoked potentials. Despite the assumption of 2000 axons defining a neural pathway, the CAP increased only by a factor of 200 compared to a single axon. This loss of amplitude was due to temporal dispersion, as individual axons generated their signals in an asynchronous fashion.

Figure 8 depicts the potential of a neural pathway as a superposition of many axonal signals as defined in (4). Here, in accordance with clinical practice, negativity is plotted upwards. The width of the signal increased in relation to temporal dispersion in the volley. In the frequency spectrum, increased width was reflected by a shift of the maximum towards a lower frequency. The morphology of the signal again reflected the underlying quadrupole source pattern. Again, the amplitude was inversely related to the depth. However, the influence of depth was less than a power of two. This can be explained by the large spatial distribution of a volley along the axon.

(a)

(b)

Figure 8: (a) Signals of the neural pathway for a conduction velocity of 60 ms⁻¹. Negative amplitude was plotted upwards. Color coding is the same as in Figure 7. The time interval between the two positive peaks was in the range of 3 ms to 4 ms for all simulations. (b) Absolute amplitude obtained for the frequency spectrum of the three axonal signals. Amplitude was normalized by the peak value obtained at .

Table 1 depicts results obtained with the middle range of assumed parameters (i.e., ). FT delivered peak frequencies and bandwidth (% of peak amplitude), which are listed together with the estimates (as defined in Section 2.4). The signals and their spectra are shown by the green traces in Figures 7 and 8.

Table 2 depicts the parameters for the highest, lowest, and middle peak frequencies in our simulation. For investigating the correlation of the peak frequency with the estimated center frequency , we computed the coefficient of determination from all parameter combinations. For the axonal signals, we obtained , and for the CAPs, we obtained . Thus, the estimated center frequencies as defined in Section 2.4 accurately reflected the effect of the parameters on the signal spectrum.

Table 2: Simulation parameters and characteristic frequencies.

Depth
Conduction velocity
Peak frequency
Bandwidth
Frequency estimates

(mm)
(ms⁻¹)
(Hz)
(Hz)
(Hz)

35
75
170
69 to 302
169

50
60
125
50 to 227
126

65
45
90
36 to 173
94

4. Discussion

With the increasing distance from a source to an observation electrode, spatial smoothing of the lead field reduced the high-frequency content in the axonal signal. The frequency band was in the range of several hundred Hz. This reduction was caused by the Laplace term which underlies each biopotential computation [23]. Note that this was a purely geometrical effect which was also described for body surface cardiac signals [10]. From our model analysis, we interpreted the obtained frequency band as follows: the observation electrode was most sensitive for signals in an axonal segment of length underneath the electrode. The inverse of the time, which an activation pulse needed in order to pass this segment, was an estimator for the center of the frequency band. The peak-to-peak amplitude of a signal generated by a single axon on the body surface was in the order of nanovolts. Thus, our model predicted that individual activity of single axons cannot be observed from the body surface.

Any model is based on idealized assumptions and this is also true for an analytic model. In order to allow for an analytic treatment, we restricted our analysis to axons of equal diameter and a simplified half-space geometry. Thus, we limited our analysis in the first step to the traveling wave along the spinal cord—a waveform which is scientifically well described [15] but of restricted clinical use due to its low amplitude. The signals obtained from our model were in reasonable agreement with experimental data:

(i)CAP amplitudes were in the order of some tens of microvolts [15]. Here, the values chosen for conductivity and for the number of axons entered the result. There is considerable uncertainty in these parameters, and thus, only coarse estimates were obtained for the amplitude. The intracellular conductivity we used in our simulation was the inverse of the cytoplasmic resistivity as used in the cell model described in [14]. Compared to cell models, direct measurement of cytoplasmic resistivity yields similar values [24] when extrapolating data to body temperature [25]. The chosen bulk conductivity was an average value of different types of tissues which surround axons or nerves. Here, the selected conductivity ratio corresponds to a bulk conductivity of approximately 0.2 Sm⁻¹ which is a frequently used value for average tissue conductivity [21]. We assumed that 2000 axons contributed to the signal generated by a spinal tract. As can be seen from Appendix C, we estimated that a neural structure containing 2000 axons had a diameter of 0.8 mm which appeared being a reasonable dimension. The total macaque pyramidal tract contains approximately 600000 axons [22]. We assumed a comparable number of axons in the sensory tracts, of which only a small portion is activated by stimulation of a peripheral nerve. Thus, our model allows for a coarse estimate of amplitude at reasonable tract dimensions. However, variation of parameters showed that there is still a huge uncertainty in amplitude estimations. Future work is needed for reducing this variability(ii)For single axons, the “quadrupole” signal morphology qualitatively agreed with early theoretical predictions [12] and recent experimental findings [26]. For the CAP, the dominating negative peak of activation was also in agreement with experimental findings [15].(iii)Our analytical framework predicted that the width of the temporal dispersion function is the main factor determining the width of the evoked CAP on the body surface. By selecting in the range of 4 ms to 6 ms, the width of the computed CAP was in the order of 5 ms to 10 ms which is comparable to experimental observations [15]. Here, peak or central frequencies were typically in the order of 100 Hz to 200 Hz. Literature on experimentally observed dispersion functions is largely lacking. Basic theory on characteristic functions of probability density functions predicts that the low-pass filter effect widening the signal trace will be observed for any dispersion function. Due to the limitations of the analytic model, we restricted our analysis to the normal distribution as a basic model of statistical distribution. As discussed below, numerical methods should be applied in future studies for addressing the temporal dispersion in combination with variable axon diameter and conduction velocity. In the first step, the chosen parameters allowed a reasonable prediction of the width of the CAPs and their central frequency band.

The simplifications described above did not allow using our model for studying signals generated by neurophysiologically relevant structures such as the spinal gray matter or subcortical signal generators [2, 15]. However, as there was some similarity between the signal morphology and width of lumbosacral stationary potentials with the CAP predicted by our models, our results might be also useful for a first estimation of spectral properties of lumbosacral postsynaptic EPs.

Our analysis on the 50% bandwidth of the signal suggested that most of the signal was contained in a relatively narrow band of 40 Hz to 300 Hz. Bandwidth was significantly more narrow than suggested in the guidelines (20 Hz to 3000 Hz [4]) and in agreement with current scientific work, which focuses on the range of several hundred hertz for neural high-frequency signals [5–8] rather than the kilohertz range. In a routine setting, a more strictly selected bandwidth might allow reducing the influence of technical noise sources such as thermal noise (electrode impedance and amplifier). Here, future experimental work is needed. With respect to physiological background activity, the CAP spectrum displayed a strong overlap with the activity of skeletal muscles and a weaker overlap with cardiac activity [27]. Signals generated by skeletal muscles are within a broad frequency band with significant spectral components in the range of 100 Hz to 400 Hz [28, 29]. Cardiac body surface potentials contain frequencies up to 150 Hz. However, the most significant ECG frequency components are at relatively low frequencies. For a standard 12-lead ECG, the bandwidth is normally reduced to 100 Hz, and for ECG monitoring in critical care units, the bandwidth is even reduced to 40 Hz [30].

5. Conclusions and Outlook

We developed a theoretical framework which allows the computation of lead field potentials generated by a neural tract, nerve, or fascicle beyond considering each single axon or fiber. Our approach allows linking membrane potentials obtained from ionic current models with body surface potentials for multiscale modeling. The neural structure was modeled by its central axis, and individual membrane action potential pulses were considered by their statistical distribution in time. Our results suggest that the signals generated by individual axons are sufficiently small and that the number of axons is sufficiently large for obtaining a reasonable approximation by statistical analysis.

Similarly, as for cardiac tissue [21, 31] or stimulation of cortical tissue [32], our model considers two domains: an active intracellular domain (membrane potential , intracellular conductivity ) and an extracellular bulk domain (potential , conductivity ). In contrast to existing models, our approach focuses on lead field activity modeling the source field by dipole distributions. Then, the passive potential is obtained by integration. Here, convolution allows solving equations in the time and frequency domains.

The treatment in the frequency domain allowed gaining insight into the spectral properties of the generated signals. Highest frequencies occur directly at the cell membrane with peak frequencies in the order of kilohertz. The volume conductor transmitting the signal to a lead field electrode acts as a bandpass reducing frequency content with increasing distance between the source and sensor. Furthermore, temporal dispersion of individual activation pulses acts as an additional low-pass filter reducing peak or central frequency in the signal to approximately 100 Hz to 200 Hz. This suggests that the clinical bandwidth of 20 Hz to 3000 Hz [4] may be unnecessarily broad. Future experimental and modeling work is needed for potentially optimizing recommendations for bandwidth and filter types. This may significantly improve the clinical practicability of SEP recordings by reducing the number of samples needed for obtaining reliable recordings.

For establishing the concept of our theory, we applied analytic treatment. This limits analysis to a model of a simple half-space geometry. The results apply to the low-amplitude propagating wave in spinal evoked potentials.

The insight obtained from developing an analytical framework for a lead field two-domain model of a neural tract or nerve provides a sound basis for using numerical tools in the near future for generalizing the approach. This will allow the investigation of a broader spectrum of anatomically and electrophysiologically realistic modeling applications:

(i)Numerical field computation schemes such as the finite element method (FEM) and the boundary element method (BEM) [33–35] can be used for computing the lead field matrix generated by a unit dipole traveling along a central axis of a neural pathway. Such methods allow considering realistic anatomy and variation of conductivity in tissue compartments. By convolution of such a numerically computed lead field matrix with the primary source matrix the axonal potential can be computed in general geometries. For example, such methods can be used for computing the potential generated by an axon of a peripheral nerve in a limb.(ii)Our approach allows numerical computation of axonal signals for fibers of varying diameter (and thus, varying conduction velocity). The model can be generalized to neural tracts containing fibers of varying diameter, by applying numerical tools for summation of many axonal signals (instead of statistical methods as used here). This will allow predicting dispersion at the stimulation site of a peripheral nerve applying computer modeling as described in [13, 14]. By considering individual conduction properties of fibers along the nerve, the lead field two-domain model can be generalized for computing CAPs of peripheral nerves or spinal tracts.(iii)Dispersion functions may be a powerful tool for computing also the lead field potentials of stationary signal generators by direct convolution of a membrane potential derivative with a convolution function.

Thus, the lead field two-domain model has the potential for providing an important novel model-based research tool for evoked or stimulated neurophysiological signals.

Data Availability

The data contained in the article were obtained from simulations which were encoded in Matlab based on the analytical derivations described in the article. All data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Supplementary Materials

The supplementary file is named R1_Appendix_CompMathMethMed. (Supplementary Materials)

References

M. Hardmeier, L. Leocani, and P. Fuhr, “A new role for evoked potentials in MS? Repurposing evoked potentials as biomarkers for clinical trials in MS,” Multiple Sclerosis, vol. 23, no. 10, pp. 1309–1319, 2017.

M. Kofler, “Evozierte Potentiale,” in Neurologie Compact, A. Hufschmidt, C. Luecking, S. Rauer, and F. Glocker, Eds., pp. 777–787, Thieme Verlag, Stuttgart, 2017.

T. Sand, M. B. Kvaløy, T. Wader, and H. Hovdal, “Evoked potential tests in clinical diagnosis,” Tidsskrift for den Norske Lægeforening, vol. 133, no. 9, pp. 960–965, 2013.

M. R. Nuwer, M. Aminoff, J. Desmedt et al., “IFCN recommended standards for short latency somatosensory evoked potentials. Report of an IFCN committee,” Electroencephalography and Clinical Neurophysiology, vol. 91, no. 1, pp. 6–11, 1994.

T. Fedele, C. Schonenberger, G. Curio, C. Serra, N. Krayenbuhl, and J. Sarnthein, “Intraoperative subdural low-noise EEG recording of the high frequency oscillation in the somatosensory evoked potential,” Clinical Neurophysiology, vol. 128, no. 10, pp. 1851–1857, 2017.

B. E. Mouthaan, M. A. van ‘t Klooster, D. Keizer et al., “Single pulse electrical stimulation to identify epileptogenic cortex: clinical information obtained from early evoked responses,” Clinical Neurophysiology, vol. 127, no. 2, pp. 1088–1098, 2016.

B. Telenczuk, R. Kempter, G. Curio, and A. Destexhe, “Refractoriness accounts for variable spike burst responses in somatosensory cortex,” eNeuro, vol. 4, no. 4, 2017.

T. Götz, T. Milde, G. Curio et al., “Primary somatosensory contextual modulation is encoded by oscillation frequency change,” Clinical Neurophysiology, vol. 126, no. 9, pp. 1769–1779, 2015.

G. Fischer, L. Wieser, and F. Hintringer, Modeling mother rotor anchoring in branching atrial tissue, J.-I. Choi, Ed., Atrial Fibrillation, Intech Open, Rijeka, 2012.

C. N. Nowak, G. Fischer, L. Wieser, B. Tilg, A. Neurauter, and H. U. Strohmenger, “Spatial-temporal filter effect in a computer model study of ventricular fibrillation,” Biomed Tech, vol. 53, no. 4, pp. 163–173, 2008.

M. Potse, T. A. Lankveld, S. Zeemering et al., “P-wave complexity in normal subjects and computer models,” Journal of Electrocardiology, vol. 49, no. 4, pp. 545–553, 2016.

J. Clark and R. Plonsey, “The extracellular potential field of the single active nerve fiber in a volume conductor,” Biophysical Journal, vol. 8, no. 7, pp. 842–864, 1968.

P. Schier, M. Handler, L. Johnson Chacko et al., “Model-based vestibular afferent stimulation: evaluating selective electrode locations and stimulation waveform shapes,” Frontiers in Neuroscience, vol. 12, p. 588, 2018.

J. L. Gaines, K. E. Finn, J. P. Slopsema, L. A. Heyboer, and K. H. Polasek, “A model of motor and sensory axon activation in the median nerve using surface electrical stimulation,” Journal of Computational Neuroscience, vol. 45, no. 1, pp. 29–43, 2018.

S. J. Jones and D. G. Small, “Spinal and sub-cortical evoked potentials following stimulation of the posterior tibial nerve in man,” Electroencephalography and Clinical Neurophysiology, vol. 44, no. 3, pp. 299–306, 1978.

E. Jondeau, “Characteristic functions and fourier transforms,” in Financial Modeling under Non-Gaussian Distributions, E. Jondeau, S. Poon, and M. Rockinger, Eds., pp. 477–486, Springer London, London, 2007.

E. Kreszig, Advanced Engineering Mathematics, John Wiley and Sons, eighth edition, 1999.

I. Busscher, J. J. Ploegmakers, G. J. Verkerke, and A. G. Veldhuizen, “Comparative anatomical dimensions of the complete human and porcine spine,” European Spine Journal, vol. 19, no. 7, pp. 1104–1114, 2010.

S. H. Zhou, I. D. McCarthy, A. H. McGregor, R. R. H. Coombs, and S. P. F. Hughes, “Geometrical dimensions of the lower lumbar vertebrae–analysis of data from digitised CT images,” European Spine Journal, vol. 9, no. 3, pp. 242–248, 2000.

C. C. McIntyre, A. G. Richardson, and W. M. Grill, “Modeling the excitability of mammalian nerve fibers: influence of afterpotentials on the recovery cycle,” Journal of Neurophysiology, vol. 87, no. 2, pp. 995–1006, 2002.

G. Fischer, B. Tilg, R. Modre et al., “A bidomain model based BEM-FEM coupling formulation for anisotropic cardiac tissue,” Annals of Biomedical Engineering, vol. 28, no. 10, pp. 1229–1243, 2000.

L. Firmin, P. Field, M. A. Maier et al., “Axon diameters and conduction velocities in the macaque pyramidal tract,” Journal of Neurophysiology, vol. 112, no. 6, pp. 1229–1240, 2014.

P. Malmivuo, J. Malmivuo, and R. Plonsey, Bioelectromagnetism — Principles and Applications of Bioelectric and Biomagnetic Fields, Oxford University Press, 1995.

C. H. Fry, S. C. Salvage, A. Manazza et al., “Cytoplasm resistivity of mammalian atrial myocardium determined by dielectrophoresis and impedance methods,” Biophysical Journal, vol. 103, no. 11, pp. 2287–2294, 2012.

G. Fischer, M. Handler, P. R. Johnston, and D. Baumgarten, “Impedance and conductivity of bovine myocardium during freezing and thawing at slow rates - implications for cardiac cryo–ablation,” Medical Engineering & Physics, vol. 74, pp. 89–98, 2019.

C. H. Chen, E. A. McCullagh, S. H. Pun et al., “An integrated circuit for simultaneous extracellular electrophysiology recording and optogenetic neural manipulation,” IEEE Transactions on Biomedical Engineering, vol. 64, no. 3, pp. 557–568, 2017.

J. J. Bailey, A. S. Berson, A. Garson Jr. et al., “Recommendations for standardization and specifications in automated electrocardiography: bandwidth and digital signal processing. A report for health professionals by an ad hoc writing group of the Committee on Electrocardiography and Cardiac Electrophysiology of the Council on Clinical Cardiology, American Heart Association,” Circulation, vol. 81, no. 2, pp. 730–739, 1990.

M. H. Jang, S. J. Ahn, J. W. Lee et al., “Validity and reliability of the newly developed surface electromyography device for measuring muscle activity during voluntary isometric contraction,” Computational and Mathematical Methods in Medicine, vol. 2018, Article ID 4068493, 9 pages, 2018.

L. Shaw and S. Bagha, “Online EMG signal analysis for diagnosis of neuromuscular diseases using PCA and PNN,” International Journal of Engineering Science and Technology (IJEST), vol. 4, pp. 4453–4459, 2012.

D. Thys and J. Narang, “Electrocardiographic monitoring,” in Cardiac Anesthesia: Principles and Clinical Practice, F. Estafanous, P. Barash, and J. Reves, Eds., pp. 175–194, Lippincott Williams & Wilkins, 2nd edition, 2001.

A. Neic, M. Liebmann, E. Hoetzl et al., “Accelerating cardiac bidomain simulations using graphics processing units,” IEEE Transactions on Biomedical Engineering, vol. 59, no. 8, pp. 2281–2290, 2012.

E. T. Dougherty, J. C. Turner, and F. Vogel, “Multiscale Coupling of Transcranial Direct Current Stimulation to Neuron Electrodynamics: Modeling the Influence of the Transcranial Electric Field on Neuronal Depolarization,” Computational and Mathematical Methods in Medicine, vol. 2014, Article ID 360179, 14 pages, 2014.

G. Fischer, B. Tilg, R. Modre, F. Hanser, B. Messnarz, and P. Wach, “On modeling the Wilson terminal in the boundary and finite element method,” IEEE Transactions on Biomedical Engineering, vol. 49, no. 3, pp. 217–224, 2002.

J. Haueisen, C. Ramon, M. Eiselt, H. Brauer, and H. Nowak, “Influence of tissue resistivities on neuromagnetic fields and electric potentials studied with a finite element model of the head,” IEEE Transactions on Biomedical Engineering, vol. 44, no. 8, pp. 727–735, 1997.

P. R. Johnston, “The potential for Laplacian maps to solve the inverse problem of electrocardiography,” IEEE Transactions on Biomedical Engineering, vol. 43, no. 4, pp. 384–393, 1996.