How to correctly quantify neuronal phase-response curves from noisy recordings

Abstract

At the level of individual neurons, various coding properties can be inferred from the input-output relationship of a cell. For small inputs, this relation is captured by the phase-response curve (PRC), which measures the effect of a small perturbation on the timing of the subsequent spike. Experimentally, however, an accurate experimental estimation of PRCs is challenging. Despite elaborate measurement efforts, experimental PRC estimates often cannot be related to those from modeling studies. In particular, experimental PRCs rarely resemble the characteristic theoretical PRC expected close to spike initiation, which is indicative of the underlying spike-onset bifurcation. Here, we show for conductance-based model neurons that the correspondence between theoretical and measured phase-response curve is lost when the stimuli used for the estimation are too large. In this case, the derived phase-response curve is distorted beyond recognition and takes on a generic shape that reflects the measurement protocol and masks the spike-onset bifurcation. We discuss how to identify appropriate stimulus strengths for perturbation and noise-stimulation methods, which permit to estimate PRCs that reliably reflect the spike-onset bifurcation – a task that is particularly difficult if a lower bound for the stimulus amplitude is dictated by prominent intrinsic neuronal noise.

Appendices

Appendix: Conductance-based neuron models

Conductance-based neuron models were based on established models from the literature. The Hopf model was first described by Morris and Lecar; we use the version from Ermentrout and Terman (2010), p. 50-51. The SNIC and HOM models are versions of the Wang-Buzsaki model (Wang and Buzsáki 1996). The membrane voltage v of the models follows the dynamics

$$ \begin{array}{@{}rcl@{}} \frac{\text{d}v}{\text{d}t} = (I_{\text{in}} + g_{\text{L}} (E_{\text{L}} - v) + I_{\text{gates}})/C_{\text{m}}. \end{array} $$

The input current is the sum of a DC current, a time-dependent stimulus and an intrinsic noise ξ, I_in = I_DC + s(t) + σξ. For simulations without intrinsic noise, σ = 0. The model parameters are summarized in Table 1.

Table 1 Parameters of the Hopf model and the SNIC model

Table 2 Stimulus amplitudes for different models and PRC estimation methods

Gating for the Hopf model:

For the Hopf model,

$$ \begin{array}{@{}rcl@{}} I_{\text{gates}}(v,n) = g_{\text{Na}} m_{\infty}(v) (E_{\text{Na}} - v) + g_{\text{K}} n (E_{\text{K}} - v). \end{array} $$

The kinetics of the ion channel gating n are given by

$$ \begin{array}{@{}rcl@{}} \frac{\text{d}n}{\text{d}t} = \phi \frac{n_{\infty}(v) - n}{\tau_{n}(v)}, \end{array} $$

with τ_n(v) = 1ms/cosh ((v/mV − 2)/60). The ion channel activation curves are given as

$$ \begin{array}{@{}rcl@{}} m_{\infty}(v) &=& 0.5(1+\tanh((v/\text{mV} - (-1.2))/18)), \\ n_{\infty}(v) &=& 0.5(1+\tanh((v/\text{mV} - 2)/30)). \end{array} $$

2.1 Gating for the SNIC and HOM models:

For the SNIC model,

$$ \begin{array}{@{}rcl@{}} I_{\text{gates}}(v,n,h) = g_{\text{Na}} m_{\infty}(v)^{3} h (E_{\text{Na}} - v) + g_{\text{K}} n^{4} (E_{\text{K}} - v). \end{array} $$

The kinetics of the ion channel gating are given by

$$ \begin{array}{@{}rcl@{}} \frac{\text{d}h}{\text{d}t} &=& \phi (\alpha_{\text{h}}(v) (1-h) - \beta_{\text{h}}(v) h),\\ \frac{\text{d}n}{\text{d}t} &=& \phi (\alpha_{\text{n}}(v) (1-n) - \beta_{\text{n}}(v) n), \end{array} $$

with

$$ \begin{array}{@{}rcl@{}} \alpha_{\text{h}}(v) &=& 0.07 \exp(-\frac{v/\text{mV}+58}{20})/\text{ms},\\ \beta_{\text{h}}(v) &=& \frac{1}{ 1 + \exp(-0.1 v/\text{mV}-2.8)} /\text{ms},\\ \alpha_{\text{n}}(v) &=& \frac{0.01 v/\text{mV}+0.34}{1 - \exp(-0.1 v/\text{mV}-3.4)}/\text{ms},\\ \beta_{\text{n}}(v) &=& 0.125 \exp(-\frac{v/\text{mV}+44}{80})/\text{ms},\\ \end{array} $$

The ion channel activation curve for the gating variable m is given as

$$ \begin{array}{@{}rcl@{}} m_{\infty}(v) &=& \frac{\alpha_{\text{m}}(v)}{\alpha_{\text{m}}(v) + 4 \exp(-\frac{v/\text{mV}+60}{18})},\\ \alpha_{\text{m}}(v) &=& \frac{0.1 v/\text{mV}+3.5}{1 - \exp(-0.1 v/\text{mV}-3.5)}. \end{array} $$

The HOM model is equivalent to the SNIC model, but with ϕ = 1.5 and I_in = 0.166μ A/cm², for more details see Hesse et al. (2017).

Conductance-based neuron models were simulated numerically with the simulation environment brian2 (Stimberg et al. 2014). To record 500 perturbed spikes, the recording duration is 50 seconds for the noise-stimulation methods, and 100 seconds for the perturbation method, where only every second spike is perturbed. We observed that the interspike interval depends on the time resolution dt of the simulation. We have chosen dt = 0.001ms for all models, to ensure that a three times smaller time step changed the deterministic interspike interval by less than 5%.

To investigate the tolerance of PRC measurements towards an unknown noise source, PRCs were in a second step measured for neuron models that included intrinsic noise, implemented as an additive zero-mean white-noise current (the brian2-implemented variable xi). The noise levels chosen in this study correspond to relatively strong intrinsic noise, with a phase noise of $\tilde {\sigma } = 2\sqrt {\text {ms}}$ or $\tilde {\sigma } = 3\sqrt {\text {ms}}$ standard deviation (Fig. 4). The next section shows how to translate the phase noise into the standard deviation of the brian2-implemented noise variable xi. The chosen phase noise results in a CV of around 0.2 and 0.3, respectively, which is within a biologically realistic range.

As a side-note, while the stimulus amplitude appropriate for PRC estimation is larger for perturbations than for a noise stimulus, the latter induce more spike jitter (larger CVs) compared to the perturbations. It seems that as the perturbation is temporally precise, even small spike deviations are sufficient to estimate PRCs, while the continuous input delivered by the noise stimulus requires larger deviations to estimate the PRC, probably because every spike informs about the full phase range, instead of just one particular phase.

2.2 Adaptation of intrinsic noise levels for brian2 models

Comparable levels of spike jitter between models can be ensured by balancing the effective phase noise. Let η(t) be a zero-mean white noise current with standard deviation 1. As intrinsic noise current of the neuron model, we chose ση(t) with standard deviation σ. The phase Eq. (1) is $\dot {\varphi }(t) = 1 + Z s(t)$ with current input s(t) and PRC Z. With the intrinsic noise as input, s(t) = ση(t), and under the assumption that the intrinsic noise is weak, and thus influences φ much less than the deterministic dynamics, one can replace the temporally precise filtering of the noise by its averaged effect in one interspike interval (Schleimer and Stemmler 2009). We get $\dot {\varphi }(t) = 1 + \tilde {\sigma } \eta (t)$, with the variance of the white noise as $\tilde {\sigma }^{2} = \sigma ^{2} {\int _{0}^{1}} Z^{2}(\varphi ) \text {d}\varphi$.

Setting the phase noise $\tilde {\sigma }$ to the same value in all models, we find the appropriate current noise strength as $\sigma = \tilde {\sigma } \frac {C_{\text {m}}}{T} \left [\int _{0}^{1} Z^{2}(\varphi ) \text {d}\varphi \right ]^{-0.5}$, where C_m is the membrane capacitance of the model. For the simulations, we evaluate this formula with the theoretical iPRCs gained from backwards integration of the adjoint equation as Z, see Section 2.3.

2.3 Error estimation by bootstrapping

In order to evaluate the quality of the PRC, we use a bootstrap approach to estimate errors for the phase response curves. Two different kinds of phase-dependent errors are considered, the baseline error which results from PRC estimates on shuffled data, and the error on the PRC estimate.

In order to provide an error for the PRC, we repetitively estimate PRCs using only a restricted amount of spikes. We estimated PRCs in 100 repetitions from a set about 250 spikes, randomly chosen from the total set of about 500 perturbed spikes. The standard deviation of the 100 PRC estimates was used as PRC error.

In order to provide a baseline for the PRC estimate – above which a significant PRC should rise – we measure PRCs for random combinations of noise snippets and spike advances. We measure the PRC in 100 repetitions for a random permutation of the numbering i of spike deviations {Δφ_i}, and calculate the standard deviation for the resulting 100 estimates. The mean of these estimates results in a PRC close to zero, such that the resulting standard deviation corresponds to the range within which a zero, i.e., non-significant PRC estimate is to be expected. This error around zero provides a lower bound to PRCs that are significantly different from zero.

2.4 Amplitude scaling for estimated phase-response curves

The appropriate unit for most experimentally derived PRCs is Hz cm²/μ A, which naturally arises for PRCs that measure the phase advance in response to a current per surface area, such as a noise stimulus or perturbation in μ A/cm². For this study, however, we aim at comparing experimentally measured PRCs with theoretical iPRCs. The theoretical iPRCs gained from backwards integration of the adjoint equation have a unit of 1/mV, as they quantify the phase advance in response to a voltage perturbation in units of mV. For comparison, we transform the estimated PRCs into the units of the theoretical iPRCs by multiplying the current PRC by the membrane capacitance, Z_Δv = Z_ΔIC_m