Simple model to incorporate statistical noise based on a modified hodgkin-huxley approach for external electrical field driven neural responses

Noise activity is known to affect neural networks, enhance the system response to weak external signals, and lead to stochastic resonance phenomenon that can effectively amplify signals in nonlinear systems. In most treatments, channel noise has been modeled based on multi-state Markov descriptions or the use stochastic differential equation models. Here we probe a computationally simple approach based on a minor modification of the traditional Hodgkin-Huxley approach to embed noise in neural response. Results obtained from numerous simulations with different excitation frequencies and noise amplitudes for the action potential firing show very good agreement with output obtained from well-established models. Furthermore, results from the Mann–Whitney U Test reveal a statistically insignificant difference. The distribution of the time interval between successive potential spikes obtained from this simple approach compared very well with the results of complicated Fox and Lu type methods at much reduced computational cost. This present method could also possibly be applied to the analysis of spatial variations and/or differences in characteristics of random incident electromagnetic signals.


Introduction
Noise activity and oscillations in the brain are ubiquitous and many studies have probed their effects on the stochastic dynamics of neuronal ensembles [1][2][3][4][5][6][7][8][9][10][11][12][13].It is well accepted that noise can affect many complex behaviors such as stochastic resonance, synchronization, and deterministic chaos.It can also bring about changes in rhythm (e.g., variability in respiratory activity [14]), and enhance the response of a nonlinear dynamical neural system to weak external inputs.Channel noise has been studied in the context of electrical stimulation of the auditory nerve by cochlear implants [15,16] and hippocampal neurons [17].Noise and voltage fluctuations in excitable cells can also affect spike generation [18,19], propagation along axons [20], and spontaneously generate action potentials in the absence of external stimulation [21].
The role of noise in enhancing the system response to weak external signals has been observed not just in biological systems [22], but also in other fields involving physics and chemistry [23][24][25].As first discovered by Benzi and co-workers [26] and known as stochastic resonance (SR), noise was shown to amplify or enhance weak signals in nonlinear systems.Though not studied in great detail thus far, such processes could also conceivably be at work in the context of weak incident signals in an electromagnetic environment.On the other hand, the role and utility of moderate to strong external electrical inputs for biomedical applications, especially in the context of excitable cells, has long been recognized [27][28][29][30][31].For example, nerve conduction block induced by electrical stimulation was studied to reduce urethral pressure during micturition towards treatment of chronic pain, or to stop unwanted motor activity, such as muscle spasms, spasticity, and choreas.External electrical stimulation presents an effective, fast, and reversible route to nerve conduction blockage that would have many advantages over clinical methods for surgically or pharmacologically blocking nerve impulses.For instance, non-electrical methods suffer from issues such as nonspecificity, low success rates and potential nerve destruction.In this overall context, it may also be mentioned that the use of external electromagnetic radiation (whether from microwave sources, cellphones, or other wireless devices) can give rise to bioheating, and membrane temperature increases have been predicted to produce hyperpolarizing currents [32].Such ion currents, coupled with changes in membrane channel conductivities, could then contribute to transient shifts in transmembrane potentials and lead to changes in the firing rates of excitable cells.In other research focusing on the interaction of external radiofrequency fields with biological systems [33,34], thermoelastic expansion and build-up of acoustic pressure waves leading to microwave hearing [35], have been reported and discussed.
Traditionally on the modeling side, the Hodgkin-Huxley (HH) deterministic approach has widely been used to model the initiation and propagation of action potentials in neurons.Action potentials (essentially electrical spikes) generated by the opening and closing of ion-channels embedded in the neural cell membrane, are thought to carry information via their timing, as opposed to their magnitude or the duration.At the single-cell level, stochastic behavior in neurons and variability in action potential timing arises from the random gating of ion channel populations.Such stochastic noise associated with the random opening and closing of individual ion channels can alter the firing threshold [36], spike timing [19,37], interspike interval statistics [38], and the amount of stochastic resonance [39].Based on the HH approach, channel noise was first modeled as a system of gating variables consistent with a multi-state Markov process description [40], with other numerical studies providing insights into fluctuations in discrete ion channel populations [41][42][43].However, the Markov chain (MC) algorithms are numerically expensive, since a large number of channels (thousands or more [44]) are needed to realistically simulate the stochastic behavior.The reason for requiring a large number of coupled ion-channels arises from the essence of the MC models, wherein each state of an ion channel is described by a discrete-state, continuous-time MC and represents a particular configuration of the ion channel.The Markov property requires that a channel's transition from one state to the next depends on its present/current state alone, and the transition rates for any state of the channel are dictated by the membrane potential.As a consequence, all channels are coupled due to their common dependence on the membrane potential and need to be considered collectively.
To overcome the computational complexity of the MC methods, there have been many research efforts, starting from the work of Fox and Lu [45,46], to use models based on stochastic differential equations (SDEs) for including channel noise.Computationally, the SDE approach can be orders of magnitude faster than the MC model [47].However, simulation studies have shown that the SDE approach does not accurately replicate the MC results [48,49].More recently, Goldwyn et al [50] demonstrated that MC approaches that describe the states of ion-channels could capture the distribution of channel noise and its effects on spiking in a Hodgkin-Huxley neural model, while the MC methods that were based on describing the states of ion channel subunits, were unable to accurately represent noise and its downstream consequences.
The aim of the present contribution is to attempt a simple modeling approach for inclusion of noise in the context of the HH model [51] which remains a popular approach for neural analysis.Thus, the approach relies on slight modification of the familiar HH model by adding noise terms.One aim is to determine whether such a modified HH model can result in generation of action potentials similar to those obtained from state space models with noise [52].Furthermore, the present contribution attempts to ascertain if the simple modified HH approach can yield noise generated action potentials (APs), and if so, whether the average number of APs are statistically comparable to those from more complicated Markov chain models.It is shown later in the results section that if the noise parameter is kept within a specified range, the modified HH model described here can mimic the results of the more complicated state-space approaches.Finally, in the context of Markov chains, it may be mentioned that strictly, use of such MC models is applicable to memoryless, stationary systems.Yet neural signaling in electrically excitable cells and the triggering of action potentials can have inherent memory, as has been described in terms of memristors [53,54].The memristors are electrical elements whose resistances depend on the history of electrical voltages or charges, that have crossed it.A more recent development has even been the notion of quantum memristors [55] with decoherence mechanisms controlled by a continuous-feedback scheme that accounts for the memory.Additionally, in the context of external signals, the bio-system is likely to become non-stationary as it is driven by external stimuli.For example, molecular dynamics simulations and experimental studies suggest that terahertz (THz) or mid-infrared electromagnetic radiation may appreciably alter the permeability properties of ion channels [56][57][58][59].Electric fields or mechanical stress are also known to cause membrane poration [31,60,61] that alters the permeability.In such situations, the Markov-based models may no longer be appropriate for an accurate analysis.
For completeness it may be stated that other simpler models, such as the leaky integrate and fire (LIF), could have been used instead of the HH approach.Such models rely on the use of a summation (or integration) process, combined with a mechanism that triggers action potentials above some critical voltage.The APs are then reduced to 'events' that happen at a precise moment in time, without regard to the shape of an action potential.They have been shown to serve as a useful route to practically implement electronic device-based circuits to mimic neural networks and the human brain functionalities [62].The HH model used here, on the other hand, perhaps places our analysis on a firmer footing and is more reflective of the many approaches and studies reported in the literature on neural response behavior.

Materials and methods
The present analysis of neural firing and action potential launches was based on the original onedimensional, distributed HH cable model [51].The numerical implementation was similar to that discussed in the literature [63][64][65][66][67][68].Thus, the membrane potential of the HH neuron was given by the following expression: where m(t), h(t), and n(t) are three gating variables.In equation (1), V(t) is the time-dependent membrane potential, C represents the membrane capacitance per unit area, G K , G Na , and G L denote the maximal potassium, sodium, and leakage conductances per unit area, respectively, E K , E Na , and E L denote the corresponding reversal potentials, and I(t) is the external periodic driving current.The generally accepted values for the various constants were chosen, with C = 1 μF cm −2 , G Na = 120 ms cm −2 , E Na = 115 mV, G K = 36 ms cm −2 , E K = 12 mV, G L = 0.3 ms cm −2 , and E L = 10.6 mV.
Conventionally, the gating variables m(t), n(t) and h(t) are taken to obey three well-known voltage-controlled, time-dependent Langevin equations as first outlined by Hodgkin and Huxley [51].Here the evolution of the gating variables was slightly modified to include random changes.This was carried out in a two-step manner.First, the conventional approach was continued as the governing set of equation, explicitly given in eqns. ( where r m denotes a uniformly generated random number between [−1, 1], while Δ is an amplitude whose value had to be chosen to yield strong agreement with the state space models.In order to restrict g (t + dt) within the [0, 1] range, a new random number (r m ) was generated to recalculate equation (2d) in case the gating variable exceeded unity or had a negative value.In this manner, the noise term was naturally added in through equation (2d).The equation set above was solved using an iterative finite difference backward Euler scheme as reported in the literature [69].
The six voltage-dependent opening and closing rate functions α i and β i have been experimentally determined and given explicitly by the following expressions: Thus here, fluctuations in the gating variables within this scheme were taken to be independent random variables with a zero-mean value.The addition of random noise is similar to the idea first proposed by Fox and Lu [45] where statistically independent Gaussian white noise with zero mean was included.The Box-Muller approach can easily be used for such calculations to generate random noise [70], since it is difficult to obtain the inverse of a Gaussian cumulative distribution function in a quick and accurate manner.However, using a Gaussian random variable implies that some negative values of noise could emerge, requiring successive generation of random numbers to assure that gating variables appropriately remain within the [0, 1] range.Such addition of noise within a differential equation framework of the HH approach can be incorporated either in the current expression, or as subunit noise associated with the gating variables, or as conductance noise [71].Here, the noise was assumed to be embedded in the gating variables.Inherent to equation (2d) are the max and min functions which preclude any possibility of the gating variables (m, h, or n) taking on negative values or values greater than one.In this context, it may be mentioned that various values of the time step (dt) were tested to gauge the likelihood of the gating variables taking on any negative values over long simulations.It was found that a choice of dt less than 10 −7 s was sufficient to avoid any such issues, and so this 0.1 μs value was used here.
The choice for the amplitude Δ in equation (2d), on the other hand, was obtained by running many simulations and comparing the predictions of this modified HH model against the results of one of the accepted models.Specifically, the published stochastic version of the HH model [46,52] was used as a comparative test against results obtained from the proposed modified HH approach.Inclusion of random noise in the stochastic version was based on the Box-Muller approach to generate Gaussian noise (Δn) from uniformly distributed random numbers (r a and r b ) over the range [0, 1] in the following way [46]: where N i denotes the density of ion channels per unit area, A is the area of a membrane patch with values available in the literature.For the Pu-Thomas (PT) implementation [52] of the Fox and Lu model [45], the Box-Muller approach and the actual numerical code was obtained from their Github repository [72].For comparisons between the proposed modified HH approach and the Fox and Lu type schemes [46,52], the Mann-Whitney U test [73] was applied by running many simulations with different assigned Δ values for various input excitation scenarios.The objective was to extract those values of the amplitude Δ that yielded statistically insignificant results between the two schemes.In addition, the number of action potentials launched upon running simulations for a long time for the proposed strategy and a representative Markov Chain (MC) scheme were also compared.The results obtained are discussed at length in the next section, and conclusively show that values of Δ in the range 2.2 × 10 −3 < Δ < 2.3 × 10 −3 yield very good agreement between the present fast and efficient approach, and the more complicated MC models.Finally, as a separate measure of robustness of the present implementation, the distributions of the interspike time durations from the present model were compared to those obtained from the PT formulation.Close agreement was obtained, thus bolstering the case for use of the simple strategy for noise implementation as discussed here.
In terms of computations, some crude comparative estimates can be made between the Pu-type models and the modified HH implementation discussed here.The Pu-type models involve matrix multiplications and singular value decomposition.On the other hand, the modified HH scheme uses multiplication and exponentiation.Based on this simple consideration, the modified HH approach is computationally more efficient, and also allows for the inclusion of noise and random stochastic changes with reasonable accuracy.

Results and discussion
Simulation results comparing the modified HH model, discussed here with the PT implementation of the Fox and Lu model without any external current stimuli, are presented and discussed first.Figure 1 shows the membrane potential [V(t) as given in equation (1)] as a function of time over a 300 ms span with the various spikes.While an exact match is not possible given the stochastic nature of the AP firing, the following features are evident from the results of figure 1.First, despite there not being any excitatory inputs, the launch of three separate AP spikes were predicted by both the modified HH and the Fox and Lu models during the 300 ms interval.Without the addition of noise in the traditional HH model, no action potentials could possibly be triggered despite running simulations for a long time.Furthermore, as might be expected, with inclusion of noise, the AP firings were random without any periodic or repetitive structure.There was, however, a close resemblance between the two sets of results.Other parameters such as the time-dependent evolution of the gating variables [m(t), n(t), h(t)] were also obtained and displayed closely matching behavior, but these have not been shown here for brevity.
Next, the simulations were repeated for a number of current inputs to gauge the responses and compare them to the outputs of the HH and the PT implementation of the Fox and Lu models.Again, very good agreements were obtained, and for concreteness, a representative result with a 100 Hz excitation input is shown in figures 2(a)-(c).The noise amplitude for the HH model [as defined in equation (2d)] was taken to be 2.25 × 10 −3 .Of these various results, figure 2(a) for the time-dependent potentials show a total of six AP launches over a 100 ms time frame for both the modified HH and the Fox and Lu models.Thus, the same number of neural firings were obtained for both cases.Curves for the associated time-dependent gating variables are given in figure 2(b) for the modified H scheme, and in figure 2(c) for the Fox and Lu model.Their behaviors are similar and though some minor differences can be seen, the overall performance of the gating variable transitions is very comparable.
For completeness, a few other simulations were also run at different excitation frequencies of the external current.Specifically, results of the timedependent potential for an excitation current of 0.1 μA for two different excitation frequencies of 50 and 500 Hz are shown in figures 3(a) and (b).As shown previously, the number of AP launches predicted by the two models were identical.For example, for the 50 Hz excitation, four AP launches can be seen in figure 3(a).On the other hand, figure 3(b) reveals fifteen AP firings for the 500 Hz excitation over the 300 ms time frame for each of the two models.The amplitude of the random noise in the HH model was again set at 2.25 × 10 −3 .Thus, looking collectively at figures 1 and 3, a fairy robust replication of the results and trends is seen to be achieved by the current modified HH model for a range of excitation frequencies.
Next, simulation tests were carried out to probe the effect of the noise amplitude Δ as given in equation (2d) on the predicted response of the AP firing.A range of different values were used for the amplitude, and more discussion on the changes resulting from variations in Δ were explicitly probed and discussed later based on a Mann-Whitney U-test [73] comparison.However, the U-test does not produce a time-dependent result, and so for completeness, a specific outcome for the potential V(t) is first given in figure 4. Though several other simulations were also carried out, only one result is shown for brevity, since the trends were qualitatively similar.Thus here, simulation results for a 0.1 μA current excitation at 100 Hz for the time-dependent potential V(t) obtained from the modified HH model for two different noise amplitudes (Δ) of 2 × 10 −3 and 2.5 × 10 −3 ; together with results from the PT implementation of the Fox and Lu approach are shown in figure 4. The difference in figure 4(a) between the Δ = 2 × 10 −3 and Δ = 2.5 × 10 −3 is negligible.Both amplitudes were predicted to produce eighteen peaks over the 300 ms timeframe.Marginal shifts in the potential waveform between the two cases can be discerned in the 175-275 ms range.Furthermore, the results from the modified HH mode for both noise amplitudes appear to match the output for V(t) from the Fox and Lu approach as shown in figure 4(b).The results of the PT implementation also produce eighteen AP peaks with an overall shape that is comparable to the HH predictions.There was no input excitation for the first case, while inputs of 0.1 μA were assumed for excitation at four frequencies (50 Hz, 100 Hz, 250 Hz and 500 Hz) studied.Hence, it can be concluded that the modified HH method would not have any significant statistical difference from the results predicted by the Fox and Lu approach, as long as the noise parameter was suitably chosen.
For greater clarity, values of the number of AP firings obtained from the modified HH approach for

Conclusions
Noise in neural networks is pervasive and can affect many complex behaviors such as stochastic resonance, action potential firing rates, synchronization, and even lead to enhancement of weak external signals.Given the presence of external electromagnetic environment which is often noisy and can influence the neural dynamics, it becomes important to develop simple yet robust models which incorporate the stochastic deviations into a nonlinear neural network analysis.Though different models have been proposed, including the use of different ion currents to incorporate bifurcation [74], or the addition of equivalent currents associated with electromagnetic radiation in neuronal loops [75], the detailed mathematics is often quite complicated, or the parameters for the modified circuit are not easily quantified.Additionally, the inclusion of intricate behaviors through multi-state neural descriptions often leads to computationally complex models that may become prohibitive for analyzing large interacting, nonlinear systems in the presence of noise.
Here, we have proposed a relatively simple scheme to incorporate random noise that relies on the wellknown HH model.In a sense, it is a simplified engineering approach, but has been shown to amply mimic the behaviors of more complex Markov Chain models.It also does not require the system to be stationary, and so any dynamic changes such as evolution in the rate constants or shifts in triggering threshold could easily be incorporated without excessive computational burden.A large number of simulations based on our modified HH model with inclusion of noise, were performed at different excitation frequencies and noise amplitudes.Our results from the Whitney-Man U test also showed the present approach to have statistically insignificant difference from the results predicted by established models reported in the literature.It was also shown that the distribution of the time interval between successive potential spikes obtained from the simple modified HH approach compared very well with the results of Fox and Lu type methods at much reduced computational cost.Finally, it may be mentioned, that though not explored in this contribution, the possibility of assigning different noise amplitudes to a network of neurons to mimic the spatial variations and/or differences in characteristics of random incident electromagnetic signals, could be incorporated in an efficient manner.Similarly, effects of ephaptic coupling could be studied, and will be reported elsewhere.

Figure 1 .
Figure 1.Comparison of the time-dependent potentials with the various spikes as predicted by the present modified HH model and the PT implementation of the Fox and Lu approach without any external stimulus for both cases.

Figure 2 .
Figure 2. Results obtained by running the modified H model and the PT implementation of the MC-based Fox and Lu approach with a 0.1 μA stimulus current at a 100 Hertz frequency.(a) The time-dependent potentials from the two models showing a total of six AP launches for both cases during the 100 ms time span, (b) behavior of the gating variables [m(t), n(t), h(t)] for the modified HH model, and (c) predicted evolution of the gating variables for the PT implementation of the Fox and Lu approach.

Figure 3 .
Figure 3.Comparison of the time-dependent potentials as predicted by the modified HH model and the PT implementation of the Fox and Lu approach for excitation frequencies of: (a) 50 hertz, and (b) 500 hertz.

Figure 4 .
Figure 4. Simulation results for a 0.1 μA current excitation at 100 hertz.(a) Time dependent potential V(t) obtained from the modified HH model for two different noise amplitudes (Δ) of 2 × 10 −3 and 2.5 × 10 −3 .(b) Corresponding time dependent potential V(t) obtained from the PT implementation of the Fox and Lu model.

Figure 5 .
Figure 5. Histogram showing the p-values obtained from the Mann-Whitney U Test comparing the modified HH model with the PT implementation of the Fox and Lu approach for amplitudes Δ ranging from 2 × 10 −3 to 2.5 × 10 −3 .As evident from the histogram, the p-values corresponding to the noise parameter Δ in the range 2.15 × 10 −3 < Δ < 2.3 × 10 −3 are above the 0.05 level for all five cases.There was no input excitation for the first case, while inputs of 0.1 μA were assumed for excitation at four frequencies (50 Hz, 100 Hz, 250 Hz and 500 Hz) studied.

Table 1 .
(a): Results obtained from running the modified HH model for 300 ms with eleven different values of the noise parameter (Δ).An input frequency of 100 Hz was used.Fifty simulation were carried out for each case.The average number of AP firings over the fifty simulation runs for each Δ value are given, along with the standard deviation.(b) Results obtained from running the modified HH model for 300 ms with eleven different values of the noise parameter (Δ).An input frequency of 250 Hz was used.Again, fifty simulation were carried out for each case.The average number of AP firings over the fifty simulation runs for each Δ value are given, along with the standard deviation.

Figure 6 .
Figure 6.Histograms showing the distribution of the time interval between successive potential spikes obtained from two different methods.For all cases a 100 Hz, 0.1uA stimulus current was applied as the input.(a) The results from the modified HH model with noise for Δ = 2.25 × 10 −3 , and (b) the inter-spike time interval distribution obtained from the PT implementation of the Fox and Lu approach.