Stochastic analysis of the electromagnetic induction effect on a neuron’s action potential dynamics

The presence of time-varying electromagnetic fields across a neuron cell may cause changes in its electrical characteristics, most notably, in the action potential dynamics. This phenomenon is examined by simulating electrophysiology of a single cortex neuron. Magnetic flux is captured by using a polynomial approximation of a memristor embedded into Hodgkin–Huxley model, equivalent electrical circuit of a neuron cell. Bifurcation analysis is carried out for two different electrical modes associated with the nature of the external neuronal stimulus. Aiming to determine the true influence of the variability of ion channels conductivity, the stochastic sensitivity analysis is undertaken post hoc. Additionally, numerical simulations are enriched with uncertainty quantification, observing values of ion channel conductivity as random variables. The aim of the study is to computationally determine the sensitivity of the action potential dynamics with respect to changes in conductivity of each ion channel so that the future experimental procedures, most often medical treatments, may be adapted to different individuals in various environmental conditions.


Introduction
The non-ionizing electromagnetic radiation and its influence on human tissue and organs is most often observed through two frameworks: low frequency effects on the nervous system (a change in the action potential dynamics of the neuronal membrane) and thus on brain metabolism, and radiofrequency thermal effects. Recently, published guidelines for limiting exposure to electromagnetic fields by International Commission on Non-Ionizing Radiation Protection (ICNIRP) [1] state that the only effect at frequencies of 10 MHz or higher is human or organ heating, though some authors claim that even on such high frequencies, given that the incident power density is sufficiently low, it is possible to observe ultrasonic vibrations in the membrane of a neuron as a result of interaction with electromagnetic radiation [2], which is supported experimentally [3]. However, the most of the literature defines clear boundaries between the influence of low-frequency radiation (up to 100 kHz) and radiofrequency radiation (100 kHz to 300 GHz), where for the range from 100 kHz to 10 MHz, the possibility of influence of neuronal stimulation and tissue heating both exists [4]. In summary, as the radiation frequency increases, heating effects predominate and the likeli-hood of nerve stimulation decreases. A detailed review of neurobehavior during exposure to low-frequency electric and magnetic fields is given in [5]. In this paper, we are concentrated on the effects of time-varying magnetic fields where the main interaction is the Faraday induction of electrical fields and associated currents, the value of which depends on the tissue configuration, orientation and conductivity. The electro-stimulation of a neuron, the electrical polarization of presynaptic processes leading to a change in post synaptic cell activity, is examined using a computational neuron model extended with an additive variable which is capable of capturing the effect of electromagnetic induction.
Extensive research has been carried out w.r.t. magnetic stimulation of a neuron and neuronal network through clinical studies and computer simulations, respectively [6][7][8]. The treatment of the neurodegenerative diseases in the context of electromagnetic induction is also discussed, e.g., for Alzheimer disease [9], and for epilepsy [10]. The effectiveness of electromagnetic stimulation in the treatment of mental illnesses is computationally assessed in [11] through the regulatory abilities of the external electromagnetic stimulation on the pattern dynamics of neuronal networks.
In addition to an external source, the change of concentration of ions in a neuron cell could also lead to the occurrence of electromagnetic induction. This phenomena is inspected in [22], where authors introduce the four-variable Hindmarsh-Rose neuron model with an additive variable used in order to obtain feedback memrisitve current, which directly depends on the variation of magnetic flow. The previously outlined work is further improved in [23], where the electrical neuronal activity mode transition is investigated by applying external electromagnetic radiation on the improved three-variable Hindmarsh-Rose neuron model. Furthermore, the effect of electromagnetic induction, described by the modulation of magnetic flux on the membrane potential of a neuron, realized using memristor coupling together with an additive phase noise imposed on the neuron in order to detect the dynamical response and phase transition mode, is investigated in [24]. Authors in [25] have carried out an extensive review of all relevant literature considering dynamics in a neuron and neuronal network, paying special attention to realization of a computational model under the effect of electromagnetic induction and/or radiation in order to better understand the electromagnetic effect in treating neurodegenerative diseases. The dynamics of the action potential of a neuron under the electromagnetic induction is shown to be greatly dependant on the environmental surrounding temperature in [26] using the extended inductionbased Hodgkin-Huxley neuron model, similar to that used in this paper. The temperature influence is previously verified in the clinical study, where temperature effects are examined in rat hypothalamic tissue slices [27]. The action potential of a neuron largely depends on the stimulus itself, especially if the noise is present [28]. Recently, many studies have been carried out, where both the presence of noise and the electromagnetic induction is observed through dynamical analysis framework for a single neuron [29][30][31] and neural networks [32][33][34]. Oftentimes, realistic neurons hold complex anatomical structure, for instance, autapse connection to some internuncial neurons. The autaptic regulation of electrical activity modes in a neuron under the effect of electromagnetic induction is investigated in [35]. The complexity of a neuron can also be observed through its electrical properties-a neuron is a charged body where any changes in fluctuation of either magnetic flux or electric charge can cause field variation and thus the variation in electrical properties of a neuron. A novel neuron model, able to capture the effect of the total electromagnetic field, is developed and described by dimensionless dynamical system and is verified on analog circuit platform [36]. Another important factor, the energy, is investigated in [37], where authors analyze the role of the energy supply in brain metabolism and signal transmission by using Izhikevich neuron model in the presence of electromagnetic field. Ion channels conductivity is large impact factor on the overall electrical behaviour of a neuron as well [38,39]. A detailed analysis of the hidden bursting patterns and bifurcation mechanisms of a memristive Hindmarsh-Rose neuron model stimulated by a current generated as a result of the external magnetic flux coupling via a memristive circuit element is given in [40]. Authors show that this model, much like a biological neuron, exhibits complex dynamics, which is further analyzed by implementing the proposed model to the hardware breadboards. It is shown that the proposed memristive neuron model mimics the dynamics of the electrical activities of biological neurons with threshold electromagnetic induction. All previously outlined model variants of a single neuron and neuronal network have a common purpose-reliability and interpretability of the output where the greatest attention is usually paid to nonlinear dynamical analysis through bifurcation or phase graphs to determine the influence of various factors (radiation, induction, noise, temperature) on the dynamics of the action potential on the membrane and consequently to gain insight on the coding and information transfer, brain metabolism, potential treatments for neurodegenerative disorders, etc. Extensive overview of bifurcation analysis of neuronal systems is available in [41] where it has been shown how the type of bifurcation determines the neuro-computational properties of the cells. More specifically, a comprehensive bifurcation analysis of a single and coupled Hodgkin-Huxley neuron model with delayed feedback with or without the presence of inhomogeneity is given in [42,43], respectively. However, a large impact of error, i.e., uncertainties, in parameter space and initial conditions is often neglected and all relevant papers published in the field of neural nonlinear dynamics and computational neuroscience, to the best of our knowledge, take the values of input parameters as constants. In this paper, by considering conductivity values of each ion channel in Hodgkin-Huxley neuron model as a random variable uniformly distributed around the expected value, the original deterministic system is recasted to stochastic system [44]. The expected values of ion channels conductivity are taken from the 2017 review study on dynamics of single neuron and neuronal networks [25]. The stochastic collocation method is used as a method of choice for uncertainty quantification because of its non-intrusive nature-no requirements for intrinsic system formulation perturbation and manual interventions are needed. By computationally assessing the sensitivity of the action potential dynamics, the interspike interval (ISI) duration and firing rate w.r.t. the changes in the conductivity of ion channels of a neuron exposed to the electromagnetic induction, experimental procedures with magnetic fields, such as medical treatments of neurodegenerative diseases with transcra-nial magnetic stimulation or deep stimulation [45], will be possible to adapt to different individuals in various environmental conditions. The outline of the paper is as follows: The formulation of the extended Hudgkin-Huxley neuron model is given in Sect. 2, numerical simulations and bifurcation diagrams for two different types of the external neuronal stimulus are given in Sect. 3, the sensitivity analysis of a neuron model electrophysiology by considering ion channels conductivity as random variables is given in Sect. 4, finally, concluding remarks are outlined in Sect. 5.

Formulation
The model of choice in this study is Hodgkin-Huxley neuron model [16], originally described as a set of four coupled nonlinear ordinary differential equations. The greatest significance of Hodgkin-Huxley neuron model is that it is bio-physically meaningful (describing the dynamics of membrane potential under the influence of some external stimulus, the dynamics of ionic currents, synaptic integration, etc.), measurable and readily available for bifurcation analysis [46]. Although Hodgkin-Huxley neuron model is prohibitive and computationally expensive, which additionally causes implementation difficulties and extended simulation time, since the focus of the study is the electromagnetic effect on a single neuron, rather than a neural network, these shortcomings can be overlooked. In order to capture the electromagnetic induction in a neuron cell, an additive variable in a membrane potential equation of Hodgkin-Huxley model is introduced. This variable further extends the original model with an additional equation describing the change in the magnetic flux, ϕ. Hodgkin-Huxley neuron model is then given in the following form [26]: Equation (1) describes the change of the membrane potential, V m , in time, where C m stands for the lipid bi-layer capacitance, andḡ K ,ḡ Na andḡ L are maximum values of potassium, sodium and leakage ion channel conductivity, respectively. E K , E Na and E L are associated reversal potentials for potassium, sodium and leakage ion channel, respectively. I ext represents the external neuro-stimulus by means of the electric current source. The dynamics of gating variables, denoted as y, is outlined in Eq. (2), where α y and β y represent an ion channel gate's opening and closing rate, respectively. Each gating variable, y = n for the potassium ion channel, y = m for the sodium ion channel or y = h for the leakage ion channel, has its activation and inactivation steady-state represented by the Boltzmann equation as a function of the membrane potential, V m , and the environmental temperature, T , as follows: where the function governing the environmental temperature is defined in the following manner: The influence of the environmental temperature is discussed in detail in [26], where it is shown that as the temperature increases, the membrane potential amplitude slightly decreases, while the spike frequency greatly increases. Arising from the Faraday's law in the Maxwell's framework, the presence of a time-varying magnetic field across an electrical conductor will induce a spatially varying non-conservative electric field. Since the neuron cell is an electrical conductor, its dynamics is well-captured by the biologically justified neuron model realized as the equivalent electrical circuit. Thus, the presence of a time-varying magnetic field across the neuron cell give raise to the electric field. This electric field causes a change in the membrane currents resulting in depolarization (or hyper-polarization) of neurons [47]. The feedback current of magnetic flux is realized by parallel placement of the memristor in the equivalent electric circuit as in Fig. 1.
Memristor (abbreviation of memory resistor) is a nonlinear resistor with memory introduced in [48] as the forth fundamental electric circuit element which compensates for the missing link between the charge, where M(q) is the memristance (abbreviation of memristor's resistance) measured in . Memristor, memristive systems and possible applications are further elaborated in [49]. The first physical realization of the memristor is achieved by constructing a nanoscale-thin film device which is equivalent to a time-dependant resistor whose value at time t is linearly proportional to the amount of charge that passed through it up until t [50]. This device can be modelled as serially connected doped resistor, R on , and undoped resistor, R off , and its iv characteristic is given in [51] as follows: R 0 stands for the total resistance at time t = 0 and is formulated as: where w 0 is the initial length of the doped region and D is the total length of the memristor. The polarity of the memristor can be either η = +1, when the dopant drift expands the doped region, or η = −1, when the contraction of the doped region occurs. Furthermore, ΔR is the difference in resistance of two regions: Finally, Q 0 is the total charge value that is required to pass through the memristor for the dopant boundary to move through the entire length of the memristor, D: 1 Hodgkin-Huxley neuron model equivalent electrical circuit [16] with the magnetic flux-controlled memristor, used to bridge the gap between the magnetic flux and the membrane potential and to enable the trans-membrane feedback current flow where the mobility of a dopant is represented as μ.
Equation (12) describes the memristor as a purely dissipative electric circuit element in a sense that it does not store, supply or transmit energy. However, in order to simplify the numerical implementation, distance ourselves from the physical realization of the memristor, and ensure seamless simulation process, we assume memristor as a flux-controlled, smooth, continuous, cubic and monotone-increasing function as presented in [52]: where a and b are arbitrary positive parameters. The memductance (abbreviation of memristor's conductance), measured in S, is defined as: and is embedded in Eq. (1) after being scaled by factor k associated with the suppression modulation on the membrane potential [22]. The product term kρ(ϕ)V m is the additive induction current posed on the neuron membrane. Parameters a and b of a nonlinear polynomial approximation of the memristor are chosen to be 0.4 and 0.02 as in [26]. To match the physical unit of the additive memristive current with the rest of Eq. (1), the parameters a and b should be set to have the physical units of C/Wb and C/Wb 3 , respectively. To physically justify the choice of a and b, free parameters in Eq. (12) are fitted via Nelder-Mead simplex algorithm [53] and the compatibility between hysteresis loops of the nonlinear polynomial approximation of the memristor and the physical realization of the memristor is shown in Fig. 2.
The full list of free parameters, details of the fitting procedure and the implementation nuances are elaborated in "Appendix A".
Finally, Eq. (3) in Hodgkin-Huxley neuron model outlines the dynamics of the magnetic flux, ϕ, where the term k 1 V m stands for the membrane potential-induced flux changes and the term k 2 ϕ describes the flux leakage [23].

Numerical simulations
Expected maximum values of ion channel conductivity are set as in [54] as follows:ḡ K = 36.0 mS/cm 2 , g Na = 120.0 mS/cm 2 andḡ L = 0.3 mS/cm 2 , with the coefficient of variation (CV), the percentage used to create the upper and lower boundary, i.e., the range of uniform distribution for parameters of interest, of 5, 10 and 20 per cent, respectively. The rest of the parameters and the initial condition setup of Hodgkin-Huxley neuron model are given in Table 1.
Numerical solution for all calculations in the study of Eqs. (1)-(3) is carried out by using the MATLAB implementation of nonstiff differential equation solver based on the fourth-order Runge-Kutta integration method [55].
The simulation of an action potential dynamics is carried out within two different electrical modes in Fig. 2 Steady periodic orbits display compatibility between the nonlinear polynomial approximation of the memristor [52] and the actual physical realization [51]   the following two subsections. The first mode, outlined in Sect. 3.1, is so-called tonic spiking electrical mode [46], where the action potential, represented as a time series of discrete, nearly identical membrane potential spikes, encodes information and carries it throughout a neuron's cell. The distribution of ISIs is close to an exponential distribution [56], where the total entropy is found to be ∼ 2.0. On the other hand, the second mode, outlined in Sect. 3.2, is so-called tonic bursting electrical mode [46], where periods of rapid action potential spiking (bursting) are followed by quiescent periods typically of longer duration than the mean ISI. Bursting usually occurs either depending on the type of a neuron and its intrinsic configuration or depending on the type of the stimulus. This study examines the case of the latter-bursting is the response of a neuron w.r.t. the external noisy stimulus current. However, it is important to note that by changing the ion channels' conductivity in repeated simulations, the internal con-figuration of a neuron is also changed, and thus different types of activity may occur.

dc input current-tonic spiking electrical mode
The tonic spiking behaviour is achieved by using the external stimulus direct current (dc) set to I ext (t) = 10 µA/cm 2 for a complete simulation duration of Δt sim = 1500 ms. Induction coefficient, k, is set to a constant value of k = 0.1, while the parameter in the potential-based change of the magnetic flux term in Eq. (3) is set to k 1 = 0.001 and the parameter in the magnetic flux leakage term in Eq. (3) is set to k 2 = 0.01. This kind of tonic spiking behaviour represents the optimal behaviour of a neuron exposed to the electromagnetic induction effect, where optimality is viewed as the ability to transfer maximum information, expressed through the entropy measure. Even though, in computational neuroscience literature, entropy is expressed in different forms for a variety of measured data, more details on the topic can be found elsewhere, e.g., [57], we assume the entropy definition as in the original work by Shannon [58]: where observable x i is the ith ISI in the simulation. For the previously outlined simulation configuration, action potential dynamics along with the histogram of ISIs, captured post-simulation, are shown in Fig. 3. The effect of the electromagnetic induction is investigated via bifurcation analysis where the induction coefficient, k, is taken as the bifurcation parameter. Analysis consisted of a repeated simulation of a fixedconfiguration Hodgkin-Huxley neuron model at a constant temperature of T = 10 • C for the induction coefficient ranging from 0 to 2. Instead of displaying all ISIs, the expected value of ISI, mISI, is evaluated. Furthermore, by using a stochastic collocation method, details of which are available in "Appendix B", the output uncertainty by means of 95 per cent confidence interval is determined. The uncertainty comes from the fact that we observe the conductivity of each ion channel in Hodgkin-Huxley neuron model as a random variable uniformly distributed around its expected value. The stochastic analysis is performed by using 3, 5, 7 and 9 collocation points in each stochastic dimension, thus leading to the total of 27, 125, 343 and 729 deterministic simulations respectively. Collocation level is changed from 3 to 9 in order to test the convergence of the stochastic collocation method and, as shown in Fig. 4, the convergence is satisfying for 5 collocation points, which is used for further stochastic analysis. Therefore, a total of 125 deterministic simulations is performed in order to carry out the uncertainty quantification and sensitivity analysis by means of stochastic collocation method.
Mean ISI over range of k ∈ [0, 2] in simulation duration of Δt sim = 300 ms using 5 collocation points for CV = 5, 10 and 20 per cent is shown in Fig. 5a-c.
All three scenarios, depicted in Fig. 5a-c, show the same exact behavior of a neuron exposed to the effect of electromagnetic induction: an increase of the induction coefficient leads to a decrease in neuronal activity with the notable values of k in range from k ∼ 1 to k ∼ 1.5, in which the greatest influence of the uncertainty of ion channels conductivity is detectable, after which the transition from a spiking state to a quiescent state occurs. It is also shown that the uncertainty for this critical region, where neuron changes its activity mode, is unchanged no matter what the value of CV is.
To test this even further, the stochastic analysis on the spike frequency bifurcation diagram, where k is again defined as the bifurcation parameter, is performed. Number of collocation points is changed from 3 to 9 in order to test the convergence of the stochastic collocation method and it is found to be at the satisfying level for 5 collocation points. Spike frequency (firing rate per 1 s) over range of k ∈ [0, 5] in simulation duration of Δt sim = 1000 ms using 5 collocation points for CV = 5, 10 and 20 per cent is shown in Fig. 6a-c. Figure 6 clearly indicates that the change in electrical activity of a neuron occurs for the same values of k as in the previous simulation, shown in Fig. 5, with no change in uncertainty for this critical region.

Noisy input current
The tonic bursting behaviour is achieved by using the noisy periodic current set to I ext (t) = 10 sin(ωt) + w.n. µA/cm 2 as an external stimulus of a neuron for a complete simulation duration of Δt sim = 1500 ms. The angular frequency, ω, is set to ω = 100/Δt sim Hz, and w.n. stands for white noise. Since the input is characterized by a large signal-to-noise ratio, compared to the ideal input for the tonic spiking electrical mode, the distribution of ISIs is less like an exponential distribution with, as expected, lower value of the entropy. All parameters of Hodgkin-Huxley neuron model along with induction parameters are set as in the previous simulation in which the tonic spiking is achieved. The simulation is depicted in Fig. 7.
Similar to the previous subsection, bifurcation analysis of mean ISI is performed in order to investigate the effect of electromagnetic induction on a neuron, in this case, stimulated by a noisy periodic current. In this simulation, environmental temperature is set to T = 6.3 • C. Induction coefficient, ranging from 0 to 5 is again set as the bifurcation parameter. Remaining parameters are not changed and the simulation is performed for 3, 5, 7 and 9 collocation points, respectively. Convergence graph is shown in Fig. 8.
Mean ISI over range of k ∈ [0, 5] in simulation duration of Δt sim = 300 ms using 5 collocation points for CV = 5, 10 and 20 is shown in Fig. 9a-c. Figure 9a [59] c show the exact behavior of a neuron exposed to the effect of electromagnetic induction: an increase of the induction coefficient leads to a decrease in neuronal activity. Here, a neuron undergoes two activity transitions. The first transition occurs at k ∼ 0.5, where a bursting state remains unchanged up until k ∼ 0.88. The second transition occurs after a sharp decrease of expected ISI when the neuron passes into a quiescent state at k ∼ 1.5.
For the sake of completeness, the stochastic analysis of spike frequency bifurcation diagram is also performed. Spike frequency over range of k ∈ [0, 5] in simulation duration of Δt sim = 1000 ms using 5 collocation points for CV = 5, 10 and 20 is shown in Fig. 10a-c.

Sensitivity analysis
Sensitivity analysis ranks input parameters from the most to the least significant ones with respect to the impact their variability has on the total variance of the output of interest. In this work, sensitivity analysis is carried out according to a so-called analysis of variance (ANOVA) principle [60]. ANOVA is based on a Hoeffding function decomposition and computation of Sobol indices [61]. Namely, the total variance due to the variability of all random input parameters is decomposed into terms corresponding to variances due to variability present in all possible subsets of random input parameters. For details on the method, see "Appendix C".
ANOVA sensitivity analysis results are presented in form of the first-and higher-order sensitivity indices. The sensitivity analysis in this work is threedimensional from stochastic point of view since random input parameters are:ḡ Na ,ḡ K andḡ L . Thus, there are three first-order sensitivity indices, S 1 (ḡ Na ), S 1 (ḡ K ) and S 1 (ḡ L ), each corresponding to the impact ofḡ Na ,ḡ K andḡ L . Higher-order sensitivity indices reflect the impact of mutual interactions of input parameters. There are three second-order sensitivity indices, S 2 (ḡ Na ,ḡ K ), S 2 (ḡ Na ,ḡ L ) and S 2 (ḡ K ,ḡ L ), reflecting     the impact of the interaction between the associated conductivity subsets. Higher stochastic dimensions yield even higher-order sensitivity indices, which is caused by the growth of the number of subsets corresponding to mutual interactions of input parameters.
If stochastic dimension is too large, a total effect sensitivity index may be computed in order to avoid the computation of higher-order sensitivity indices to ultimately reduce the computation time. The total effect index measures the contribution to the output variance of a single input parameter including all variances caused by its interactions, of any order, with any other input variables. For example, the total effect index S T (ḡ Na ) measures the total effect sensitivity thatḡ Na has, together with its interaction toḡ K andḡ L , on both the mean ISI and spike frequency distribution over k. If the first order and total effect indices have virtually the same value, the mutual interactions have a negligible impact on the variation of the output value.
Sensitivity indices for mean ISI distribution over k are shown in Fig. 11 for the case of tonic spiking electrical mode, and in Fig. 12 for the case of tonic bursting electrical mode.
Both modes are exposed to the electromagnetic induction effect where induction coefficient is in range [0, 2] and [0, 5] for the tonic spiking and tonic bursting electrical mode, respectively.
Interestingly, for CV = 5 and CV = 10 per cent, the mean ISI hard drop in the critical range of k values for the tonic spiking electrical mode is affected by the variability in the interaction of the sodium and potassium ion channels, and the leakage channel variability exerts practically no impact. For CV = 20 per cent, the output, that is mean ISI over k in range [0, 2], is most sensitive w.r.t. the variability in sodium ion channel from the starting point of the critical induction coefficient range, For all CVs, the variability of the sodium channel beginning at k ∼ 1 onward has the largest impact on the mean ISI, for the tonic bursting electrical mode. For k < 1, the first order and the total effect sensitivity indices revolve around similar values and the significance of the impact is difficult to discern. Both sodium and potassium ion channel, as well as the respective second-order combinations, achieve quite similar behavior and influence on the output.
Virtually the same observation can be reported for the sensitivity indices for spike frequency distribution over k for the case of tonic spiking electrical mode and  Fig. 11 Variance-based sensitivity analysis for mean ISI distribution over k, where external neuronal stimulus is constant dc current. Upper row holds the first and the second-order sensitivity indices of the output givenḡ Na ,ḡ K andḡ L respectively, for the case of 5 collocation points. Lower row holds the total effect sensitivity indices of the output. Each column shows simulations for a different coefficient of variation   Fig. 12 Variance-based sensitivity analysis for mean ISI distribution over k, where external neuronal stimulus is noisy periodic current. Upper row holds the first-and the second-order sensitivity indices of the output givenḡ Na ,ḡ K andḡ L respectively, for the case of 5 collocation points. Lower row holds the total effect sensitivity indices of the output. Each column shows simulations for a different coefficient of variation tonic bursting electrical mode shown in Figs. 13 and 14, respectively.

Conclusion
An extended variant of Hodgkin-Huxley neuron model under the electromagnetic induction is presented. In the presence of an induced electric field governed by the Faraday's law of induction, associated electrical feedback current is embedded in the model using a theoretical model of a flux-controlled memristor. Both numerical simulation and stochastic sensitivity analysis are performed in order to assess the influence of the uncertainties on the output, where the output has been represented as the dependence of both the mean inter-spike Fig. 14 Variance-based sensitivity analysis for spike frequency distribution over k, where external neuronal stimulus is noisy periodic current. Upper row holds the first-and the second-order sensitivity indices of the output givenḡ Na ,ḡ K andḡ L respec-tively, for the case of 5 collocation points. Lower row holds the total effect sensitivity indices of the output. Each column shows simulations for a different coefficient of variation interval and spike frequency for various induction coefficient values, which ultimately indicate the strength of the magnetic field to which the isolated neuron has been exposed. The uncertainty of the output arises from the propagated uncertainties of the conductivity of ion channels, often neglected in the literature. It is shown that for the characteristic range of induction coefficient values at which a neuron realizes a transition in electrical behavior -from spiking/bursting to quiescent state, no change in uncertainty occurs, regardless of the prior distribution of parameters defined as random variables. Additionally, stochastic sensitivity analysis is performed, where it is shown that the output depends mostly on the interaction of variability in the conductivity of sodium and potassium ion channel (second order sensitivity). Insights gained from this paper can greatly contribute to the development of experimental applica-tions of induction for therapeutic purposes, especially in medical treatments of neurodegenerative disorders.

A Fitting the iv characteristic of the memristor
To physically justify the choice of a and b in Eq. (17), parameters D, μ and R on are set to be fitted. Parameter R off is assumed to be 20 times greater than R on and parameter w is assumed to be 0.2 times the value of the parameter D. The physical memristor's polarity is assumed to be η = −1. The formulation for the memristance of the physical memristor is defined as follows: where R 0 , ΔR and Q 0 are given in Eqs. (13)- (15). Loss function is defined as the L 2 -norm between the memristance of the actual physical realization Eq. (A.1) and the memristance of the nonlinear approximation given as the reciprocal of Eq. (17) for various time-dependent charge values. The optimization is performed using the Python scientific module SciPy [62] or, more precisely, its implementation of the Nelder-Mead iterative minimization method [63]. The convergence graph is shown in Fig. 2, where free parameters have taken on the following values: R on = 1.589 · 10 −1 μ = 6.237 · 10 −14 m 2 Vs D = 1.698 · 10 −7 m and are consistent with the available experimental data [50].

B Stochastic moments computation by stochastic collocation method
Stochastic collocation (SC) method is a non-intrusive method for uncertainty quantification of the output of interest for mathematical models whose input parameters exhibit random nature. The non-intrusive property implies that the formulation and computational code are not interfere with, rather the model is used in a black-box manner. Much like in Monte Carlo methods, SC is used for pre/post processing of a certain number of deterministic simulations, but it requires significantly less deterministic simulations, and thus can be regarded as computationally effective.
The principle of SC method lies in the polynomial approximation of the output of interest, Y , over the domain defined by input parameters modelled as random variables (RVs). Assuming that d input parameters exhibiting random nature are organized in a vector X = [X 1 , X 2 , . . . , X d ], the polynomial approximation for Y is expressed as follows: where (X) is a multivariate basis function and Y i is the outcome of the ith out of n deterministic computations for Y . The polynomial expression in Eq. (B.1) is plugged into well-known formulas for stochastic moments given in the integral form [64] and the final expressions are given as follows: -Expected value (expectation or mean): Here, ω i denotes the weight of the ith multidimensional simulation point and is computed as: where pd f (X) = k pd f (X k ) it the joint probability density function [65].
Combining different types of basis functions and integration rules for computation of ω results in different variants of SC method. Furthermore, tensor or sparse grids of the univariate basis functions are used to construct multivariate basis functions [66]. In this paper a full tensor grid with univariate Lagrange basis functions and Gauss-Legendre integration rule is chosen. Namely, previous work demonstrated that Lagrange basis functions are good choice for interpolation in SC approach [67].

C Analysis of variance (ANOVA) approach for sensitivity analysis
Sensitivity analysis (SA) is a discipline with an aim to rank input parameters from the most to the least significant ones with regards to the impact certain parameter's variability has on the total output variance. Let us denote the total number of input random variables as d. A very simple approach in SA is a one-at-a-time (OaT) approach, in which d univariate stochastic scenarios are considered and the corresponding variances of the output are compared [68]. Another approach is the ANOVA in which Sobol indices arise from the following formulas [61]: -First-order sensitivity index, S 1 (X k ), gives information on the impact of the kth input parameter [60] as follows: , k = 1, 2, . . . , d.
(C.1) -Second-order sensitivity index gives information on the impact of the interaction between the ith and the jth input parameters: , k = 1, 2, . . . , d.
(C.2) -Total effect sensitivity index determines overall impact of the kth input parameter along with its own interaction with other input parameters: In Eqs. (C.1)-(C.3) the term Y |X (.) X ∼(.) stands for a conditional expectation where the tilde symbol, ∼, is read "all except". If S 1 (X k ) and S t (X k ) have the same value, the interaction between input parameters does not impact the overall output variance. Note also that all variances and expectations in Eqs. (C.1)-(C.3) are computed from the same set of deterministic simulations that are utilized in uncertainty quantification as outlined in appendix B, i.e., the stochastic collocation method is used for both the uncertainty quantification and the sensitivity analysis.