Peripheral Nerve Stimulation Thresholds Based on Waveform Shape and Implications for Guideline Limits

Authors: Gregory B. Gajda

Abstract

The objective of this paper is to derive basic restrictions for induced internal electric field and reference levels for external magnetic flux density for a class of periodic non-sinusoidal waveforms as multiples of the existing limits applicable to sinusoidal waveforms in current exposure standards. The Law of Electrostimulation and the Spatially Extended Nonlinear Node computational model were used to derive peripheral nerve stimulation thresholds of the internal electric field for both non-sinusoidal and sinusoidal waveforms. Threshold ratios (non-sinusoidal to sinusoidal) permitted basic restrictions and reference levels to be derived as multiples of the sinusoidal ones. Intercomparisons of threshold ratios from both models suggest that they are in agreement for flat-topped flux density waveforms with fast rise-times relative to the period but showed a discrepancy for the continuous sinusoid. Results from the computational model were used to establish the threshold ratios used in the conversion. Resulting non-sinusoidal basic restrictions and reference levels were found to have the same functional relationship with frequency as the sinusoidal ones, consisting of two a flat rheobase and a frequency-dependent (basic restriction) or inverse frequency-dependent (reference level) portion that intersects the rheobase at a transition frequency that is waveform-dependent. Above the transition frequency, the non-sinusoidal basic restriction was found to be inversely related to the flux density rise-time, resulting in an increased limit for fast-rising waveforms. The transition frequencies of fast-rising waveforms were found to be lowered relative to the sinusoidal one. Above the same transition frequency, the non-sinusoidal reference level is flat with frequency and was found to be approximately 79% lower than the sinusoidal one.

INTRODUCTION

Exposure guidelines in the low-frequency range up to 10 MHz exist to prevent undesirable or painful stimulation of peripheral nerves (peripheral nerve stimulation or PNS). Limits are on the external magnetic flux density field strength (referred to as a reference level or RL) and the internal dose metric (DM), which is the induced or in-situ electric field (termed a basic restriction or BR) (ICNIRP 2010). The RL and BR values are usually given in terms of the temporal root-mean-square (RMS) amplitude under the implicit assumption of a sinusoidal magnetic flux density waveform, and it is usually assumed that the spatial distribution of the external flux density is uniform. Because PNS thresholds are a function of the peak amplitude of a waveform, comparisons with the limits should be made based on the peak. Throughout the remainder of this study, all amplitudes were referenced to temporal peak values.

In general, PNS thresholds are different for non-sinusoidal waveforms as compared to sinusoidal ones (Reilly et al. 1985; Havel et al. 1997). The existing guideline BRs and RLs are based on these sinusoidal thresholds, which have been extensively characterized and vetted (IEEE 2019; ICNIRP 2010). Compliance assessment methods for non-sinusoidal and multi-frequency fields have been proposed (e.g., Jokela 2000; Reilly and Diamant 2002; Heinrich 2007). Of the different ones developed, the most capable are methods based on the peak of a spectrally weighted waveform in the frequency-domain (weighted by the inverse RL in terms of amplitude and phase with frequency) that is converted back into the time domain (De Santis et al. 2013). These methods are suited for highly complex or multi-frequency fields. If the weighting function is implemented as a multi-pole filter, a time domain implementation can be incorporated directly into field measuring instruments (Keller 2017).

The purpose of this study was to develop simple, easily implemented compliance criteria for a class of periodic flux-density waveforms that can be described as “flat-topped” (i.e., square waves) with either ramped (linear) or exponential leading edges (Fig. 1a). (The latter waveform occurs when an inductive coil is excited with a square-wave voltage.) Although not expressly developed for the more complex field waveforms encountered in some applications, the criteria developed in this work are simple to implement for the waveforms investigated and give insight into aspects of the non-sinusoidal field that are important from a compliance viewpoint.

This is achieved by determining the BRs of the waveforms in question, which are assumed to be a constant fraction of their stimulation thresholds. The non-sinusoidal BRs are expressed as multipliers of the existing sinusoidal BRs, where the ratio of the two is assumed to be equal to the ratio of non-sinusoidal to sinusoidal stimulation thresholds, both at the same cyclic frequency, f. Once the non-sinusoidal BRs are established, the associated non-sinusoidal RL is determined using the known relationship between the external flux density and the internal dose electric field. The non-sinusoidal RL is expressed in terms of the sinusoidal one.

Two different approaches were used to determine the relative stimulation thresholds for the non-sinusoidal and sinusoidal waveforms. They are the Law of Electrostimulation (LOE) (Irnich and Hebranck 2009) and Spatially Extended Nonlinear Node (SENN) models (Reilly et al. 1985; Reilly and Diamant 2011). The former is an empirical model consisting of a pair of mathematical conditions that must be satisfied for stimulation to occur. In it, the thresholds for an arbitrarily shaped pulse of the internal electric field are related to two parameters in the hyperbolic strength-duration the chronaxie time and rheobase electric field threshold. In contrast, SENN is a computational model of a myelinated nerve fiber exposed to an in-situ electric field along its axis. The SENN model was used to calculate sinusoidal thresholds, followed by thresholds for the non-sinusoidal waveform of interest, from which a threshold ratio (at the same frequency) can be formed.

Existing sinusoidal limits

The induced tissue electric field is spatially non-uniform, even for a uniform external flux density (Kavet 2015). Because it is spatially non-uniform, the dose metric (DM) is the highest spatial value of the electric field vector occurring in the peripheral nerve tissue. The time function of the DM is denoted as E(t) and is related to the time-rate-of-change of the magnetic flux density through Faraday’s law.

The relationship between the DM and flux density can be determined using electromagnetic analysis methods such as scalar potential finite difference (Dimbylow 1998; Dawson et al. 1999), finite difference time domain (Gandhi and Chen 1992), or finite element methods (Wang and Eisenberg 1994) used in conjunction with anatomically detailed, heterogeneous conductivity human models (phantoms). The temporal peak amplitude of the DM, designated as Epk = peak{E(t)}, is related to the external uniform flux density B(t) through a scalar constant, termed the induction or coupling coefficient Ci. Based on Faraday’s induction law, peak DM is given by (Kavet 2015):

where peak{B’(t)} is the temporal peak of dB(t)/dt.

For a sinusoidal flux density waveform,

where Bpk is the temporal peak of the external flux density, f = 1/T the frequency, and T the period. The coupling coefficient depends on the geometrical arrangement of the externally applied flux density and body tissues as well as the tissue conductivity distribution but is independent of the waveform shape.

The frequency behavior of sinusoidal PNS thresholds can be characterized in two asymptotic frequency ranges. Up to a specified transition frequency, the asymptotic threshold is flat, whereas beyond that frequency, the asymptote is linearly dependent on frequency. Existing sinusoidal BRs of up to 5 or 10 MHz are typically fractions of the threshold asymptotes in the same two frequency ranges (the fraction being the reciprocal of the assumed safety factor). In terms of the temporal peak DM, the sinusoidal BR, EBR,sine(f), is expressed as

where EBR,rheo is the rheobase value of the sinusoidal BR in units of V m^−1^, and fe is the transition frequency where the two threshold asymptotes intersect.

Sinusoidal magnetic flux density RLs (in terms of their peak amplitude), BRL,sine(f), usually have the same transition frequency and are in the

where BRL,rheo is a rheobase value of the sinusoidal RL.

In principle, the same relationship between the sinusoidal BR and RL as in eqn (2) should exist. Thus, the sinusoidal BR is related to the sinusoidal flux density RL

where the same coupling coefficient in eqn (1) applies.

In either frequency range (i.e., above or below fe), the ratio of the sinusoidal BR to RL is

This ratio has a linear frequency dependency and from eqns (5) and (6), implies that the coupling coefficient is

From eqn (7), the coupling coefficient for the IEEE 2019 guideline (unrestricted exposure to the head or torso) is 0.163 V m^−1^ T^−1^s. According to the ICNIRP 2010 guideline, the coupling coefficient (exposure to the head and body of the general public) is 0.80 V m^−1^ T^−1^s. The large disparity between the two is largely due to the smaller RL in ICNIRP (2010) relative to IEEE 2019 (despite their similar BRs), which is due to the adoption of additional safety factors by the International Commission on Non-Ionizing Radiation Protection (ICNIRP).

METHODS

The periodic B(t) waveforms investigated in this study are shown in Fig. 1a. The two non-sinusoidal waveforms consist of flat-topped or square waves with ramped or exponential leading edges, whereas the third is a continuous sinusoid. They produce trains of biphasic pulses of E(t), as shown in Fig. 1b, where the term “biphasic” refers to pulses occurring in both amplitude directions. (The term “monophasic” is used for a single pulse of E(t) in the positive amplitude direction.) The term “phase duration” was introduced in IEEE (2019) and is defined as the time between zero crossings of the E(t) waveform. In this context, the term “phase” implicitly describes the excursion of E(t) above or below the zero line.

All flux density waveforms are assumed to (1) continue indefinitely in time, (2) have peak amplitude Bpk, and (3) begin in time so that the resulting E(t) waveform forms a complete, positive-going pulse starting at t = 0. All three waveforms can be described analytically in their first half-period, both in terms of flux density and DM, through the coupling coefficient Ci.

It will be seen that the threshold for a train of periodically repeating biphasic pulses is highly dependent on the width of the pulse relative to the half-period (i.e., T/2). For convenience, unitless multipliers can be formed as a measure of the narrowness of the pulse relative to the half-period. For the square-wave flux density with ramped leading edges in Fig. 1a (simply abbreviated to “ramp”), the multiplier

where d is the duration of the ramp; and for the square-wave with exponential leading edges (abbreviated to “exponential”) in Fig. 1a, the multiplier

where τ denotes the exponential time constant. For narrow pulses of E(t) in relation to the half-period, these multipliers become large. In keeping with the terminology defined previously, they will be termed “phase factors” since they are a measure of the duration of a stimulatory phase at a given frequency or period.

The analytical expressions for the magnetic flux density, B(t), of each waveform are given below with the DM, E(t), for each waveform obtained from E(t) = Ci dB(t)/dt. For the sinusoid, the waveforms are given

The flux density and E(t) for the ramp in the first half-period from t = 0 to t = d is given

The flux density and E(t) for the exponential in the first half period are given

The analytical relationships between B(t) and E(t) in eqns (10) through (12) allow the threshold values of Bpk to be determined from the thresholds of E(t) computed from the SENN and LOE analyses. The peak amplitudes of the E(t) pulses are shown as unity in Fig. 1b because this is the starting amplitude in the SENN threshold calculations. However, conversion factors relating E(t) and Bpk are given for each waveform in the Figure.

SENN model

The SENN model is available for download as a Fortran executable file and associated input/output files.^2^ The model, as downloaded, pertains to a 20-μm-diameter myelinated nerve with 51 nodes of Ranvier, all of which are set as nonlinear. It allows determination of electric field thresholds that result in a propagating action potential through a process described in Reilly and Diamant (2011) as a “titration.” The executable file has a range of built-in waveforms, such as continuous sinusoids, and allows for arbitrary, user-defined waveforms in the form of a time-sampled, waveform input file.

For all calculations reported on, the default nerve-fiber physical parameters that were downloaded were used. (The default fiber parameters intranodal gap = 2.5 μm, membrane capacitance per unit area = 2.0 μF cm^−2^, membrane conductance per unit area = 30.365 mS cm^−2^, axoplasm resistivity = 110.0 Ωcm and medium resistivity = 300.0 Ωcm.) When the SENN-derived thresholds for a single rectangular pulse stimulus are modeled as a hyperbolic strength-duration response, the default parameters result in a chronaxie time of 0.1248 ms (this procedure will be described subsequently). This value of chronaxie will be shown to approximately align with the nominal 3-kHz transition frequency in exposure standards.

LOE model

The LOE model for a pulse of E(t) consists of two simultaneous conditions relating the pulse amplitude and duration at threshold (Irnich and Hebranck 2009):

where ERh is the rheobase electric field and τe is the chronaxie time constant. The times t1 and t2 (t1 < t2) are the time instants where E(t) crosses the rheobase line, and ts is the duration where E(t) (at threshold amplitude) exceeds the rheobase. In effect, it is the portion of the pulse that contributes to stimulation (i.e., generation of an action potential). For brevity, it is referred to as the “stimulus duration” but should not be confused with the temporal dynamics of the action potential, which occurs sometime after the actual stimulus (Reilly 1998). Thresholds found using the LOE model are dependent on assumed values of the rheobase and chronaxie, which is discussed in more detail later.

Examples of single monophasic rectangular, exponential and sine pulses are shown in Fig. 2 to illustrate how the LOE model is applied. In the case of an exponential pulse, t1 = 0 and ts = t2. For the rectangular pulse, t1 = 0, and ts is the entire pulse width. Because the two LOE conditions apply at threshold, they form a system of two equations in the two unknowns; the pulse amplitude at threshold, denoted as ETh, and the stimulus duration, ts.

For trains of biphasic pulses, the LOE model is applied in the same manner as for a single monophasic pulse; that is, it is applied to the first pulse in the train (Irnich and Hebrank 2009). For the pulse trains in Fig. 1b, it is applied to the pulse appearing in the first half-period. For this reason, it is helpful to examine application of the LOE to single monophasic pulses because it results in the same sets of equations. For a rectangular pulse (shown in Fig. 2) of width d, the hyperbolic condition (eqn 13) results in the well-known hyperbolic strength-duration (SD)

where ETh is the peak amplitude of the pulse at threshold and ts = d. With the threshold amplitude given by eqn 15, the rheobase condition is automatically satisfied.

Application of the LOE to a single sine pulse (shown in Fig. 2) with a half-period of the same duration, d, as the rectangular pulse results in the two equations (details of the derivation can be found in Gajda and Bly 2017):

In this case, the rheobase condition of eqn (14) is applied as an equality, resulting in the right-hand equation in eqn (16).

The two equations in eqn (16) form a system in the two unknowns, ETh and ts, as a function of the parameters ERh and τe (those parameters define the hyperbolic relationship of a rectangular pulse given by eqn 15). Factoring out ERh results in a transcendental equation in the unknown ts that can be solved numerically. Once ts is known as a function of d, it can be inserted into the right-hand equation of eqn (16) to solve for the threshold ETh. The threshold amplitude vs. d forms the SD relationship for the sine pulse. An important feature of the LOE model is that the same values of rheobase and chronaxie apply to all pulse shapes. Thus, if the rheobase and chronaxie can be determined from the SD curve for a rectangular pulse, in principle, the SD curve for any arbitrarily shaped pulse can be obtained.

The single exponential pulse in Fig. 2 is given by E(t) = Eo exp(−t/τ), where τ is the time constant, and the pulse extends to infinity. The LOE applied to the exponential pulse (with the rheobase condition applied as an equality) results in the set of

where ETh is the threshold value of Eo.

RESULTS

Single pulses

The SD curves for rectangular, sine, and exponential pulses are plotted in Fig. 3. For all three LOE SD curves, a rheobase of 5.90 V m^−1^ and a chronaxie of 0.1248 ms were used. These two hyperbolic parameters were obtained from curve fitting the results of a SENN analysis of a single rectangular pulse (the results of which are also plotted in Fig. 3). The curve fitting was carried out using a non-linear least squares method, and the resulting RMS relative error between the SENN threshold data and the curve fit to eqn (15) was 1.2%.

The results of the SENN threshold calculations are compared with the LOE-derived SD curves for the sine and exponential pulses in Fig. 4 (for the SENN analysis of the exponential pulse, a total time base of 5τ was used). The close correspondence between the LOE and SENN SD curves in Fig. 4 suggests that the LOE and SENN models are in close agreement for single monophasic pulses.

The correspondence between the SENN and LOE-predicted thresholds is particularly close for short pulse widths or time constants. In this range, the stimulus duration ts approaches d for the sine pulse and τ for the exponential pulse (Gajda and Bly 2017). For d < < τe, the LOE sine threshold approaches ETh,sine = (π/2) ERh τe/d (this can be observed by letting ts → d in the hyperbolic equality of eqn 16). For a rectangular pulse [ramp of B(t)], the threshold approaches ETh,ramp = ERh τe/d. Thus, the ratio of the thresholds for sine and rectangular pulses for short pulse durations is π/2.

For a single exponential pulse, the asymptotic LOE threshold for τ < < τe approaches ETh,exp = ERh τe/τ (which can be shown by using the approximation exp(−ts/τ) ≈ 1 − ts/τ and letting ts → τ in the hyperbolic equality of eqn 17). Thus, when plotted against τ in Fig. 3, the exponential threshold converges to the asymptotic curve for the rectangular pulse (i.e., for d < < τe). For the exponential pulse, the time constant, τ, is equivalent to the pulse duration, d, for a rectangular pulse in terms of its effect on asymptotic thresholds (τ can be considered the “effective” pulse or phase duration).

For long pulse durations (i.e., for d or τ > > τe), the LOE SD curves in Fig. 3 all converge to the rheobase ERh. For a sine pulse in this asymptotic range, the stimulus duration ts becomes a small fraction of the pulse duration and is centered on the peak (Gajda and Bly 2017). Under the condition ts/d → 0, the rheobase condition in eqn (16) yields ETh = ERh because cos(πts/2d) ≈ 1. A similar conclusion is reached for the exponential pulse when ts/τ → 0, that is, exp(−ts/τ) ≈ 1, and the rheobase condition of eqn (17) becomes ETh = ERh. This asymptotic threshold behavior can be generalized to any pulse shape with a distinct peak. If the peak is high enough to satisfy the rheobase condition, then the brief stimulus duration will occur only around the peak. The threshold amplitude will equal the rheobase regardless of the detailed shape of the waveform.

A notable difference between the two models for the sine pulse at long pulse durations is the upturn in the SENN threshold with increasing pulse duration. This behavior is present in many other SENN analyses of sinusoids and is attributed to the slow rate of change of the sine pulse (in contrast to the immediate rise in the rectangular and exponential pulses, which is independent of the pulse width) (Reilly et al. 1985; Reilly 1998). In view of the upturn, the long-pulse asymptote (or rheobase) can be considered as the horizontal tangent to this part of the SD curve.

Continuous waveforms

The magnetic flux density waveforms for which exposure standards are applied are typically continuous and periodic. They produce internal electric fields consisting of alternating or biphasic pulses or phases. For periodic waveforms, the threshold amplitude can be plotted against frequency to produce a strength-frequency (SF) curve. This curve generally has two asymptotes, one flat with frequency in the lower frequency range and the other having a linear frequency dependence in the higher frequency range (this is illustrated in Fig. 5). The threshold transition frequency is the frequency at which the two asymptotes intersect.

The periodic waveforms that were the subject of this study differed from the single monophasic pulses described previously. First, the periodic sinusoidal pulses in Fig. 1b extend over the entire half-period T/2, whereas the periodic ramp and exponential pulses (also shown in Fig. 1b) are assumed to be fractions of the half-period (i.e., α > 1 and β > 1). In addition, the periodic version of the exponential pulse terminates at each half-period before it flips in polarity. Its leading edge of B(t) (Fig. 1a) is assumed to be short relative to the half-period, resulting in a time constant τ < < T/2. When the leading edges are short in duration, the E(t) pulses of the flat-topped flux-density waveforms are separated by a zero interval or a dead space. This is not present in the sinusoid, where the positive phase is immediately followed by the negative phase. This is an important feature that will be shown to affect the agreement between the LOE and SENN models.

Because the LOE analysis is carried out on the first positive phase of each E(t) waveform (Irnich and Hebranck 2009), the defining LOE equalities for the sinusoid, ramp, and exponential are given by eqns (16), (15), and (17), respectively. (Note that for the continuous sinusoid, d is replaced with T/2 in eqn 16.) However, because the focus is on the relative thresholds and BRs of these three waveforms, their asymptotic, linear frequency-dependent thresholds are of primary interest (since all thresholds converge to the same rheobase at low frequencies). Relative thresholds are characterized in this study in the form of threshold ratios (e.g., the ratio of asymptotic threshold of ramp to sinusoid). In this way, the rheobase ERh and chronaxie τe do not explicitly appear in the LOE threshold ratios because they cancel through division.

The thresholds from the LOE analysis of any of the waveforms have two asymptotes that cross over at the threshold transition frequency. The asymptotic threshold for the ramp, ETh,ramp(f), can be expressed

where ERh is the rheobase value appearing in the LOE equalities (eqns 15, 16, or 17), and fTh,ramp is the ramp threshold transition frequency. Like the case of single monophasic pulses, the rheobase threshold for continuous periodic waveforms, as in eqn (18), is identical for all continuous periodic waveforms.

For a continuous sinusoid, the high-frequency asymptotic LOE threshold approaches the quantity (π/2) ERh τe/(T/2). When written in the form ETh,sine(f) = ERh f/fe, the sinusoidal transition frequency fe from the LOE analysis is related to the chronaxie

Because this transition frequency is explicitly given in the guideline, knowledge of the chronaxie is unnecessary. The IEEE (2019) and ICNIRP (2010) guidelines list values of fe as 3,350 Hz and 3,000 Hz, respectively.

The high-frequency asymptotic threshold for the ramp is the quantity ERh τe/d, which can be expressed as ETh,ramp(f) = ERh f/fTh,ramp. This gives a threshold transition frequency for the ramp in terms of fe

The high-frequency asymptotic threshold for the exponential (for τ < < T/2) is the quantity ERh τe/τ, which can be expressed as ETh,exp(f) = ERh f/fTh,exp. Thus, its transition frequency is given in terms of fe

The SENN model was used to calculate the thresholds of all three E(t) waveforms in Fig. 1b to establish the rheobase threshold and their high-frequency asymptotic thresholds for a range of phase factors α and β. The high-frequency asymptotes were characterized by a single coefficient (e.g., for the ramp and for f > > fTh,ramp, ETh,ramp(f) = Krampf where Kramp = ERh/fTh,ramp from eqn 18). Similar coefficients were defined for the sinusoid and exponential as ETh,sine(f) = Ksine f and ETh,exp(f) = Kexp f, respectively. In principle, this only requires a single SENN calculation at a frequency at which the threshold is known to display a linear frequency dependency. However, at least two SENN threshold calculations were performed in the linear region to confirm asymptotic behavior.

For SENN calculations, all E(t) waveforms were assigned an initial peak amplitude of 1 V m^−1^, and the titration process allowed to evolve to the threshold peak amplitude. For each frequency, a waveform data file (sampled in time) was created for multiple periods. The number of periods was incremented until the change in the titrated threshold was less than 3%. For frequencies of 50, 20, 10, and 5 kHz, the number of required periods was eight, six, four, and three periods, respectively. At 2 kHz and below, only a single period was necessary to produce the requisite threshold resolution. In the case of the sinusoid, the executable file contained an internally generated waveform where the frequency and number of periods could be set in the input file. This eliminated the need to set up sampled waveform files for the sinusoid.

The results of the SENN-calculated thresholds versus frequency for the ramp (for phase factors of β = 1, 3, 5, and 7) and sinusoid are shown in Fig. 5. In the figure, the computed thresholds are plotted as discrete markers, whereas the rheobase threshold and frequency-dependent asymptotes are plotted as lines. The coefficients of the asymptotic straight lines are obtained from the linear portion of the threshold curve at the high-frequency end. The rheobase line is the tangent of the threshold at lower frequencies. In Fig. 5, the low-frequency thresholds for the four ramp waveforms and sinusoid all appear to converge to a common rheobase, as predicted by the LOE.

Table 1 summarizes the asymptotic data in Fig. 5 for the LOE and SENN models. Because the threshold transition frequency is the intersection of the rheobase with the linear asymptote, the SENN values are calculated from (e.g., ramp): fTh,ramp = ERh/Kramp. The SENN and LOE-derived transition frequencies are given in rows 2 and 3 for comparison. The LOE transition frequencies are based on the value of the SENN sinusoidal transition frequency, fe = 2.64 kHz (since LOE transition frequencies are relative to an assumed value of fe).

The equivalent chronaxie that produces the SENN sinusoidal transition frequency of fe = 2.64 kHz is τe = 0.121 ms (computed using eqn 19). As mentioned previously, the chronaxie found from curve-fitting the SENN-derived SD data for a single monophasic rectangular pulse was τe = 0.1248 ms. The correspondence in chronaxie is close (within 6%) considering the wide differences between the two waveforms. In addition, the rheobase threshold for a single monophasic pulse from the SENN analysis (5.90 V m^−1^) is close in value to the SENN-derived one pertaining to continuous periodic waveforms (6.30 V m^−1^). This is not surprising since only a single period was needed in the titration process to arrive at a convergent threshold for the low frequency range.

The fourth and fifth rows in Table 1 provide asymptotic threshold ratios for the SENN and LOE models, respectively. The former is computed by taking the ratios of the asymptote coefficients for the two waveforms (e.g., for the ramp and sinusoid, Kramp/Ksine). The latter is obtained using high-frequency asymptotic thresholds from the LOE:

where eqns (19) and (20) have been used to replace fe and fTH,ramp, respectively.

The results of the SENN-calculated thresholds versus frequency for the exponential for phase factors α = 1, 5, 10, 15, and 20, are shown in Fig. 6. The lowest phase factor, α = 1, violates the requirement τ < < T/2 but is included to observe the approach to the rheobase for long pulse durations. To economize the effort for the two highest phase factors, α = 15 and 20, SENN threshold calculations were performed only at 20 and 50 kHz, as this was sufficient to establish the asymptote coefficient. For the three lower phase factors, SENN threshold calculations for frequencies as low as 0.1 kHz were made. Below this frequency, limitations on the record length of the waveform file and the minimum required sampling (time) interval precluded the determination of the threshold.

The asymptotic approach to the rheobase can be clearly observed for the lowest phase factor, α = 1. For α = 5 and α = 10, the trend of the data toward the rheobase suggests that it is reached for frequencies well below 0.1 kHz. A summary of the linear asymptote coefficients and a comparison of the SENN and LOE transition frequencies and asymptotic threshold ratios are given in Table 2. The LOE asymptotic threshold ratio was obtained as

where eqns (19) and (21) have been used to replace fe and fTH,exp, respectively.

DISCUSSION

As shown in Fig. 5 and Tables 1 and 2, except for the ramp with β = 1, all other thresholds lie above the rheobase and the linear asymptote for the sinusoid. (A ramp with β = 1 constitutes a triangular wave of flux density since the ramped portion occupies the full half period, i.e., d = *T/*2.) Thus, the sinusoidal BR, which is a multiple of the sinusoidal rheobase threshold and linear asymptote, forms a lower bound for the BRs that would be applicable to the flat-topped waveforms in this study (except for β = 1). In terms of compliance with an exposure guideline, the sinusoidal BR can be considered conservative, except in the case of a triangular flux density waveform. However, as will be seen subsequently, the use of the exact form of the BR is important in the derivation of RLs for flat-topped flux densities in both asymptotic frequency ranges.

An insight into why the high-frequency asymptotic thresholds for the flat-topped waveforms (e.g., ramp), are higher than those for the sinusoid can be seen in Fig. 1b). Here the “pulse” duration for the ramp occupies a shorter portion of the half-period than does the sinusoid, which occupies the entire half-period. From the LOE analysis of single monophasic pulses, it was observed that for short pulse durations (relative to the chronaxie), the product of the stimulus duration and threshold amplitude is equal to a constant (e.g., for the ramp, ETh d = ERh τe). Thus, shorter pulse durations of the ramp or exponential in Fig. 1b result in higher thresholds. In addition, in the short-pulse regime, the peak amplitude at the threshold is closely tied to the total area under the E(t) pulse (in eqn 13) because the stimulus duration ts approaches the entire pulse duration. This is important when considering the effect of the pulse shape.

Using the terminology of IEEE (2019), the sinusoid has a “phase duration” of T/2, whereas the phase duration of the ramp is d (in the high-frequency asymptotic range, the phase duration is equivalent to the stimulus duration ts in the LOE). Thus, for large ramp phase factors, β, the phase duration is much smaller than the half-period for the same cyclic frequency. This would produce a threshold that is proportionately higher than for a waveform with a longer phase duration (e.g., the sinusoid), even though the two waveforms have the same cyclic frequencies. Based on the same reasoning, two waveforms with the same phase duration but different cyclic frequencies should have similar thresholds.

To illustrate this point, consider two ramp one with β = 3 and the other with β = 5. At 33.3 kHz, the ramp with β = 3 had a phase duration of d = 0.005 ms, whereas the ramp with β = 5 had the same phase duration at 20 kHz. From the SENN results in Table 1, the asymptotic threshold for the first ramp is 3.87 × 33.3 = 128.9 V m^−1^ while that of the second is 6.28 × 20 = 125.6 V m^−1^. Thus, the threshold depends on the phase duration rather than on the actual cyclic frequency.

There is also a dependence on the pulse shape in the LOE results. For single monophasic pulses that are much shorter in duration than the chronaxie, the ratio of the rectangular to sinusoidal thresholds for equal pulse widths is 2/π = 0.637. This same result is obtained from the LOE for the periodic ramp and sinusoid whose asymptotic threshold ratio is given by eqn (22), where for equal phase duration of both d = (T/2) or β = 1. For the periodic exponential vs. sinusoid, the LOE threshold ratio is also 2/π in the short-pulse regime since τ is equivalent to the pulse width, d, and α = 1 for an equal phase duration.

The SENN results can be tested for waveforms of equal phase duration in the high-frequency asymptotic range by noting that (e.g., for the ramp) ETh,ramp(f) = Kramp f = Kramp/T and T = 2 β d (from eqn 8). Thus, the ratio of the asymptotic ramp threshold to the sinusoid threshold at equal phase duration (i.e., d = T/2)

The threshold ratio of exponential-to-sinusoid for equal phase duration (i.e., τ = T/2) is given by a formula similar to eqn (24), except that β is replaced by α. The SENN threshold ratio of ramp-to-exponential for an equal phase duration (i.e., d = τ)

Threshold ratios for equal phase duration are computed from the SENN asymptote coefficients in Tables 1 and 2 and are given in Table 3. It shows that the asymptotic threshold ratios for the ramp-to-sinusoid converge to approximately 0.5 at β = 7, while the same is true for the exponential-to-sinusoid for α = 20. The SENN threshold ratios are lower than the ones predicted by the LOE (i.e., 0.637) by 18% and 24%, respectively. Alternatively, the SENN threshold ratio for the ramp-to-exponential (β = 7 and α = 20) approaches the LOE value of unity.

The results in Table 3 suggest that the LOE can accurately predict the thresholds of the two flat-topped waveforms but not for the continuous sinusoid. Assuming that the ramp and exponential thresholds are in agreement between the LOE and SENN models (for large β and α), this implies that the asymptotic thresholds of the sinusoid derived from the SENN analysis exceed the LOE thresholds by a factor of 0.637/0.5 or 4/π. This factor is again observed when non-sinusoidal RLs are derived from non-sinusoidal BRs.

A possible explanation for the difference between the LOE and SENN sinusoid thresholds is the presence of dead space (or flat-line period) between adjacent pulses in the ramp or exponential E(t) waveforms, which does not occur in the continuous sinusoid. This hypothesis was tested by calculating the SENN thresholds for a flat-topped flux density waveform with sinusoidal leading edges (SLEs), as shown in Fig. 7. The pulse duration, d, of the SLE E(t) waveform has the same phase factor, β = (T/2)/d, as the ramp, only in this case the pulse is sinusoidal instead of rectangular. The SENN thresholds were determined at cyclic frequencies of 50 and 20 kHz for a phase factor of β = 7. At 50 kHz, a total of eight periods were required in the waveform file to produce a constant threshold of within 3%.

The asymptotic threshold coefficient obtained from these SENN calculations was KSLE = 13.26 V m^−1^ kHz^−1^. Comparing SLE with the ramp, both with a phase factor of β = 7, results in an asymptotic threshold ratio of ETh,SLE(f)/ETh,ramp(f) = 13.26/8.66 = 1.531. This threshold ratio is very close to the value π/2 = 1.571 predicted by the LOE for the ratio of a sine pulse to a rectangular pulse for the same pulse width. This result and the asymptotic threshold ratio for the ramp-to-exponential in Table 3 imply that the LOE and SENN models are consistent for trains of biphasic E(t) pulses when the pulse is followed by a relatively long dead space. A similar conclusion about the LOE was reached by Irnich and Hebranck (2009), except in their study, the LOE was tested against experimental human threshold data for two flat-topped flux density waveforms, one having ramps and the other SLEs, both with a constant flat-top (or dead space) period of 0.5 ms.

BRs for flat-topped flux density waveforms

The information in Tables 1 to 3 was used to derive the BRs for flat-topped flux densities, starting with the guidance given in IEEE 2019 for the case of non-sinusoidal E(t) waveforms. For the asymptotic linear frequency range, the guidance recommends calculating an effective frequency, feff, which is equal to the inverse of double the phase duration, tp, that is, feff = 1/(2tp). (The phase duration is defined in IEEE 2019 as the interval between zero crossings of a pulse, or, in the case of an exponential, equal to the time constant, τ.) In the IEEE guidance, the effective frequency is inserted into the formula for the sinusoidal BR (eqn 3). Using the example of a ramp with a phase factor β = (T/2)/d, this guidance results in an effective frequency of feff = 1/(2d) = β f where f is the cyclic frequency. Inserting the effective frequency into the sinusoidal BR formula for f > fe results

This procedure results in a ramp-to-sinusoid BR ratio equal to β. From Table 1, it can be seen that the SENN asymptotic threshold ratio is approximately half this amount for large β (e.g., for β = 7, the SENN threshold ratio is 3.63). Because the ratio of the BRs should equal the asymptotic threshold ratio, the IEEE 2019 guidance allows for a ramp BR that is approximately twice the value suggested by the SENN threshold data.

The same procedure, performed for an exponential, results in a BR ratio (for f > fe) of EBR,exp(f)/EBR,sine(f) = α. From the SENN asymptotic threshold ratios in Table 2, this is again seen as an over-allowance by a factor of approximately two times for large α. While correctly identifying the phase duration as the main determinant of the threshold, the IEEE procedure can be revised to align with the SENN-derived threshold ratios (relative to a continuous sinusoid) by halving the resulting BRs. Using the same examples as above, the BRs in the linear asymptotic frequency range under the revised procedure are EBR,ramp(f) = (β/2) EBR,sine(f) and EBR,exp(f) = (α/2) EBR,sine(f), where in both cases, EBR,sine(f) is given by eqn (3) (for f > fe). Additionally, this revised procedure is conservative for small β and α values, as seen from the SENN asymptotic threshold ratios in Tables 1 and 2, which all trend above the one-half value.

An important aspect of the non-sinusoidal BRs derived using the revised IEEE (2019) procedure is that the non-sinusoidal BR transition frequency is no longer equal to the threshold transition frequency. Because all waveforms have the same rheobase thresholds, all rheobase BRs are identical and equal to the sinusoidal one, EBR,rheo. Using the ramp as an example, the BR transition frequency fBR,ramp is the frequency at which the linear frequency part of the BR intersects the rheobase BR and is obtained by equating the two

A similar BR transition frequency is obtained for the exponential, except that β is replaced with α. Thus, the BR for a ramp is given

while a similar BR is obtained for the exponential except with β replaced by α.

Up to this point, the IEEE definition of the phase duration presents no problem for the ramp or exponential (which is considered to be equal to d or τ, respectively). In the case of the SLE, the asymptotic threshold differs from that of the ramp for the same phase factor (i.e., KSLE = 13.26 V m^−1^ kHz^−1^ and Kramp = 8.66 V m^−1^ kHz^−1^, both for β = 7). These different coefficients lead to different BRs in the linear frequency range for the same phase factor and are due to the different pulse shapes of the ramp and the SLE. Thus, the definition of the phase duration requires further refinement for cases where the pulse of E(t) has no zero crossings (as exemplified by an exponential) or has a well-defined pulse duration but is not rectangular (as in the case of the SLE).

An insight can be obtained from the LOE in the high-frequency asymptotic range. In this range, ts < < τe and the integral in the hyperbolic LOE condition (eqn 13) is equal to the constant ERhτe. If, at threshold, the integral is carried out over the entire half-period (i.e., from t = 0 to T/2), then for the ramp, the resulting integration still yields ETh d, where d is the obvious phase duration and ETh is the peak of E(t) at threshold. Similarly, for τ < < T/2, integration of the exponential over the entire half-period results in ETh τ, where τ is known to be the effective phase duration. In both cases, the phase duration (or its effective value), tp,eff, can be determined

where peak{E(t)} is the temporal peak of E(t) during the half-period.

If eqn (29) is applied to the SLE pulse depicted in Fig. 7, the resulting effective phase duration is tp,eff = 2d/π. This is less than the phase duration of the ramp (i.e., d), because of the smaller area under the SLE pulse relative to the rectangular pulse of the ramp. In this example, the effective frequency is feff = π/(4d). For a phase factor of β = (T/2)/d, the effective frequency becomes feff = π β f/2. Substituting this effective frequency into the revised BR formula, the SLE BR becomes (for f > fBR,SLE):

This is equivalent to EBR,SLE(f) = (β/2)(π/2) EBR,sine(f), which, when compared to the BR for a ramp with the same cyclic frequency and pulse duration, is π/2 times higher. This is confirmed by both the SENN and LOE analyses. Finally, eqn 30 permits the transition frequency, fBR,SLE, to be determined

RLs for flat-topped flux density waveforms

The threshold for the magnetic flux density can be derived from the threshold of E(t) and the analytical relationship between B(t) and E(t); that is, eqns (10) through (12). The ratio of flux density RL to its threshold is assumed to be identical to the ratio of the BR to its threshold of E(t). For example, the RL for a ramp waveform BRL,ramp(f), is given

The relationship between peak{E(t)} and Bpk for the ramp (eqn 11) allows the threshold ETh,ramp(f) to be written in terms of the flux density threshold (i.e., the value of Bpk at threshold), BTh,ramp(f),

For f < fBR,ramp, the ramp RL can be written starting with eqn (32), followed by substitution of BTH,ramp(f)/ETH,ramp(f) using eqn (33), substitution of EBR,ramp(f) using eqn (28), and substitution for EBR,rheo using eqn (7) to

where eqn (27) was used to replace fe.

For f > fBR,ramp, the ramp RL is obtained using the same

For the exponential (τ < < T/2), the peak of E(t) occurs at t = 0, and the relationship between ETH,exp(f) and BTH,exp(f) from eqn (12) is identical to that for the ramp (except that β is replaced with α). Using the same steps as for the ramp, the RL for the exponential is the same as eqns (34) and (35), except that the transition frequency is given by fBR,exp = 2fe/α.

For the SLE waveform, the relationship between peak{E(t)} and Bpk at threshold allows the thresholds to be written

where ETh,SLE(f) is the threshold of the DM and BTh,SLE(f) is the associated flux density threshold.

For f < fBR,SLE, the SLE BR is equal to EBR,rheo and using similar steps as followed for the ramp (i.e., using eqns 32, 36, and 7), the SLE RL, BRL,SLE(f),

where eqn (31) is used to replace fe.

For f > fBR,SLE, the SLE BR is given by eqn (30), and the RL is found using steps similar to those used for eqn (37) and is given

These results show that the RLs for the considered flat-topped flux-density waveforms have a similar form to that of the sinusoidal RL. The effective RL rheobase for the flat-topped waveforms is reduced by π/4 (79%) from that of the sinusoidal waveform, and the transition frequency is the BR transition frequency, unique to each waveform. The factor 4/π was previously observed to be the ratio of the SENN asymptotic threshold to that of the LOE for the continuous sinusoid. When only the LOE model is considered (Gajda and Bly 2017), the rheobase RL for flat-topped flux densities (i.e., ramp and exponential) is identical to the sinusoidal one, and the transition frequencies for the BRs and RLs are the same as for the thresholds in eqns (20) and (21).

Even more significant than the small reduction in the effective RL rheobase from that of the sinusoidal RL is the lower RL and BR transition frequencies for flat-topped waveforms. From eqn 27 and its analog for the exponential, a sinusoidal transition frequency of nominally 3 kHz results in a transition frequency of 1.2 kHz for a phase factor (β or α) of 5. Thus, for fast-rising flux density waveforms, the frequency range over which the RL is flat significantly extends downwards.

CONCLUSION

Within the existing framework of sinusoidal limits in standards, the BRs and RLs for a class of flat-topped flux-density waveforms have been derived. Their magnitudes are multiples of the sinusoidal BR and RL in two asymptotic frequency ranges that are delineated by a waveform-specific transition frequency. To accomplish this task, two nerve stimulation threshold models were LOE and SENN.

The computation of relative thresholds (i.e., the ratios of thresholds of one waveform to another) for single monophasic pulses of the DM [i.e., E(t)] using both models demonstrated a close correspondence between the two. For periodic trains of biphasic E(t) pulses, the LOE model was applied to the first pulse in the train as if it were a single isolated pulse, whereas the SENN computations were carried out over a sufficient number of periods to produce an unchanging threshold. Two different flat-topped flux density waveforms were evaluated having ramped and exponential leading edges with varying widths of their leading edges over a wide range of frequencies. The third waveform, evaluated as a function of frequency, was a continuous sinusoid. Finally, thresholds for an additional flat-topped waveform with sinusoidal leading edges were evaluated only in the high-frequency asymptotic range to test a hypothesis.

In both models, the thresholds converged to a common rheobase threshold at low cyclic frequencies. At high cyclic frequencies, thresholds from both models approached asymptotic linear frequency dependence. Threshold ratios between the different waveforms and the sinusoid, in the linear asymptotic frequency range, were used in the derivation of the BRs and RLs.

A multi-way intercomparison of asymptotic threshold ratios between the ramp, exponential, SLE, and continuous sinusoid suggested that the LOE and SENN models are in agreement for flat-topped flux densities, where E(t) is a train of pulses separated by dead spaces. Differences in thresholds relative to the sinusoid between the LOE and SENN models suggested that the former undervalued the threshold of a continuous sinusoid (where there are no dead spaces between pulses) by a factor of 79% (or π/4). The applicability of the LOE for periodic biphasic waveforms with a dead space was also confirmed experimentally by Irnich and Hebranck (2009); however, they did not test it against a waveform with adjoining pulses such as the continuous sinusoid.

Using the threshold ratios obtained from the SENN calculations, BRs applicable to flat-topped waveforms were derived. For this purpose, the guidance found in IEEE (2019) for converting existing sinusoidal BRs to non-sinusoidal waveforms was modified to take into account the SENN thresholds. Similar to the sinusoidal BR, the non-sinusoidal BR is divided into two frequency ranges delineated by a waveform-specific transition frequency. In the low-frequency range, the non-sinusoidal BR is equal to the existing sinusoidal rheobase BR. Above the transition frequency, in the frequency-dependent range, it is equal to a multiple of the sinusoidal BR, which is inversely proportional to the pulse (phase) duration relative to the half-period. This implies that for fast-rising leading edges, the non-sinusoidal BR in the frequency-dependent range is significantly higher than the sinusoidal BR at the same cyclic frequency.

The transition frequencies of the non-sinusoidal BRs were found to be related to the sinusoidal BR transition frequency and inversely proportional to the phase factor. This means that they are proportional to the pulse or phase duration relative to the half-period of the waveform, implying that for fast-rising leading edges, the BR transition frequencies are significantly lower than the sinusoidal transition frequency.

Once the non-sinusoidal BRs were found, the corresponding RLs were derived. The RLs have the same general form as the existing sinusoidal RLs; that is, an inverse frequency dependence at low frequencies and a flat portion (or rheobase) at high frequencies. The transition frequencies dividing the two ranges were identical to that of the non-sinusoidal BRs. The rheobase for the RL was found to be lower than that for continuous sinusoids by a factor of π/4 or 79%.

The information obtained in this study reaffirms the validity of the LOE model as a means of estimating the relative thresholds of E(t) waveforms with relatively long dead spaces (arising from flat-topped flux densities). The procedures for deriving BRs and RLs for such waveforms, as multiples of the sinusoidal BRs and RLs in the standards, can be used as guidance in such standards. Finally, the RLs for non-sinusoidal external electric fields can be similarly adapted using the methods in this work because the internal DM is related to the time-rate-of-change of the external electric field (Kaune and Gillis 1981).