As in Soudry and Meir (2014) we examine a single, synaptically isolated, excitable neuron under “spike” stimulation. In this stimulation regime, the stimulation current, I(t), consists of a train of identical short pulses arriving at times t_m and amplitude I₀. The intervals between the stimulation times are denoted T_m ≜ t_{m + 1} − t_m (Figure 1A, top). We assume that the stimulation is sparse, i.e., T_m ≫ τ_AP, with τ_AP being the timescale of an AP (Figure 1B). Since the neuron is “excitable” it does not generate APs unless stimulated, as in Gal et al. (2010) (i.e., the neuron is neither oscillatory nor spontaneously firing). However, after a stimulation the neuron can either respond with a detectable AP or not respond. We denote AP occurrences a Y_m, where Y_m = 1 if an AP occurred immediately after the m-th stimulation, and 0 otherwise (Figure 1A, bottom). Note also that Y_m is not generally the same as the common count process generated from the APs by binning them into equally sized bins (Appendix B.1) – unless T_m is constant and equal to the bin size.

We assume both Y_m and T_m are wide-sense stationary (Papoulis and Pillai, 1965). We denote p_* ≜ 〈 Y_m 〉 to be the mean probability to generate an AP and T_* ≜ 〈 T_m 〉 as the mean stimulation period. Furthermore, we denote Ŷ_m ≜ Y_m − p_* and T^_m = T_m − T_* as the perturbations of Y_m and T_m from their means. An important tool in quantifying the statistics of signals is the power spectral density (PSD), namely the Fourier transform of the auto-covariance (Papoulis and Pillai, 1965). For analytical convenience, in this work we will use a PSD of the form

with 0 ≤ f ≪ T⁻¹_* in Hertz frequency units. Note that this PSD is proportional to the PSD of the common binned AP (Equation 70), under periodical stimulus and for low frequencies – which is the regime under which we will investigate the PSD (similarly to the experiment Gal et al., 2010). We similarly define the PSD

and the cross-PSD

Typically, CBMs (Equations 4–6) contain many unknown parameters, and are highly non-linear. Therefore, it is quite hard to fit them using a purely simulation based approach, especially over long timescales, where simulations are long and models have more unknown parameters. Therefore, we developed a reduction method that simplifies analysis and enables fitting of such models. We refer the reader to Soudry and Meir (2014) for full mathematical details.

In this method, we semi-analytically² reduce the full model (Equations 4–6) to a simplified model, under the assumption that the timescales of rapid and slow variables are well separated. Given another assumption, that the neuron dynamics are sufficiently “noisy,” we can linearize the model dynamics, so that (Soudry and Meir, 2014, Equation 12)

where e_m is a white-noise signal with zero mean and variance σ²_e ≜ p_* − p²_* (recall p_* is the mean probability to generate an AP) and s_* (the excitability fixed point) and w_j (an “effective weight” of component s_j) can be found self-consistently together with p_* as a function of T_* (Soudry and Meir, 2014, Equation 10). After these quantities are found, an expression for the output PSD S_Y (f) in this model can be written explicitly. We let X₊, X₋ and X₀ denote the averages of the quantity X_s during an AP response, a failed AP response and rest, respectively. Also, we denote

as the steady state mean value of X_s [this would be X(p_*, T_*) in Equation 7 in Soudry and Meir, 2014]. For example, A_* and D_* are the respective steady state means of A_s and D_s. Additionally, we denote (definition below Equation 12 in Soudry and Meir, 2014)

as a “feedback” vector (see Figure 1C in Soudry and Meir (2014) to understand this interpretation), and (Soudry and Meir, 2014, Equation 14)

as the “closed loop transfer function” (including the effect of the feedback), with I being the identity matrix. Using the above notation, we can derive the PSD of the response. Given a periodical stimulation (T^_m = 0) we obtain (Soudry and Meir, 2014, Equation 13)

Though Equation (10) relies on two simplifying assumptions, extensive numerical simulations (Soudry and Meir, 2014, Figures 3–5) showed that this expression is rather robust and remains accurate in many cases even if these assumptions do not hold. Therefore, in this work we will always assume that Equation (10) is accurate.

In the neuron, the slow excitability variables s affect the response of the neuron, which, in turn, affects the dynamics of the slow excitability variables. To simplify analysis, it is desirable to “isolate” this feedback effect. In order to do this, we apply the Sherman Morrison lemma to Equation (9),

with

being the “open loop” version of H_c (f) (i.e., if a is set to zero). Using this in Equation (10) we obtain

with S^o_Y (f) being the “open loop” version of S_Y (f) (i.e., S_Y (f) with a set to zero),

and κ (f) determines the effect of the feedback

Note that κ (f) depends on the feedback through the variable a. If a → 0, for example, the kinetic rates of s are not sensitive to AP occurrences³. In that case κ (f) → 1 and S_Y (f) → S^o_Y (f).

In order to simplify analysis, we decompose the vector expressions in Equations (13, 14) to partial fractions.

If A_* is diagonalizable, than we can write Equation (13) as (Appendix A.1)

where the poles λ_k are the inverse timescales of the slow variables (the eigenvalues of A_*), arranged from large to small according to their magnitudes (0 < |λ_M| < |λ_{M − 1}| < … < |λ₁|) and

being the amplitude of these poles, with D_kj and w_k being the respective components of D_* and w in a basis in which A_* is diagonal. Note that, ∀k, Re [λ_k] < 0 (from the properties of A_*).

Using a similar derivation for κ (f), we obtain

with a_k and w_k being the respective components of a and w in a base in which A_* is diagonal.

For concreteness, we demonstrate our results on a simple model in which A_* is a diagonal matrix and, as a result, D_* (which depends on A_*) is also diagonal. In this “diagonal” model all the components of s are uncoupled (i.e., belong to different channel types), Equation (6) can be written as (Soudry and Meir, 2014, section 4.1)

∀ k∈{1, …, M}, where σ_s,k(V, s) = [(δ_k(V)(1 − s_k) + γ_k(V)s_k) N ^{− 1}_s,k]^1/2 and N_s,k are the number of slow ion channels of type k. Similarly as before, γ_+,k,γ_−,k and γ_0,k denote the averages of the kinetic rate γ_k(V) during an AP response, a failed AP response and rest, respectively. In addition

is the average γ_k(V) in steady state. We use a similar notation for δ. Therefore

with zero on all other (non-diagonal) components and

Therefore, in Equations (15, 16) we have,

Importantly, by tuning the parameters M, γ_k(V), δ_k(V), N_s,k and w_k we seem to have complete freedom in determining λ_k, c_k and a_k (Equations 19–21). This, in turn, would give complete freedom in tuning S^o_Y (f) and κ (f). Therefore, it seems that for any CBM (i.e., not only diagonal models) we can find an equivalent diagonal model – which produces exactly the same PSD of the response.

The only caveat in the previous argument is that in non-diagonal models λ_k can be complex, but not in a diagonal model, since the kinetic rates γ_k(V) and δ_k(V) must be real numbers. How would the situation change if some of the poles had complex values? Complex poles (i.e., for which Im[λ_k] > 0) always come in conjugate pairs. These pairs behave asymptotically (i.e., for 2π f ≫ |λ_k| or 2π f ≪ |λ_k|) very similarly to two real poles, with an additional “resonance” (either a bump or depression) in a narrow range in the vicinity of these poles (i.e., 2π f ~ |λ_k|) (see Appendix A.2, or Oppenheim et al., 1983).

In this section we derive general conditions on the parameters of a CBM (section 2.2) so it can generate robust f^−α statistics in S_Y (f). Here, we focus on the case of sparse periodical input T_m = T_* ≫ τ_AP (as in Gal et al., 2010).

This analysis is based on the decomposition of the PSD S_Y (f) as a ratio of S^o_Y (f) and the feedback term |κ (f)|². Recall that S_Y (f) ∝ f^−α is robustly observed for all stimulation parameters – even when p_* is near 0 or 1 (see section 3.1). Note that one can arbitrarily vary p_* by changing the stimulation parameters (such as I₀ or T_*). It is straightforward to show that when p_* → 0 or p_* → 1, the effect of feedback is negligible⁸, and therefore S_Y (f) ≈ S^o_Y (f). This implies that, at least for some simulation parameters, S^o_Y (f) ∝ f^−α. For this reason, and for the sake of analytical simplicity, we first develop general conditions so that S^o_Y (f) ∝ f^−α, and later we discuss the effects of the feedback κ (f).

Note from Equation (15) that if M (the dimension of s – the number of slow processes) is finite, one can have S^o_Y (f) ∝ f^−α exactly if and only if α = 0 or 2. However, these values are far from what was measured experimentally (Equation 42). Therefore, S^o_Y (f) ∝ f^−α can be generated exactly only in some limit (in which M is infinite), or approximately (if M is finite). Also, note that if 2π f ≫|λ₁|, then S^o_Y (f) − T_* σ²_e ∝ f⁻². Additionally, if 2π f ≪ |λ_M|, we have S^o_Y (f) ≈ constant. Therefore, Equation (15) can generate S^o_Y (f) ∝ f^−α with 0 < α < 2 only for |λ_M| < 2π f < |λ₁|.

Next, we explain when this becomes possible. For simplicity assume that in Equation (15) T_* σ²_e is negligible and all the poles are real (the effect of complex poles will be discussed below). We define the following pole density

where δ(·) is Dirac's delta function. Using Equations (15 and 23) we obtain

For |λ_M| ≪ 2π f ≪ |λ₁| and 0 < α < 2, Equation (24) becomes

if and only if (Appendix A.3)

in the range |λ₁| > |λ| > |λ_M|, with ρ₀ = 2π⁻¹C sin(πα/2). Therefore, ρ(λ), the distribution of the poles, must scale similarly to S^o_Y (f) (but with a different exponent).

Several comments are in order at this point.

In the previous section we found general conditions under which Equation (13) gives S^o_Y (f) ∝ f^−α. In this section, we aim is to generate S^o_Y (f) ∝ f^−α over f_min < f < f_max in a minimal model, in which M (the dimension of s) is as small as possible. As explained in section 2.3.4 we do not lose any relevant generality if we restrict ourselves to the case where A_* is diagonal (Equation 18). From Equation (26), we know that |λ_k| must “cover” the frequency range f_min < f < f_max. In order for M to be small, we choose λ_k to be uniform over a logarithmic scale (similarly to Keshner, 1982), so λ_k ∝ ϵ^k with ϵ < 1. The “simplest” way to achieve this is to have (see Equation 18)

so

In order for λ_k/(2π) to cover the range [f_min,f_max] we require that

Given M, this sets a constraint on ϵ. In order to have scaling in ρ(λ), as in Equation (26), we also require that c_k ∝ |λ |^1−α dλ ∝ ϵ^(2−α)k, since d λ = λ_k − λ_{k − 1} ∝ ϵ^k. Therefore, from Equations (21) and (20) we have

so that S^o_Y (f) ∝ f^−α. Therefore, we require that w_k ∝ ϵ^{−μ k}, N_s,k ∝ ϵ^νk with 2μ + ν = α − 1. For μ > 0 the slower processes (larger k) have larger weight. For ν > 0 slower processes have a smaller number of ion channels (therefore, they are more “noisy”).

In Appendix A.4, we investigate what type of scaling will generate also S_Y (f) ∝ f^−α, taking into account the effects of feedback (through κ (f)). We conclude that, because of the feedback, a value of μ > 0 would not change the exponent of S_Y (f) over a “reasonable” range of parameters (i.e., assuming ν > −2). Therefore, the simplest way to generate S_Y (f) ∝ f^−α would be to take μ = 0. In this case, we have (Equation 59), for −1 < ν < 1 and |λ_M| ≪ 2π f ≪ |λ₁|,

where the logarithmic correction arises from the effect of feedback κ (f). A few comments on Equation (30) are in order at this point.

Recall that T^_m ≜ T_m − T_*, with T_* ≜ 〈 T_m 〉 and S_T (f) the PSD of T_m (Equation 2). As explained in Soudry and Meir (2014), for a general CBM⁹ we can decompose Ŷ_m, the fluctuations in the neuronal response, to a linear sum of the history of the input and internal noise, i.e.,

with the filter h^ext_k used to integrate external fluctuations in the inputs, and the filter h^int_k used to integrate z_m, a zero mean and unit variance white noise representing internal fluctuations (e.g., ion channel noise). It is easier to analyze this I/O in the frequency domain, where Equation (32) becomes (Soudry and Meir, 2014, Equation 20)

where we define X(f) to be the Fourier transform of X(t). Together, H^ext(f) and H^int(f) describe the T^_m → Ŷ_m neuronal I/O at very long timescales.

Note that these filters are related to the PSDs, in the following way

where we recall that S_YT (f) is the cross-PSD (Equation 3). Notably, from Equation (35), if the input to the neuron is not periodical (so, S_T (f) ≠ 0), then the PSD S_Y (f) should be the same as calculated previously, except for the addition of |H^ext(f)|²S_T (f).

For a general CBM, we can derive semi-analytically the exact form of the filters in Equation (33) from its parameters, as we did for S_Y (f). For example, if T^_m = 0 (periodical input), then also S_T (f) = 0, and so

where S_Y (f) is the PSD we derived previously (Equation 10). Additionally, we obtain (Soudry and Meir, 2014, Equation 17)

with

Next, we find both filters for the minimal model described in section 3.2.2. Recall that in this model

with a₁ and d₁, respectively given by Equations (8) and (38). To simplify analysis, we derive an asymptotic form for both filters, for the cases |λ_M| ≪ 2π f ≪ |λ₁| and 2π f ≫ |λ₁|. First, from Equation (36), and Equation (59), we find

Similarly, from Equation (37), we find (Appendix A.6) that for the minimal model the interpolation between the two asymptotic cases is monotonic, so we can approximate

where q ≜ (1 − ϵ)⁻¹T⁻¹_*w₁. A few comments on Equations (40, 41) are in order at this point.

The experiment from Gal et al. (2010), where a single synaptically isolated neuron, residing in a culture of rat cortical neurons, is stimulated periodically with a train of extracellular short current pulses with constant amplitude I₀. The observed neuronal response was characterized by different modes (Gal et al., 2010, Figure 2). We focus on the “intermittent mode” steady state, in which 0 < p_* < 1 (i.e., sometimes the stimulation evokes an AP, and sometimes it does not). The patterns observed in Y_m, the AP occurrences timeseries, are rather irregular (Gal et al., 2010, Figure 2E), span multiple timescales (Gal et al., 2010, Figure 5) and variable (i.e., patterns are not repeatable Gal et al., 2010, Figure 9A). More quantitatively, as indicated by the analysis (Gal et al., 2010, Figure 6), for all intermittently firing neurons, the patterns in Y_m fall into the category of “f^−α noise” where the value of α varied significantly between neurons – with

(mean ± SD). As we explained in section 3.1, this f^−α behavior is true only in some limited range f_min < f < f_max. From the experimental data, (Figure 6C in Gal et al., 2010) it can be estimated that f_min < 10⁻⁵Hz and f_max ~ 10⁻²Hz. Also, since α > 1, then 0 < f_min (see section 3.1).

In our previous work (Soudry and Meir, 2012) we already fitted a model that fits many of the “mean” properties of the neuronal response (e.g., firing modes, transients and firing rate). This model is an extension of the original Hodgkin–Huxley model which includes Slow sodium inactivation (Chandler and Meves, 1970; Fleidervish et al., 1996) (The HHS model, see Appendix C.1). In order to maintain this fit with the experimental results, we extend the HHS model with additional slow components, obeying Equation (18). We denote this as the HHMS model (Hodgkin Huxley model with Many Slow variables, Appendix C.2). The equations are identical to the HHS model, except that in the voltage equation (Equation 73) g_Nas is replaced by g_NaM⁻¹ ∑i = 1Msi, where s₁ has the same equation as s in the HHS model (Equation 77). By symmetry, this gives identical weights to component s_i (i.e., ∀k: w_k = w₁). The remaining rates (for k ≥ 2) are chosen according to our constraints, so γ_k(V) = γ(V) ϵ^{k − 1}, δ_k(V) = δ(V) ϵ^{k − 1} (as in Equation 27), where γ(V) and δ(V) are taken from the HHS model (Equation 77) and also N_s,k = N_s ϵ^{ν k}. Therefore, the only free parameters are ϵ, M, N_s, ν and I₀ (I₀ is the current amplitude of the stimulation pulses).

This model can be used to fit the experimental results for any α ∈ [0,2). We performed a numerical simulation of the full equations (Equations 4–6) of the HHMS model under periodical stimulation with T_* = 50 ms. We aimed to fit an experiment from Gal et al. (2010), which had a similar stimulation and exhibited S_Y (f) ~ f^−α, with α = 1.4 (which is approximately the average α value measured in Gal et al., 2010). The current amplitude I₀ was set to I₀ = 7.7 μ A so that the model would have the same mean response probability p_* ≈ 0.4 as in the experimental data (using the self consistent equations for p_* from Soudry and Meir, 2014). We chose M = 5 and ϵ = 0.2 in order to satisfy constraint Equation (29) with a minimal M. We chose ν = 0.5 to satisfy Equation (31). Lastly, we chose N_s = 10⁴ in order to fit the magnitude of the S_Y (f). This reproduced all the scaling relations observed experimentally (Figure 3).

After fitting the HHMS model to the experimental results, we can examine its resulting linearized input–output relation, described by the filters H^ext(f) and H^int(f) (Equation 33). The H^int(f) filter integrates internal fluctuations, while the H^ext(f) filter determines how external fluctuations (in the input) affect its response.

In accordance with the asymptotic forms in Equations (40) and (41), we find that H^ext(f) is a low pass filter with a pole f_ext ~ 10⁻² Hz (Figure 4, green) while H^int(f) ~ f^−α/2 for f_min < f < 10⁻² Hz (Figure 4, red) with f_min < 10⁻⁵ Hz. Therefore, as explained in section 3.3.2 this model implies that the response of the neuron is affected by internal fluctuations over the scale of days (~ f⁻¹_min) or more, generating the f^−α behavior we observe in Figure 3. However, external input is remembered only for minutes (~ f⁻¹_ext).

Next, we examine two methods which allow us to probe H^ext(f) directly and examine these predictions.

First, a simple method to probe the external input filter H^ext(f) is through Equation (34). Allowing reliable estimation of H^ext(f) in a certain frequency range requires a random process stimulation for which |H^ext (f) |²S_T (f) ≫ |H^int(f)| in that range, as explained in Appendix B.2. To demonstrate this method we estimate S_TY (f) from the existing experimental data taken from Gal and Marom (2013), in which S_T (f) ~ f^−β (above some lower cutoff). In Figure 5A we compare this estimation with S_TY (f) in the HHMS model in a limited range where S_T (f) is sufficiently high for estimation to be accurate. It is similar to S_TY (f) from our fitted HHMS model, validating our estimate of H^ext(f) for that frequency range.

Second, The filter H^ext(f) could be probed more accurately and at lower frequencies – by sinusoidally modulating the input (the internal-stimulus intervals), analogously to the sinusoidally modulated input current used in Lundstrom et al. (2008, 2010) and Pozzorini et al. (2013),

As we explain in Appendix B.3, in this case the output of the neuron would be

This allows us to easily estimate |H^ext(f)| using the peaks of T_amp⁻¹ Ŷ (f) (the Fourier transform of T⁻¹_amp Ŷ_m) at frequencies f_l, as we demonstrate in Figure 5B, using our fitted HHMS model.

As explained in Gal et al. (2010) and Soudry and Meir (2012) the latency of the AP can serve as an indicator of the cell's excitability. Specifically, this is true in the HHMS model, for periodical stimulus and p_* = 1, where the PSD of the latency, S_L (f), is a shifted and scaled version of S_Y (f) with p_* → 1 (see section 4.4.6 in Soudry and Meir, 2014). Therefore, in the HHMS model we also have S_L (f) ∝ f^−α approximately (neglecting logarithmic factors).

Next, suppose we vary some measurable stimulation parameter, such as the mean stimulation rate T⁻¹_*. How would this affect the shape of the filters we derived? The analytical results allow us to calculate this explicitly in the HHMS model.

First, we consider the gain of the external input filter H^ext(f) (i.e., H^ext(0)). As we explain in Appendix A.7, if f ≪ f_cutoff, than

which is the mean firing rate of the neuron – a simply measurable quantity.

Second, how would H_int (f) change if T_* is varied? Since H_int (f) is directly measurable only through S_Y (f) (Equation 36), we consider S_Y (f) instead. From Equation (59) it is clear that if S_Y (f) ~ f^−α approximately at low frequencies then the exponent α should not depend much on any external parameter (assuming 0 < p_* < 1). This was observed experimentally when the stimulation rate (T⁻¹_*) was varied, as can be seen in Figure 1G in Gal and Marom (2013).