4.2 Conductance-based refractory density (CBRD) approach for multiple state variable neurons

Because of their multiple state variables, evaluating noisy neural populations with non-instantaneous synapses, multiple compartments and/or active HH channels with a KFP-based method, as described above, cannot be evaluated analytically, and thus a KFP (or Monte-Carlo) approach requires significant computational resources [11, 69]. To address this problem, we employ a CBRD approach that can evaluate a broad class of models in a one dimensional phase space parameterized by the refractory time, or time since the last spike, t* [17, 25, 70]. In the present case, this approach is used to evaluate all the two population models, including the HH models A and B, and the LIF models C and D.

In conjunction with a given neuron model, the CBRD approach considers a hazard function H which defines the probability density of firing, and describes the dynamics of a population as a whole by the probability density of neurons distributed according to their refractory times, ρ(t, t*). In particular, the CBRD model substitutes the phenomenologically derived hazard function developed in the original refractory density approach [68] by its rigorously derived expression [17, 70]. The population model for the neuron model is expressed by a set of transport equations, one for each state variable, in partial derivatives with t* and time t as independent variables, and thus the number of equations is directly proportional to the complexity of the model.

In order to obtain an expression for H in a probabilistic framework, e.g. in the presence of noise, the CBRD approach requires that the state variables be well described at any given t* by a linear model around their mean value. This constraint holds for most of state variables of the HH model given their voltage-dependence and a typical dispersion of the membrane voltage of several millivolts. For the LIF model, this constraint is automatically satisfied.

In contrast, the strong nonlinear properties of sodium channels of the HH model near threshold do not meet this criterion. To address this problem, and as mentioned in Section 2.1.3, we use a threshold HH cell model without sodium channels, with spikes defined by an explicit dynamic threshold and renewal function. This model is based on three important assumptions: first, that the probability of firing can be well predicted by the dynamics of a neuron without an explicit sodium current, that the influence of strongly non-linear channel dynamics, namely the sodium current, is only significant during the spike, and third, that there is a fixed (i.e. renewal process) impact of each spike on the state variables. The last assumption allows imposed conditions on the state variables following a spike, for example reset values for the voltage and fast potassium channel gating variables (e.g. for I_DR), and increments for slow potassium channel gating variables (for I_AHP and I_M), as appropriate. We have shown previously that these approximations are valid in a range of neuron models [17]. The associated threshold model compares quite well to the full conductance model with respect to interspike potentials and spike time moments.

The dynamics of a given population in the CBRD model are thus described by a set of equations of the form dY(t, t*)/dt = F(), where Y() includes the distribution ρ(t, t*) of neurons with a given value of t*, the average soma voltage U(t, t*), the average dendrite voltage U_D(t, t*) for two compartment models, and finally the average gating particle states x(t, t*) for HH cell models. Since the refractory time t* between spikes is proportional to t, the total time derivative dY(t, t*)/dt may be replaced by a sum of partial derivatives in t and t*, i.e. d/dt = ∂/∂t + ∂/∂t*.

The source term in the right-hand side of the equation for ρ(t, t*) is the fraction of neurons crossing the threshold, given by ρ(t, t*), multiplied by the hazard function H:

The remaining terms consist of the model details. The expressions corresponding to the cell voltages derive from those presented in Section 2.1.3, specifically the sets of Eqs 13 and 14, with the following replacements of the average voltages U, U_D and ξ, for V, V_D and 0, respectively:

The boundary conditions for ρ (which also defines the population firing rate, ν(t)) and the model voltages are:

where ΔtjAP is the duration of the action potential for cell type j (= 0 for the LIF model; for the HH cell models see Sections 4.3 and 4.4); likewise, the somatic reset voltage VjReset is defined also for each cell type.

The CBRD expressions for the gating particle states x(t, t*) of the HH channel models, including their boundary conditions, are presented in the next sections (Sections 4.3 and 4.4), which in turn reference the average somatic voltage U.

The dynamics of the entire neural population are found by the integration of Eqs 13 and 14, with the replacements indicated above for the average voltages U and U_D, which then defines the distribution of cell voltages across t* at each time t (for the HH models similar equations are solved for the gating particle state variables). The effect of threshold crossing in response to the input and noise is then taken into account by the H function, with the integration of Eq 55 giving the distribution of ρ across t* as well as the firing rate ν at time t.

As the source term in the density equation, the H function is the solution of the classical first-time threshold crossing problem for arbitrary history of stimulation. The H function has been semi-analytically derived from the KFP equation for voltage fluctuations, based on the following assumptions:

(i) Away from threshold, voltage fluctuations due to individual noise realizations can be described by a linear equation given the mean voltage U and the mean membrane conductance, which are in turn given by the associated transport equations and the synaptic input. Furthermore, the evolving probability density distribution of voltage fluctuations about U can be described by a KFP equation.

(ii) The flux across threshold described by H are due to two additive underlying processes: diffusion along the voltage axis due to noise, and transfer in response to the input.

(iii) Diffusion due to noise dominates the firing when U is constant, e.g. when the input is constant. In this case the problem is described by the KFP equation with a leak and a constant threshold, denoted here as A. For white noise the governing equation is reduced to an ordinary differential equation with an analytical solution for the spike generation probability expressed in Kummer functions [17]. In the case of correlated (colored) noise, A is obtained numerically and tabulated for a range of values for U and the time constant of the noise correlations τ [70].

(iv) Transfer based firing dominates when excitation starts abruptly, i.e. U increases with infinite rate. The initial condition is a stationary Gaussian distribution of neurons in terms of their potential V(t) − U(t). This implies that the state of zero distribution immediately following the previous spike is forgotten at the time of the next threshold crossing. As U(t) increases, the moving boundary at the threshold of fluctuations, V^Th − U(t), crosses the distribution. The flux through the boundary determines the probability of spike generation, which we denote in this case as B, and is obtained algebraically.

In our previous work we formulated an approximation for H for white Gaussian noise, as a function of the time-varying quantities that characterize the cell, including U (and its derivative with respect to time t at a given t*), σ_V, V^Th and the effective membrane time constant τ_m = C/g_tot, where g_tot is the total membrane conductance:

Note that the hazard function above has no free parameters, because they reflect the approximation of the first time passage solution. The dynamics of a neural continuum are thus evaluated at each time step by solving the set of one dimensional partial differential equations in terms of the refractory time, including the probability density of neurons, and the population averages of the membrane potential and channel gating variables. The CBRD model well approximates the firing rate of an infinite set of biophysically detailed neurons for an arbitrary stimulus, e.g. oscillatory input, and recurrent connections, when compared with Monte-Carlo simulations of individual neurons (see [17] and [70]). Fig 13 compares the response of a population of noisy LIF neurons, evaluated by Monte-Carlo simulations, a KFP approach (Section 4.1), and the CBRD method, with the response of the FR model with and without a time constant. Fig 14 compares the steady-state f/I characteristics as a function of synaptic conductance for a population of LIF neurons, with and without noise (Eq 5), with the result from the CBRD evaluation.