Uncertainty Propagation in Nerve Impulses Through the Action Potential Mechanism

Journal of Mathematical Neuroscience (2015) 5:3

DOI 10.1186/2190-8567-5-3

S H O RT R E P O RT Open Access

Aldemar Torres Valderrama Jeroen Witteveen

Maria Navarro Joke Blom

Received: 2 September 2014 / Accepted: 27 November 2014 / Published online: 12 January 2015 2015 A. Torres Valderrama et al.; licensee Springer. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0

Web End =http://creativecommons.org/licenses/by/2.0 ), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract We investigate the propagation of probabilistic uncertainty through the action potential mechanism in nerve cells. Using the HodgkinHuxley (H-H) model and Stochastic Collocation on Sparse Grids, we obtain an accurate probabilistic interpretation of the deterministic dynamics of the transmembrane potential and gating variables. Using Sobol indices, out of the 11 uncertain parameters in the H-H model, we unravel two main uncertainty sources, which account for more than 90 % of the uctuations in neuronal responses, and have a direct biophysical interpretation. We discuss how this interesting feature of the H-H model allows one to reduce greatly the probabilistic degrees of freedom in uncertainty quantication analyses, saving CPU time in numerical simulations and opening possibilities for probabilistic generalisation of other deterministic models of great importance in physiology and mathematical neuroscience.

Keywords Neurodynamics Uncertainty Quantication Sparse grid quadrature

HodgkinHuxley model Neuronal Noise

1 Findings

To provide sensible mechanistic models of neuronal circuits in the brain and to facilitate principled interpretations of data from experiments in vitro and in vivo, biophysical models of neurodynamic systems must incorporate variability as a dening feature [16]. Neuronal variability cannot be accounted for by dynamic models

A. Torres Valderrama (B) J. Witteveen M. Navarro J. Blom

CWI, Amsterdam, Science Park 123, 1098 XG Amsterdam, The Netherlands e-mail: [email protected]

A. Torres Valderrama

Department of Neurology and Neurosurgery, Brain Center Rudolf Magnus, University Medical Center Utrecht, Heidelberglaan 100, 3584 CX Utrecht, The Netherlands

Uncertainty Propagation in Nerve Impulses Through the Action Potential Mechanism

Page 2 of 9 A. Torres Valderrama et al.

formulated as deterministic differential equations. However, such models can be extended to do so by incorporating probabilistic degrees of freedom. For this purpose, the versatile technique of Stochastic Collocation on Sparse Grids (SCSG) [79], from Uncertainty Quantication (UQ) theory [7, 1013], has been recently suggested as a computationally efcient strategy [14].

In this report, we deploy SCSG to investigate uncertainty propagation through the action potential mechanism in the H-H model [15], namely

x(t, ) =

gNaC m3h(v ENa) gKC n4(v EK) gLC (v EL) + I

0.01(10v)

exp[(10v)/10]1 (1 n) 0.125 exp[v80 ]n

0.1(25v)

exp[(25v)/10]1 (1 m) 4 exp[v18 ]m 0.07 exp[v20 ](1 h)

1

exp[(30v)/10]+1 h

, (1)

with x(t, ) = [v(t, ), m(t, ), n(t, ), h(t, )]T , where v(t, ) is the electrical po

tential across a neurons membrane, m(t, ), h(t, ) are the gating variables associated with the activation and inactivation of Na+ ion currents, respectively, n(t, ) is the gating variable associated with the activation of K+ ions current and where I represents current injected through the cell membrane.

The ODE system (1) exhibits a combination of excitability and highly nonlinear dynamics. We will show how these features can render the model output extremely sensitive to parameter uctuations under realistic physiological conditions, demonstrating that Uncertainty Quantication might be indispensable to provide sensible biophysical models of nerve impulses. The H-H model [15] includes a vector of 11 parameters, which in state space form reads

= [v0, m0, n0, h0, gNa, gK, gL, ENa, EK, EL, C], (2) consisting of four initial conditions v0, m0, n0, h0, three parameters describing the maximum conductances gNa, gK, gL, corresponding to Na+, K+ and leakage ion currents respectively, three parameters describing their equilibrium Nernst potentials

ENa, EK, EL, and the membrane capacitance C. These parameters are uncertain, for they must be determined experimentally and are therefore subject to errors and uctuations. Their nominal values measured in squid axon preparations are listed in Table 2 [15].

In the probabilistic generalisation of the H-H model, a Banach space describing such 11 probabilistic degrees of freedom is required for uncertainty analysis. Without appropriate discretisation and dimensionality reduction procedures, the computational cost of such analysis depends exponentially on the number of parameters, challenging the capacities of current hardware.

An efcient numerical strategy to deal with such a curse of dimensionality, is to nd a selection of points in probability space that yields a good approximation of the models response surface. This can be achieved by constructing a Sparse Grid in the multi-dimensional probability space for each time point following Smolyaks algorithm [8, 9, 1618].

Using Sparse Grid quadrature in conjunction with Sobol indices, we unravelled the parameters whose uctuation causes most of the observed variability in neuronal

Journal of Mathematical Neuroscience (2015) 5:3 Page 3 of 9

Fig. 1 a The top panel shows the probabilistic membrane potential, distributions are indicated by box-and-whiskers plots at each time instant. The middle panel shows the rst order Sobol indices, revealing the sources of the uctuations in the membrane potential for each uncertain parameter. The black line is the sum of all indices; deviation from one indicates variance due to parameter interactions. The bottom panel shows the variance difference w.r.t. the model with 11 uncertain parameters for four parsimonious models. Only two parameters (gK and ENa) are required to estimate a variability with less than 10 % error. b Average of the Sobol indices for all uncertain parameters during the time interval shown in a for the membrane potential and gating variables

responses. Such an analysis hinted at a parsimonious H-H model with two probabilistic degrees of freedom instead of 11, the use of which can achieve tremendous savings in CPU time as discussed below.

The uncertain dynamics of the membrane potential during neuronal discharge is shown in Fig. 1(a). The probability distribution at each time point was obtained by Monte Carlo sampling of a third level Sparse Grid piecewise linear interpolant [8, 19] of the response surface for the H-H model, using 1000 samples per time point from a uniform probability distribution and assuming 20 % of variability in the nominal values of the parameters listed in Table 2. The mean and deterministic solutions are also shown. Notice that they differ when the membrane depolarisation reaches its acme, which is often a signature of high nonlinearity.

In Fig. 1, the mean was obtained at the fourth level of a Sparse Grid constructed using the GaussPatterson (GP) quadrature rule, which requires 18591 function calls.

Page 4 of 9 A. Torres Valderrama et al.

These deterministic sequences are nested by construction, and for smooth integrands they can achieve the maximum degree of exactness from all nested rules [9, 18, 20].

In Fig. 1 neuronal variability is reported on the top panel by box-and-whisker plots, with one box per time point. On each box, the central mark is the median, the edges of the box are the 25th and 75th percentiles and the whiskers extend to the most extreme membrane potentials observed.

The middle panel of Fig. 1(a) shows the dynamics of the rst order Sobol indices [18, 2124], computed by Sparse Grid quadrature. The Sobol indices in Fig. 1(a) allow one to determine the relative contribution of each uncertain parameter to the variability of the model output observed in the top panel at each instant in time. The black line shows the sum of all indices; deviation from one indicates uctuations in the model output due to parameter interactions not accounted for by the rst order indices.

The variance and Sobol indices were obtained at the third level of a GP Sparse Grid in 22 dimensions, requiring 17249 function calls. Notice that only 21 dimensions are needed to obtain the rst order Sobol indices via the conditional variance. However, adding one extra dimension allows one to estimate the total variance with the same Sparse Grid, which often increases accuracy [18].

The RMS errors for the mean, variance and Sobol indices shown in Fig. 1 are summarised in Table 1. The most accurate numerical solution is taken as a reference for error analysis, thus the error estimates listed in the table are conservative, since they approximate the numerical error at the last but one quadrature level. Table 1 shows that the numerical errors are sufciently small, validating the choice of quadrature level for further analysis.

In Fig. 1(b) all uncertainty sources are represented in a colour map, showing the average value of all Sobol indices during the time interval shown in Fig. 1(a). All outputs of the model, namely the membrane potential and the three gating variables are shown.

Interestingly, this colour map reveals that only a small subset of parameters causes most of the uncertainty in the model outputs. Therefore, all other parameters can be xed to their nominal values preserving the probabilistic interpretation of the model. This is illustrated in the bottom panel of Fig. 1(a), which shows the difference in standard deviation of four effective models with respect to that of the model with the 11 uncertain parameters listed in Table 2. Notice that only two parameters, gK and

ENa, are required to estimate a variability with less than 10 % difference w.r.t. the original model for all time points.

The singular value decomposition (SVD) of the (11 4) data matrix M displayed

in Fig. 1(b) as a colour map, namely M = USV T , provides a quantitative ranking of

the average Sobol indices. The rst four columns of U after normalising the patterns so that each column has maximum entry 1 are shown in Table 2. Notice that the rst mode (rst column of U) contains more than 90 % of the energy, given by the square of the singular values. Notice also that in this maximum energy mode, two probabilistic degrees of freedom are dominant, namely gK and ENa. The inuence of all the other parameters in the variability of the model output is small. Therefore, they can be xed to their nominal values. This greatly reduces the dimensionality of the model retaining most of its probabilistic features.

Table 1 RMS error of the Sobol indices show in Fig. 1(a) for all model outputs

E Sn0 Sm0 Sh0 v0 SgK SgNa SgL SEK SENa SEL SC

n(t, ) 0.000163 0.000350 0.016969 0.016969 0.016710 0.016970 0.012563 0.009419 0.016891 0.015160 0.011181 0.016949 0.015736

m(t, ) 0.001494 0.004814 0.090669 0.090668 0.089543 0.090673 0.097756 0.057230 0.090059 0.078554 0.062516 0.090535 0.082084

h(t, ) 0.000112 0.000271 0.046205 0.046205 0.045488 0.046202 0.035574 0.027616 0.045989 0.042473 0.017223 0.046155 0.042037

v(t, ) 0.004669 0.052106 0.087004 0.087002 0.084896 0.086993 0.079195 0.055766 0.086538 0.073698 0.060872 0.086917 0.078468

JournalofMathematicalNeuroscience(2015)5:3Page5of9

Page 6 of 9 A. Torres Valderrama et al.

Table 2 Normalised ranking of the average Sobol indices displayed in the colour map in Fig. 1(b) using SVD

Mode 1 Mode 2 Mode 3 Mode 4

Mode energy: 0.944690 0.042433 0.011498 0.001379

Nominal values

n0 = 0.0003; 0.297 0.001064 0.000372 0.012061 0.034643 m0 = 0.0011; 0.00616 0.001064 0.000372 0.012060 0.034648 h0 = 0.9998; 0.12 0.140569 0.436071 0.555783 1.000000 v0 = 10 mV; 0.001 0.001206 0.000072 0.013534 0.059476 gK = 36 S/cm

2 1.000000 0.627724 0.197631 0.384254 gNa = 120 S/cm

2 0.172271 0.300390 0.305118 0.084833 gL = 0.3 S/cm

2 0.001741 0.001124 0.016312 0.091209 EK = 12 mV 0.232445 0.225011 1.000000 0.978672

ENa = 115 mV 0.532895 1.000000 0.190240 0.759104 EL = 10.613 mV 0.001293 0.000515 0.013721 0.048586

C = 1 F/cm

2 0.099635 0.342083 0.410624 0.931327

The symbol indicates the initial values used in Fig. 2.

This feature of the H-H model is important, since it shows that probabilistic neurodynamic simulations can be performed parsimoniously, in an effective probability subspace of dimension far smaller than that of the original model. This greatly reduces the computational cost and facilitates the analysis of neuronal systems beyond individual cells.

As a non-trivial illustrative application of our parsimonious model Fig. 2 shows the impact of parameter uctuations on two relevant neurocomputational properties, namely, neural response to current input and membrane refractoriness.

In Fig. 2(a) two 1 ms current pulses are applied to a deterministic H-H neuron whose membrane potential is initially at the resting potential. The rst pulse of 8.5 A/cm2 increases the membrane potential, but it does not elicit neuronal discharge. A second stronger pulse of 20 A/cm2 applied later, immediately triggers a spike. The probabilistic counterpart of this experiment is shown in (c). The large error bars in response to the rst pulse show that a small current pulse can indeed trigger action potentials contrary to deterministic predictions.

In panels (b) and (d) the strong current pulse of 20 A/cm2 is applied rst, immediately eliciting a spike. In the deterministic model shown in panel (b), a second small pulse of 10.5 A/cm2 fails to ignite an action potential due to membrane refractoriness. Under parameter uncertainty, panel (d) shows large variability in the neuronal response, indicating a second action potential, which is impossible on deterministic grounds.

Synaptic input responses and membrane refractoriness are neurocomputational properties of realistic nerve cells, which directly inuence collective neuronal behaviour as well as the efciency of neural codes. Thus, the examples above illustrate how incorporating variability in neurodynamic models might be crucial to our understanding of the nervous system and the behaving brain.

Journal of Mathematical Neuroscience (2015) 5:3 Page 7 of 9

Fig. 2 In panel a two current pulses are applied to a deterministic neuron with the nominal values for the parameters listed in Table 2. The rst pulse is not strong enough to elicit a spike, the second stronger pulse immediately triggers a spike. Panel c shows the probabilistic counterpart of this experiment assuming 20 % of variability in the nominal parameters gK and ENa. Large error bars show that the rst small current pulse can trigger action potentials in this instance. In panels b and d the stronger current pulse is applied rst, immediately eliciting a spike. In the deterministic model shown in panel b, the second small pulse fails to trigger a spike when applied during the refractory period. In the probabilistic model d, large variability in the output indicates a second action potential, which would not be expected from deterministic predictions

Probabilistic generalisation of other deterministic models of great importance in mathematical neuroscience might be possible using the methods and results in this report, such as those describing neuronal interactions and signal transmission in active neural media, at a feasible computational cost.

Competing Interests

The authors declare that they have no competing interests.

Authors Contributions

All authors contributed equally to the writing of this paper. All authors read and approved the nal manuscript.

Page 8 of 9 A. Torres Valderrama et al.

Acknowledgements We thank Prof. Stan Gielen for useful discussions and suggestions. This work was supported by European Commissions 7th Framework Program BioPreDyn grant number 289434. ATV gratefully acknowledges the nancial support of the Exact Sciences division of the Netherlands Organisation for Scientic Research NWO grant number: veni 639.071.905.

References

1. Fellous JM, Rudolph M, Destexhe A, Sejnowski TJ: Synaptic background noise controls the input/output characteristics of single cells in an in vitro model of in vivo activity. Neuroscience 2003, 122(3):811-829.

2. Destexhe A, Rudolph M, Fellous J-M, Sejnowski TJ: Fluctuating synaptic conductances recreate in vivo-like activity in neocortical neurons. Neuroscience 2001, 107(1):13-24. doi:10.1016/ S0306-4522(01)00344-X.

3. Par D, Shink E, Gaudreau H, Destexhe A, Lang EJ: Impact of spontaneous synaptic activity on the resting properties of cat neocortical pyramidal neurons in vivo. J Neurophysiol 1998, 79:1450-1460.

4. Destexhe A, Par D: Impact of network activity on the integrative properties of neocortical pyramidal neurons in vivo. J Neurophysiol 1999, 81(4):1531-1547.

5. H N, Destexhe A: Synaptic background activity enhances the responsiveness of neocortical pyramidal neurons. J Neurophysiol 2000, 84:1488-1496.

6. Destexhe A, Rudolph-Lilith M: Neuronal Noise. Dordrecht: Springer; 2012. [Springer Series in Computational Neuroscience.]

7. Nobile F, Tempone R, Webster CG: A sparse grid stochastic collocation method for partial differential equations with random input data. SIAM J Numer Anal 2008, 46(5):2309-2345.

8. Bungartz H-J, Griebel M: Sparse grids. Acta Numer 2004, 13:147-269. doi:10.1017/ S0962492904000182.

9. Gerstner T, Griebel M: Numerical integration using sparse grids. Numer Algorithms 1998, 18:209-232.

10. Le Matre OP, Knio OM: Spectral Methods for Uncertainty Quantication: With Applications to Computational Fluid Dynamics. Dordrecht: Springer; 2010. [Scientic Computation.]

11. Grigoriu M: Stochastic Systems: Uncertainty Quantication and Propagation. New York: Springer;

2012. [Springer Series in Reliability Engineering.]

12. Xiu D: Numerical Methods for Stochastic Computations: A Spectral Method Approach. Princeton: Princeton University Press; 2010.

13. Ghanem RG, Spanos PD: Stochastic Finite Elements: A Spectral Approach. Mineola: Dover; 2003. [Civil, Mechanical and Other Engineering Series.]

14. Laing C, Zou Y, Smith B, Kevrekidis I: Managing heterogeneity in the study of neural oscillator dynamics. J Math Neurosci 2012, 2(1):5. doi:10.1186/2190-8567-2-5.

15. Hodgkin AL, Huxley AF: A quantitative description of membrane current and its application to conduction and excitation in nerve. J Physiol 1952, 117(4):500-544.

16. Novak E, Ritter K: High dimensional integration of smooth functions over cubes. Numer Math 1996, 75(1):79-97. doi:10.1007/s002110050231.

17. Smolyak SA: Quadrature and interpolation formulas for tensor products of certain classes of functions. Dokl Akad Nauk SSSR 1963, 4:240-243.

18. Torres Valderrama A, Blom J: Uncertainty propagation in neuronal dynamical systems. Technical report CWI LS-1304; 2012.

19. Klimke A, Wohlmuth B: Algorithm 847: Spinterp: piecewise multilinear hierarchical sparse grid interpolation in MATLAB. ACM Trans Math Softw 2005, 31(4):561-579. doi:10.1145/ 1114268.1114275.

20. Patterson TNL: An algorithm for generating interpolatory quadrature rules of the highest degree of precision with preassigned nodes for general weight functions. ACM Trans Math Softw 1989, 15(2):123-136. doi:10.1145/63522.63523.

21. Homma T, Saltelli A: Importance measures in global sensitivity analysis of nonlinear models. Reliab Eng Syst Saf 1996, 52:1-17.

22. Ishigami T, Homma T: An importance quantication technique in uncertainty analysis for computer models. In International Symposium on Uncertainity Modelling and Analysis (ISUMA90). December 36, University of Maryland; 1990.

Journal of Mathematical Neuroscience (2015) 5:3 Page 9 of 9

23. Sobol IM: Sensitivity estimates for nonlinear mathematical models. Math Models 1990, 2:112-118 (in Russian).

24. Saltelli A: Making best use of model evaluations to compute sensitivity indices. Comput Phys Commun 2002, 145(2):280-297. doi:10.1016/S0010-4655(02)00280-1.