Reduction of stochastic conductance-based neuron models with time-scales separation

Abstract

We introduce a method for systematically reducing the dimension of biophysically realistic neuron models with stochastic ion channels exploiting time-scales separation. Based on a combination of singular perturbation methods for kinetic Markov schemes with some recent mathematical developments of the averaging method, the techniques are general and applicable to a large class of models. As an example, we derive and analyze reductions of different stochastic versions of the Hodgkin Huxley (HH) model, leading to distinct reduced models. The bifurcation analysis of one of the reduced models with the number of channels as a parameter provides new insights into some features of noisy discharge patterns, such as the bimodality of interspike intervals distribution. Our analysis of the stochastic HH model shows that, besides being a method to reduce the number of variables of neuronal models, our reduction scheme is a powerful method for gaining understanding on the impact of fluctuations due to finite size effects on the dynamics of slow fast systems. Our analysis of the reduced model reveals that decreasing the number of sodium channels in the HH model leads to a transition in the dynamics reminiscent of the Hopf bifurcation and that this transition accounts for changes in characteristics of the spike train generated by the model. Finally, we also examine the impact of these results on neuronal coding, notably, reliability of discharge times and spike latency, showing that reducing the number of channels can enhance discharge time reliability in response to weak inputs and that this phenomenon can be accounted for through the analysis of the reduced model.

Appendices

Appendix A: Auxiliary functions

The transition rates for the HH model are given by:

$$ \begin{array}{rll} \alpha_n(V)&=&(0.1-0.01V)/[\exp(1-0.1V)-1],\\ \beta_n(V)&=&0.125\exp(-V/80),\\ \alpha_m(V)&=&(2.5-0.1V)/[\exp(2.5-0.1V)-1],\\ \beta_m(V)&=&4\exp(-V/18),\\ \alpha_h(V)&=&0.07\exp(V/20),\\ \beta_h(V)&=&1/[\exp(3-0.1V)+1]. \end{array} $$

Accordingly, τ _x (V) and x _∞(V), for x = m, n, h, are given by: τ _x (V) = 1/(α _x (V) + β _x (V)) and x _∞(V) = α _x (V)/(α _x (V) + β _x (V)).

Appendix B: Diffusion term in the two-state model

The purpose of this appendix is to give more details about the computation of the diffusion term in the two-state model. The computation is more complicated than in the multi-state model because of the non-linearity m ³. Indeed, in the multi-state case, the vector field F _MS is linear with respect to u _Na and one only needs to compute the corrective diffusion for a single channel, and then divide it by $\sqrt{N_{Na}}$, since the average of N _Na independent Brownian motions is equal in law to a single Brownian motion divided by $\sqrt{N_{Na}}$. In the non-linear case (two-state model), one needs to work out the computation directly on the process describing the empirical measure:

• 1. First, one has to compute the law at time t of the empirical measure, with V fixed, that is the probability $P_N^V(j,j_0,t)$ of having j open channels at time t starting from a population with j ₀ open channels.

  Denote N = N _m . Starting from a proportion of open m gates $u_N(0)=\frac{j_0}{N}$ at t = 0, with 0 ≤ j ₀ ≤ N, one can show that the proportion of open gates u _N (t) at time t follows a bi-binomial distribution:

  $$\begin{array}{rll} && P_N^V(j,j_0,t)\\ &&:=\mathbf{P}\left[u_N(t)=\frac{j}{N}|u_N(0)=\frac{j_0}{N}\right]\\ &&\displaystyle{\sum_{x=\max(0,j-N+j_0)}^{\min(j,j_0)}}\mu_{x,j_0}\left(p_0(t)\right)\mu_{j-x,N-j_0}\left(p_1(t)\right) \end{array}$$

  with $\mu_{i,j}(p)=C_{j}^i p^{i}(1-p)^{j-i}$ and where p ₀(t) and p ₁(t) are the solution of

  $$ \dot{y}=(1-y)\alpha_m(V)-y\beta_m(V) $$

  (22)

  with respective initial conditions p ₀(0) = 0 and p ₁(0) = 1. Defining

  $$ \tau_m:= \alpha_m(V)+\beta_m(V) \mbox{ and } m_{\infty}(V):=\frac{\alpha(V)}{\tau_y(V)} $$

  (23)

  the variation of constant gives: p ₀(t) = m _∞(V)$\left(1-e^ {-t/\tau_m(V)}\right)$ and $p_1(t)=e^ {-t/\tau_m(V)}+p_0(t)$. If j ₀ = 0, then one retrieves the classic binomial distribution with parameters (N,p ₀(t)): $P_N^V(j,0,t)=C_N^j p_0(t)^j (1-p_0(t))^{N-j}$. When t →  ∞, one checks that the quasi-stationary distribution

  $$ \rho_N^V(j)=\displaystyle{\lim\limits_{t\to \infty}} P_N^V(j,j_0,t) $$

  does not depend on j ₀ and is given by:

  $$ \rho_N^V(j)= \mu_{j,N}(m_{\infty}) = C_N^j m_{\infty}(V)^j (1-m_{\infty}(V))^{N-j} $$

• 2. Then, integrating with respect to time the difference between $P_N^V(j,i,t)-\rho_N^V(j)$, we can express the quantity R(i,j) required to compute the diffusion term. Introducing K = (1 − m _∞(V))/ (m _∞(V)), and making the change of variable $z=e^{-t/\tau_m(V)}$, the computation boils down to integrals of the form:

  $$ \int_0^1 (1 - z)^a (1 + Kz)^b (1 + z/K)^c dz $$

  Using a formal computation software, R(i, j) can be expressed as a sum involving Appell F1 functions, defined by

  $$ F_1(a,b_1,b_2,c;x,y) = \sum\limits_{m,n=0}^\infty \frac{(a)_{m+n} (b_1)_m (b_2)_n} {(c)_{m+n} \,m! \,n!} \,x^m y^n $$

  with (q) _n  = q (q + 1) ⋯ (q + n − 1). Then we can write:

  $$\begin{array}{rll} R(i,j)&=&\tau_m\sum\limits_{x=\max(0,j-N+i)}^{\min(i,j)} C_i^x C_{N-i}^{j-x}Y_{\infty}^{(N)}(x,i,j)\\ &&\times \left(H_K^{(N)}(x,i,j)F_1^{(N)}(x,i,j)-1\right) \end{array}$$

  with:

  $$\begin{array}{rll} Y_{\infty}^{(N)}(x,i,j)&:=& m_{\infty}^x m_{\infty}^{j-x} (1-m_{\infty})^{i-x}\\ &&\times \ (1-m_{\infty})^{N-i-j+x}\\ H_{K}^{(N)}(x,i,j)&:=&\frac{K^{x-i}(1+K)^{i+j-2x}}{1+2x-N-i-j}\\ F_1^{(N)}(x,i,j)&:=& F_1\left(w,x-i,x-j,w+1,\right.\\ &&\qquad\left.\frac{1}{1+K},\frac{K}{1+K}\right)\\ w&:=&1+2x+N-i-j \end{array}$$

• 3. The complete expression for the variance is then given by plugging the above expression for R(i, j) into formula (20).