Stochastic modeling of excitable dynamics: improved Langevin model for mesoscopic channel noise
aa r X i v : . [ q - b i o . N C ] M a y Stochastic modeling of excitable dynamics:improved Langevin model formesoscopic channel noise
Igor Goychuk
Institute for Physics and Astronomy, University of PotsdamKarl-Liebknecht-Str. 24/25, 14476 Potsdam-Golm, Germany
Abstract.
Influence of mesoscopic channel noise on excitable dynamicsof living cells became a hot subject within the last decade, and the tradi-tional biophysical models of neuronal dynamics such as Hodgkin-Huxleymodel have been generalized to incorporate such effects. There still ex-ists but a controversy on how to do it in a proper and computationallyefficient way. Here we introduce an improved Langevin description ofstochastic Hodgkin-Huxley dynamics with natural boundary conditionsfor gating variables. It consistently describes the channel noise variancein a good agreement with discrete state model. Moreover, we show bycomparison with our improved Langevin model that two earlier Langevinmodels by Fox and Lu also work excellently starting from several hun-dreds of ion channels upon imposing numerically reflecting boundaryconditions for gating variables.
Keywords:
Excitable dynamics, ion channels, mesoscopic noise, Langevindescription
Hodgkin-Huxley (HH) model of neuronal excitability [1] provides a milestonefor biophysical understanding of information processing in living systems [2] interms of electrical spikes mediated by ionic currents through voltage-dependentmembrane pores made by ion channel proteins. One considers the cell membraneas an insulator with specific electrical capacitance C m per unit of area, whichis perforated by ionic channels providing generally voltage-dependent parallelionic pathways with specific conductances G i per unit of area for various sortsof ion channels. This yields the following equation for transmembrane electricalpotential difference VC m dVdt + G K ( n )( V − E K ) + G Na ( m, h )( V − E Na ) + G L ( V − E L ) = I ext . (1)Here, three ionic currents are taken into account, sodium Na, potassium K andunspecific leakage current (mainly due to chloride ions). This is nothing else Stochastic modeling of excitable dynamics the Kirchhoff current law, which takes into account the ionic and capacitancecurrents, as well as an external current I ext which can stimulate electrical exci-tations. This equation reflects assumption on Ohmic conductance of completelyopen ion channels with E i being the reversal or Nernst potentials. They emergedue to the difference of ionic concentrations inside and outside of the excitablecell, which are kept approximately constant by the work of ionic pumps, which isnot considered explicitely. Nonlinearity comes from the open-shut gating dynam-ics of sodium and potassium channels. The corresponding specific conductances G K ( n ) = g maxK n ( V, t ) ,G Na ( m, h ) = g maxNa m ( V, t ) h ( V, t ) , (2)depend on three voltage-dependent gating variables, n , m , and h , where n ( t ) isthe probability of one gate of potassium channel to be open (more precisely thefraction of open gates), m corresponds to one activation gate of sodium channel,and h is the fraction of closed sodium inactivation gates. One assumes fourindependent identical gates for potassium channel, hence its opening probabilityis n , as well as three activation and one inactivation gate for the sodium channel.Hence, m h is the fraction of open sodium channels. The maximal conductances g maxK and g maxNa can be expressed via the unitary conductances g i, of single ionchannels as g max i = g i, ρ i , where ρ i is the membrane density of the ion channelsof sort i . The gating dynamics is in turn described by the relaxation kinetics ddt x = α x ( V ) (1 − x ) − β x ( V ) x, x = m, h, n , (3)with voltage-dependent rates α m ( V ) = 0 . V + 40)1 − exp[ − ( V + 40) / , β m ( V ) = 4 exp[ − ( V + 65) / , (4) α h ( V ) = 0 .
07 exp[ − ( V + 65) / , β h ( V ) = { − ( V + 35) / } − , (5) α n ( V ) = 0 .
01 ( V + 55)1 − exp[ − ( V + 55) / , β n ( V ) = 0 .
125 exp[ − ( V + 65) / . (6)Here the voltage is measured in millivolts and rates in inverse milliseconds. Otherclassical parameters of HH model suitable to describe excitable dynamics of squidgiant axon are: C m = 1 µ F / cm , E Na = 50 mV, E K = −
77 mV, E L = − . G L = 0 . / cm , g maxK = 36 mS / cm , g maxNa = 120 mS / cm .The set of four coupled nonlinear differential equations defined by (1)- (6)presents a milestone in biophysics and neuroscience because of its very clearand insightful physical background. In the same spirit, one can build up variousother conductance-based biophysical models starting from the pertinent molec-ular background and following to the bottom-up approach. However, it assumesmacroscopically large numbers of ion channels in neglecting completely the meso-scopic channel noise effects. The number of ion channels in any cell is, however,finite, and the corresponding channel noise can be substantial [2]. Especially, oneconfronts with this question by considering the spatial spike propagation amongapproximately piece-wise isopotential membrane clusters of ion channels [3]. mproved Langevin model for mesoscopic channel noise 3 How to take stochastic dynamics of ion channels within the physical frame-work of HH model into account is subject of numerous studies [2,4,5]. The mostrigorous way is to consider the variable population of open ion channels as abirth-and-death process [6]. Consider for simplicity a population of N indepen-dent two-state ion channels (one gate only) with opening rate α and closing rate β (constant under voltage clamp). Each ion channel fluctuates dichotomouslybetween the closed state with zero conductance and the open state having uni-tary conductance g . For such two-state Markovian channels, the stationaryprobability of opening is p = α/ ( α + β ) and the averaged conductance is h g ( t ) i = p g . The number of open channels n is binomially distributed withprobability P st N ( n ) = p n (1 − p ) N − n N ! / ( n !( N − n )!), average h n i = p N , andvariance h ( n − h n i ) i = N p (1 − p ) = N αβ/ ( α + β ) . For sufficiently large N ≥ ≤ x ( t ) = n ( t ) /N ≤
1, withsmoothened binomial probability density p st N ( x ) = N p xN (1 − p ) N (1 − x ) N ! / ( Γ (1 + xN ) Γ (1 + (1 − x ) N )). Use of approximate Stirling formula n ! ≈ ( n/e ) n yields p st N ( x ) ≈ C N ( α, β ) (cid:16) αx (cid:17) Nx (cid:18) β − x (cid:19) N (1 − x ) , (7)where C N ( α, β ) is a normalization constant. We are looking for the best dif-fusional (continuous) approximation for discrete state birth-and-death processdefined by the master equation˙ P N ( n ) = F ( n − P N ( n −
1) + B ( n + 1) P N ( n + 1) − [ F ( n ) + B ( n )] P N ( n )(8)for 1 ≤ n ≤ N −
1, with forward rate F ( n ) = α ( N − n ) and backward rate B ( n ) = βn , complemented by the boundary conditions˙ P N (0) = B (1) P N (1) − F (0) P N (0) , (9)˙ P N ( N ) = F ( N − P N ( N − − B ( N ) P N ( N ) . (10) Astandard way to obtain diffusional approximation for p N ( x ) := P N ( xN ) /∆x ( ∆x = 1 /N ) with rates f ( x ) := F ( xN ) ∆x , b ( x ) := B ( xN ) ∆x is to do theKramers-Moyal expansion [6,7], like p N ( x + ∆x ) ≈ p N ( x ) + ( ∂p N ( x ) /∂x ) ∆x +( ∂ p N ( x ) /∂x )( ∆x ) / f ( x + ∆x ) ≈ f ( x )+( df ( x ) /dx ) ∆x +( d f ( x ) /dx )( ∆x ) / ∂∂t p ( x, t ) = − ∂∂x [ f ( x ) − b ( x )] p ( x, t ) + ∂ ∂x D KM ( x ) p ( x, t ) (11)with diffusion coefficient D KM ( x ) = [ f ( x ) + b ( x )] / (2 N ). This Fokker-Planckequation corresponds to the Langevin equation˙ x = f ( x ) − b ( x ) + p D KM ( x ) ξ ( t ) , (12) Stochastic modeling of excitable dynamics where ξ ( t ) is white Gaussian noise of unit intensity, h ξ ( t ) ξ ( t ′ ) i = δ ( t − t ′ ), inpre-point, or Ito interpretation [8]. This equation is quite general for any one-dimensional birth-and-death process within this standard diffusional approxima-tion. For the considered population of ion channels,˙ x = α (1 − x ) − βx + p [ α (1 − x ) + βx ] /N ξ ( t ) . (13)This is stochastic equation for a gating variable in the stochastic generalizationof Hodgkin-Huxley equations by Fox and Lu [5]. It replaces Eq. (3) with corre-sponding voltage-dependent α x ( V ), β x ( V ), and N = N Na = ρ Na S for m and h ,or N = N K = ρ K S for the variable n . S is the area of membrane patch, and ρ Na = 60 µm − , ρ K = 18 µm − within HH model [4]. Clearly, in the limit N → ∞ the channel noise vanishes, restoring the deterministic HH model. We name thismodel the second model by Fox and Lu (Fox-Lu 2) in application to stochasticHH dynamics. Linear noise approximation.
The further approximation (Fox-Lu 1 withinstochastic HH model) is obtained by D KM ( x ) → D KM ( x eq ) = const , where x eq isequilibrium point of deterministic dynamics, f ( x eq ) = b ( x eq ). It corresponds tothe so-called 1 /Ω expansion with linear additive noise approximation advocatedby van Kampen [6]. Then, with x eq = p = α/ ( α + β ) Eq. (13) reduces to˙ x = α (1 − x ) − βx + p αβ/ [ N ( α + β )] ξ ( t ) . (14) Diffusional approximation with natural boundaries.
The both diffusionalapproximations are not quite satisfactory because they do not guarantee theboundary conditions in a natural way. As a result, for a sufficiently small open-ing probability p ≪
1, and not sufficiently large number of channels the negativevalues, x <
0, become possible with appreciable probabilities p ( x, t ) >
0. Like-wise, the larger than one values, x >
1, are also possible when the openingprobability p is close to one. However, this deficiency can easily be correctednumerically by imposing reflecting boundary conditions at x = 0 and x = 1 instochastic simulations. With this correction, Langevin approximation of stochas-tic HH dynamics is widely used [9,10,11,3]. However, it is not quite clear if thisprocedure indeed delivers the correct results [12]. To clarify the issue, we con-sider a different diffusional approximation with natural reflecting boundarieswhich naturally bound stochastic dynamics to the interval 0 ≤ x ≤ P st N ( n ) = exp[ − Φ ( n )] P st N (0) in terms of a pseudo-potential Φ ( n ) = − P ni =1 ln [ F ( i − /B ( i )] [6]. Hence, in the continuous limit, p st N ( x ) ∝ exp[ − N φ ( x )],with pseudo-potential φ ( x ) = − R x ln [ f ( x ′ ) /b ( x ′ )] dx ′ = ln(1 − x ) − x ln( α (1 − x ) / ( xβ )). This indeed yields the probability density (7). The correspondingFokker-Planck equation must read ∂∂t p ( x, t ) = ∂∂x (cid:18) D ( x ) e − Nφ ( x ) ∂∂x e Nφ ( x ) p ( x, t ) (cid:19) (15) mproved Langevin model for mesoscopic channel noise 5 = ∂∂x N D ( x ) φ ′ x ( x ) p ( x, t ) + ∂∂x D ( x ) ∂∂x p ( x, t ) (16)with N D ( x ) φ ′ x ( x ) = b ( x ) − f ( x ) (17)in order to be also consistent with the deterministic limit N → ∞ . The lastequation fixes the diffusion coefficient as D ( x ) = 1 N f ( x ) − b ( x )ln[ f ( x ) /b ( x )] . (18)The Langevin equation which corresponds to this best diffusional approximationof the birth-and-death processes [7,13] reads˙ x = f ( x ) − b ( x ) + p D ( x ) ξ ( t ) , (19)in the post-point, or Klimontovich-H¨anggi interpretation [7]. In the standardIto interpretation suitable for integration with stochastic Euler algorithm [8] thecorresponding Langevin equation becomes˙ x = f ( x ) − b ( x ) + D ′ x ( x ) + p D ( x ) ξ ( t ) (20)with spurious drift D ′ x ( x ). In application to stochastic dynamics of one gatingvariable it reads ˙ x = α (1 − x ) − βx + D ′ x ( x ) + p D ( x ) ξ ( t ) . (21)with D ( x ) = 1 N α (1 − x ) − βx ln[ α (1 − x ) / ( βx )] . (22)Replacing with such equations the stochastic equations for gating variables inthe standard Langevin variant of stochastic Hodgkin-Huxley equations we ob-tain the improved Langevin description of mesoscopic channel noise, with naturalboundaries because D (0) = D (1) = 0, i.e. the channel noise (and the probabilityflux) vanishes exactly at the reflecting boundaries, in the theory. Nevertheless,in numerical algorithm one must yet additionally secure such boundaries for any finite integration time step δt . Notice also that near the equilibrium point with | f ( x ) − b ( x ) | ≪ f ( x ) + b ( x ), D ( x ) ≈ D KM ( x ), and the standard diffusional ap-proximation is almost restored, almost, if to neglect the spurious drift correction D ′ x ( x ), which still remains within the Ito interpretation.We test the best diffusional approximation for a gating variable against theearlier Langevin descriptions with reflecting boundary conditions implementednumerically. For this we use stochastic Euler algorithm with time step δt = 0 . N and the simulation software XPPAUT [14]. The resultsare shown for α = 1 and β = 9 with p = 0 . N = 100 (a) and Stochastic modeling of excitable dynamics N = 10 (b). As a big surprise, the simplest linear noise approximation actuallyseems to work best, if only to implement reflecting boundary conditions. For N = 100, it reproduces well the still somewhat skewed binomial distributionwith the exact mean h x i = 0 . h ∆x i / = 0 .
03. Even for N = 10, it gives the mean closer to the correct value of 0 . h ∆x i / ≈ .
095 larger than within two other approximations. For a sufficientlylarge N = 1000 (not shown), all three diffusional approximations give practicallyidentical results, within the statistical errors of simulations. Surprisingly, allthree work reasonably well even for N = 10! However, such a performance is apriori not guaranteed for stochastic nonlinear dynamics with voltage-dependent α ( V ( t )) and β ( V ( t )). In fact, for a multistable dynamics the best diffusionalapproximation is generically expected [13] to operate essentially better. x p s t N ( x ) Kramers-Moyal:
Stationary distributions of gating variable x for two ensembles of ion channelswith α = 1 and β = 9, (a) N = 100 and (b) N = 10. Numerics are compared withbinomial distribution (a) and distribution (7) for the best diffusional approximation We compare three different Langevin descriptions of stochastic HH dynamicsin Fig. 2, for two different membrane patches. Here, the interspike interval dis-tributions are presented, together with the corresponding mean, h τ i , standarddeviation, h [ τ − h τ i ] i / , and the relative standard variation, or the coefficientof variation, C V = h [ τ − h τ i ] i / / h τ i , which measures the spike coherence. For S = 10 µm , all three approximations agree well. However, for S = 1 µm thediscrepancies become apparent, and we prefer our improved Langevin descriptionon general theoretical grounds.The coefficient of variation C V , calculated within our Langevin variant ofstochastic HH model, is plotted vs. the patch size S in Fig. 3. It displays a typicalcoherence resonance [15] behavior revealed earlier within stochastic HH modelsin [9,16] as a system-size coherence resonance. There exists an optimal patchsize (optimal number of ion channels) with most coherent stochastic dynamicsdue to internal mesoscopic noise. mproved Langevin model for mesoscopic channel noise 7 In this paper, we presented the best diffusional Langevin approximation forexcitable cell dynamics within stochastic Hodgkin-Huxley model, with naturalboundary conditions for the channel noise implemented. It has clear theoreticaladvantages over the standard diffusional approximation in the case of transitionsinduced by mesoscopic noise as discussed for bistable birth-and-death processeslong time ago [13]. However, within stochastic HH model for a sufficiently largenumber of ion channels, the standard diffusional approximations were shown towork also very good. Hence, this work confirmes the validity of the previouswork done within the Langevin approximations of stochastic HH dynamics, fora sufficently large number of channels. This does not mean, however, that forother excitable models the situation will not be changed. Generally, the improvedLangevin description should operate better. Other stochastic models of excitabledynamics, e.g. stochastic Morris-Lecar model can easily be improved accordingly.This task, as well as comparison with discrete state stochastic models for channelnoise, is left for a future investigation.
Acknowledgments.
Support by the DFG (German Research Foundation),Grant GO 2052/1-2, is gratefully acknowledged. interspike interval, msec I S I d i s t r i bu ti on Fox-Lu 2: < τ > =39.29, < ∆τ > =24.67, CV=0.63best diffusional: < τ > =38.97, < ∆τ > =23.84, CV=0.61Fox-Lu 1: < τ > =39.98, < ∆τ > =24.30, CV=0.61 a) I S I d i s t r i bu ti on Fox-Lu 2: < τ > =16.71, < ∆τ > =10.83, CV=0.65best diffusional: < τ > =16.67, < ∆τ > =9.76, CV=0.59Fox-Lu 1: < τ > =20.02, < ∆τ > =10.55, CV=0.53 b) Fig. 2.
Interspike time interval distribution for self-excitable dynamics, I ext = 0, due tothe channels noise for two membrane patches: (a) S = 10 µm ( N Na = 600, N K = 180),and (b) S = 1 µm ( N Na = 60, N K = 18). References
1. Hodgkin, A. L., Huxley, A. F.: A quantitative description of membrane currentand its application to conduction and excitation in nerve. J. Physiol. (London)117, 500–544 (1952).2. Koch, Ch.: Biophysics of Computation, Information Processing in Single Neurons.Oxford University Press, Oxford (1999). Stochastic modeling of excitable dynamics S, µ m C V Fig. 3.
Coefficient of variation versus the membrane patch size within our variant ofstochastic HH model. Self-excitable dynamics, I extext