Unveiling Stimulation Secrets of Electrical Excitation of Neural Tissue Using a Circuit Probability Theory

Electrical excitation of neural tissue has wide applications, but how electrical stimulation interacts with neural tissue remains to be elucidated. Here, we propose a new theory, named the Circuit-Probability theory, to reveal how this physical interaction happen. The relation between the electrical stimulation input and the neural response can be theoretically calculated. We show that many empirical models, including strength-duration relationship and linear-non-linear-Poisson model, can be theoretically explained, derived, and amended using our theory. Furthermore, this theory can explain the complex non-linear and resonant phenomena and fit in vivo experiment data. In this letter, we validated an entirely new framework to study electrical stimulation on neural tissue, which is to simulate voltage waveforms using a parallel RLC circuit first, and then calculate the excitation probability stochastically.

Here, we propose a new theory, named Circuit-Probability (C-P) theory, to provide a physical framework, which is completely different from the conventional way of using H-H model with E-field modeling. Then, we show that some widelyused empirical models and rules can be intuitively derived from the C-P theory.

METHOD TO OBTAIN THE CIRCUIT-PROBABILITY FRAMEWORK
How should we analyze tissue response to an external stimulation? To answer this question, we performed a thought experiment, which ultimately led to our new framework of Circuit-Probability. When considering the electrode-tissue interaction, the first question is how the electrode is bridged to the tissue. We know the activation of action potential is induced by the gating of the voltage-dependent ion channels. Then, for electrode-tissue interaction, the key issue is how the electrical input affects the voltage on these ion channels. Considering the cell membrane is a capacitor, which is impermeable to ions, it affects its electrical response in two aspects. Firstly, the voltage changes on the capacitor, which is induced by charging and discharging procedures, will generate a different waveform in response to the input waveform. And the charging and discharging procedures are not only affected by the capacitor itself, but also affected by its peripheral circuit. Secondly, the Efield will always be perpendicular to the plate of the capacitor, which is the cell membrane surface, and the direction of the E-field is only determined by the orientation of the capacitor. Apparently, the correct voltage waveform and correct E-field direction can be both obtained with a proper circuit involving the capacitor of cell membrane. This is why we use a circuit to characterize the electric response on the cell membrane.
With a proper circuit, we can model the voltage waveform. Then, from this voltage waveform, how can we know the stimulation strength? In the in vivo testing, the number of activated action potentials shows a continuous change with the electric input. However, single channel measurement shows that an individual ion channel does not display a continuous state change in response to electric input. It acts like a digital bit, which only has two states: closed and open. Meanwhile, the gating pattern of a single ion channel also shows a stochastic behavior. Then, how can we build a bridge from the microscopic discrete to macroscopic continuity? The exclusive option is a probabilistic description of the ion channel gating, just as the situation of thermal dynamics and quantum mechanics. The ion channel gating is stochastic and can be described by the exponential distribution. When a specific electrical field is applied onto the ion channel, the time duration it takes to open the ion channel can be described by the exponential distribution. Then with the voltage waveform simulated using the circuit, it is easy to calculate the probability of activating an action potential with a certain electrical input.
Up to here, we have obtained a basic framework of Circuit-Probability based on pure physical reasoning. The proper circuit configuration can only be fitted using experimental results, which is a posteriori, while the probability equation can obtained by theoretical derivation, which is a priori.

THE METHOD TO OBTAIN THE PROBABILITY CALCULUS IN CIRCUIT-PROBABILITY THEORY
Here we firstly show how to theoretically derive the probability equation.
In the electrical stimulation of a neuron, we assume that electron transition of the protein causes the opening of sodium ion channel, which then generates an action potential (AP). Electron transition is a quantum phenomenon, which is random. Hence, the generation of APs can be described with an exponential distribution for quantum event: Here f (λ, t) represents the probability of AP to be generated within a time duration of t. 1 λ is the expected time until AP is generated. Expand the exponential distribution to a general calculus form: then the normal exponential distribution is the special form when λ (t) is a constant. Meanwhile, 1 λ is also a function of the voltage V: We have three electrophysiological considerations for g(V): Consideration 1: 1 λ to be infinitely large when the voltage, V, is more positive than the threshold voltage, V Threshold . In this condition, the AP cannot be generated. Consideration 2: 1 λ to be inversely proportional with the amplitude of |V − V Threshold |, when V is more negative than V Threshold . Consideration 3: 1 λ to approach a minimal level when |V − V Threshold | goes to infinite. So 1 λ should get saturated at a certain value. With these three considerations, one possible form of g(V) can be expressed as: The equation can be re-written as: Here, α, β, n, and c are adjustment parameters, where α > 0, β > 0, n > 1, and 0 ≤ c ≤ 1.
To simplify the equation, here we assume that n = 1 and c = 0.
Then the complete expression of λ is: Considering the voltage waveform, V(t), is a function of time, t, the complete probability calculus equation is: (7) In this equation, α, β, and V Threshold are three parameters to be determined by data fitting. For a specific voltage waveform as shown in Figure 1A, the voltage waveform can be converted to a λ waveform as shown in Figure 1B. Then the probability calculus can be further simplified as: where S λ is the area of the λ waveform. A detailed analysis of the probability calculus can be found in Supplementary S1.

METHOD TO OBTAIN A CORRECT NEURAL CIRCUIT IN CIRCUIT-PROBABILITY THEORY
Then, we build a proper circuit using the results shown in Figure 2. Its general configuration and analysis can be obtained by reasoning (Supplementary S2). Based on the general configuration, its exact configuration is obtained by fitting the experiment data. This is a parallel RLC circuit. The capacitor refers to the cell membrane. The inductor is included to explain the frequency dependent response observed in the experiments. We validated this by applying a single-frequency input (sinusoid wave) to the Common Peroneal (CP) nerve. Four different current amplitudes were used at different frequencies, spanning from 500-9,000 Hz; (D) the corresponding modeling results showing the local minima and maxima predicted by the C-P theory; (E) a detailed probability mapping showing how the shape of the probability curve changes from low current, which exhibits the resonance effect, to high current, which has monotonically decreasing trend.
Sine wave currents (red curve in Figure 2B) were applied upon the CP nerve to activate the Tibialis Anterior (TA) muscle and the resulting force was recorded. The force measured with respect to frequency forms a curve, here named as 'force mapping curve' in this study. With a specific current waveform, a resultant voltage waveform on the capacitor can be calculated, shown as the blue curve in Figure 2B for probability calculation. Similarly, a probability curve with respect to frequency calculated by modeling is defined as a probability mapping curve. The detailed experiment procedure and testing setup can be found in the Supplementary S3. The force mapping results (force generated by TA muscle) against the pulse width of single pulses (in Hz) of four different current amplitudes curves are shown in Figure 2C.
The same data plotted with error bar can be found in the Supplementary S6.1.
The shapes of these four curves are quite different, showing a complicated changing trend with increasing current amplitude. For the curves of 20 and 40 µA, a clear resonance effect can be observed. However, 80 µA curve shows an initial decline, before increasing to a resonance frequency. The curve of 200 µA shows a monotonically decreasing trend without the resonance effect. Despite these variations, C-P theory can still reproduce the general shapes of the curves via probability mapping (Figures 2D,E shows a more detailed probability mapping). The parameters for the circuit and probability calculus can be found in Table 1-1(d&e). It clearly shows how the forcefrequency curve changes from one shape to another shape with increasing current amplitude over a variety of pulse frequencies and accurately predicts the trend, particularly the existence of local minima and maxima. The probability mapping from the C-P theory reproduces the complex changing trends of the testing Frontiers in Computational Neuroscience | www.frontiersin.org

SIMULATION FITTING TO EXPERIMENTAL RESULTS BY CIRCUIT-PROBABILITY THEORY
The C-P framework and the probability calculus equation is achieved by reasoning, which is a priori, rather than a posteriori. This is very unusual for biological research. Meanwhile, the circuit is still of a preliminary configuration.
To validate the correctness of this priori theory, a series of experiments on four types of non-neural and neural tissues using a rat model were conducted: the skeletal muscles (Supplementary S6.3), the sciatic nerve (Supplementary S6.4), the cortex (Supplementary S6.5), and the pelvic nerve (Supplementary S6.6). All the testing data can be well-fitted or explained by the C-P theory: 1. Different current waveforms will generate force mapping curves with different shapes; 2. Force mapping curves generated by arbitrary current waveforms can be fitted by modeling of C-P theory; 3. The resonance frequency widely exists in nervous systems and can be measured with proper stimulation parameters. To help readers understand how various force mapping patterns are generated and affected by parameters, a general demonstration (Supplementary S5) and a detailed case analysis (Supplementary S4) of how the circuit parameters affect the probability mapping pattern are provided.
Meanwhile, C-P theory can give a unique prediction: the electrical voltage response by electrical stimulation, which is conventionally considered as the stimulus artifact, can be wellfitted by the voltage response of the circuit in Figure 2A. This voltage response will show the same voltage response as a parallel RLC circuit. The data by experiment and modeling can be found in Supplementary S6.2.

THEORETICAL EXPLANATION TO STRENGTH-DURATION RELATIONSHIP AND LNP MODEL BY CIRCUIT-PROBABILITY THEORY
This C-P theory provides a physical understanding of the electrical nerve stimulation, which is not available in previous theories and models. Thus, most of the phenomenological models and theories can be directly derived or even amended from C-P theory. Here we just show how to derive and correct two well-known phenomenological models in electrical nerve stimulation: strength-duration relationship (Lapicque, 1909;Weiss, 1990) and LNP model (Linear-Nonlinear-Poisson cascade model) (Schwartz et al., 2006).
Firstly, we will derive and amend the strength-duration relationship. Previously, it is widely believed that charge is the factor to induce nerve stimulation. In such charge based theory, there is an empirical linear relationship between the threshold charge level and pulse duration, which is called Weiss's strengthduration equation (Weiss, 1990) for negative monophasic square current pulses. This equation describes the threshold charge as a function of pulse width as follows: where I rh is the rheobase current, T ch is the chronaxie, and PW is the pulse width. The rheobase current is defined as the threshold current for infinitely long pulses. The chronaxie is defined as the pulse duration required for excitation when the current amplitude is equal to twice the rheobase current. And Lapicque reiterated Weiss's equation for the strength-duration relationship (Lapicque, 1909), but in terms of the threshold current, and introduced the rheobase current and chronaxie as the constants: Apparently, these two equations are just mathematical descriptions without explaining how I rh happen and why the curve follows a specific trend. As follows is the derivation of this relationship with physical definition of I rh . Figure 3A shows a typical voltage waveform by applying negative monophasic square current with difference SPPW (single phase pulse width). For the voltage waveform of each SPPW, the peak voltage is denoted as V P , which is a function of I and SPPW and written as V P (I, SPPW). Based on C-P theory, nerve excitation can be induced when V P (I, SPPW) ≥ V Threshold . Then both the threshold current I th and the threshold charge, Q th = I th ×SPPW, are defined as the current and charge required to make the V P reaches V Threshold .
Then the critical condition is: I th and Q th can be written as functions of SPPW and V Threshold : (13) Since V P (I th , SPPW) cannot be expressed analytically, only numerical solutions of I th and Q th , which are calculated with   1-3(b,c)] are provided in Figures 3B,C. In Figure 3B, all curves decrease to a constant value, I rh . This is because the V P will saturate at a maximum value, V P max , when SPPW ≥ SPPW P max , as shown in Figure 3A.
The physical meaning of I rh is the threshold current when V P max = V Threshold . Meanwhile, the non-linear curve of I th vs. SPPW, existence of I rh and linear curve of Q th vs. SPPW, can be directly obtained without any additional hypotheses. The exact analytical equation for this relationship is not available. The curves in Figures 3B,C are the numerical solution of strength-duration relationship. It corrects the relationship in two aspects: 1. Rather than infinitely approaching to the I rh as the case in Weiss's strength-duration equation, the threshold current curve will be equal to the I rh when SPPW ≥ SPPW P max . 2. Rather than being a completely straight line, the threshold charge curve is linear only when SPPW ≥ SPPW P max . When the SPPW is approaching zero, the slope of threshold charge curve will also approach zero, meaning that the threshold charge will converge in a constant value at low SPPW.
These two major special differences with the Weiss's equation have already be confirmed by previous research (Friedli and Meyer, 1984;Su et al., 2017) and now can be well-explained in the C-P theory. Moreover, it also explains why such relationship can only be applied for negative monophasic square current waveform. Because the voltage waveforms differs with the current waveforms, inducing a more complicated trend without a stable I rh , which was observed in other researches (Friedli and Meyer, 1984). In Figure 4, representative strength-duration curves of other waveforms including different types of square waves and sine waves are shown. For the curve of sinewave current, the threshold current curve increases at high SPPW range, this phenomenon has been observed in previous research with triangle current waveform (Rodríguez-Fernández et al., 2016). But these curves also vary with different circuit parameters.
Then, we will derive the LNP model. The LNP model is a simplified functional model of neural spike responses (Schwartz et al., 2006). It has been successfully used to describe the response characteristics of neurons in early sensory pathways, especially the visual system. The LNP model is generally implicit when using reverse correlation or the spike-triggered average to characterize neural responses with white-noise stimuli. The number of action potential generated can be described by the Poisson distribution in LNP model.
Actually, the Poisson distribution and exponential distribution describe the same stochastic process. If the Poisson distribution provides an appropriate description of the number of the occurrences per interval of time, then the exponential distribution will provide a description of the time interval between occurrences.
The Poisson distribution is as follow: P x = k; λ is the probability of the k times occurrences of the event in a unit time interval, λ is the expected times of occurrence. The exponential distribution is as follow: FIGURE 4 | (A) The relationship between the threshold current amplitude (I th ) and the SPPW for different current waveforms; (B) the relationship between threshold charge (Q th ) and SPPW for different current waveforms. P (t; λ) is the probability of the occurrences of the event with the time interval t, 1 λ is the expected time interval. These two distributions share the same λ. Apparently, in the C-P theory, if the generation of action potential can be described by exponential distribution, it surely can be described by Poisson distribution.
As follow is the derivation of LNP model. The white noise involved in LNP model can be simplified as a triangle wave series of frequency f and amplitude V w as shown in Figure 5A. Actually, any kind of periodical voltage waveform can be used. The triangle wave is used as an example of simple waveform.
Only part of the voltage can exceed the V Threshold . As explained in Figure 1B, the voltage curve can be converted to a λ curve as shown in Figure 5B. The area S λ of the λ curve within a period T = 1 f can be calculated. Since the λ implemented in the C-P theory is not a constant value while λ in Poisson distribution can only be a constant value, an equivalent λ e for Poisson distribution can be calculated based on the S λ : which is the blue straight line in Figure 5B. Apparently, the λ curve and the λ e curve are of the same area, so they will induce the same statistical results. So the probability calculus equation can be rewritten as: The corresponding Poisson distribution is: By increasing the noise amplitude V w , S λ will also increase, result in an increasing λ e as shown in Figures 5C,D. Since S λ is a function of V w , and λ e is a function of S λ , λ e is also a function of V w , shown as the non-linear curve in Figure 5E. This explains how a linear increment of V w induces a non-linear increment of λ e happened in LNP model. Because the expression of S λ is a piecewise function of V w , the exact function λ e (V w ) can only be calculated numerically with a fixed α, β, V Threshold , and f . The analytical expression of λ e (V w ) is not available.

SUMMARY
In summary, we propose a new theory, named the Circuit-Probability theory, to unveil the "secret" of electrical nerve stimulation, essentially explain the non-linear and resonant phenomena observed when nerves are electrically stimulated. In this theory, an inductor is involved in the neural circuit model for the explanation of frequency dependent response. Furthermore, predicted response to varied stimulation strength is calculated stochastically. Two empirical models, strengthduration relationship and LNP model, can be theoretically derived from C-P theory. This theory is shown to explain the complex non-linear interactions in electrical nerve stimulation and fit in vivo experiment data on stimulation-responses of many nerve experiments. As such, the C-P theory should be able to guide novel experiments and more importantly, offer an in-depth physical understanding of the neural tissue. More detailed discussion about possible issues, such as the inductance and probability calculus, can be found in Supplementary S7. As a promising neural theory, we can even further explore the more accurate circuit configuration and probability equation to better describe the electrical stimulation of neural tissues in the future.