Previously, Guet-McCreight and Skinner (2021) applied synaptic perturbations to their detailed multi-compartment model, which as noted above is essentially the same as FULL and so we will not distinguish them moving forward. These perturbations were mostly done using 30 (inhibitory or excitatory) synapses randomly distributed across the dendritic tree as a reasonable capturing of experimental observations. Synaptic weights had been optimized based on experiment (Guet-McCreight and Skinner, 2020). Since SINGLE only has one compartment, perturbation was set to be 30 times the individual synaptic weight to maintain an approximately equivalent response in SINGLE relative to FULL. Synaptic equations, implementation and parameter values are identical to Equation 1 in Guet-McCreight and Skinner (2020) except for the synaptic weight values that are 0.006 (excitatory) and 0.0054 (inhibitory) μS here.

The mathematical model structure of SINGLE may also be used in data assimilation methods to estimate parameters based on an “observation”. Application of the RAUKF method (Azzalini et al., 2022) here was developed based on the reduced OLM model cell, and the model parameters were estimated from simulated noisy observations of SINGLE or FULL. The RAUKF method was adapted so that the conductance based model used for state prediction was the same as that used in SINGLE (see details in Supplementary material). The general formulation of the Kalman Filter (KF) is:

Where the state at k, x_k depends upon the previous state x_k−1 and some control input u_k−1, i.e., injected current, and mx~N(0,Qx). The function f is the system's state model, which in this case are the set of equations governing model SINGLE. In order to estimate parameters of the system, additional trivial dynamics were added to the system for each given parameter θ via a random walk:

With mθ~N(0,Qθ). Combining x and θ we may produce a new state X with the dynamics represented by F and noise M, where I is the identity function:

The observation model used to characterize noisy membrane voltage measurements is described by

Where n_{V, k} denotes measurement noise (in the general case, yk=g(Xk)+nk, nk~N(0,R)). Here, the direct observation of the membrane voltage y_k mimics experimental recording techniques (e.g., current-clamp methods). With only a subset of X_k being measurable, the method presented in this study allows hidden states to be estimated.

The system's dynamics model F, an initial estimate of the state X₀, as well as the control input u were used to maximize the probability of an observation y afforded by the noise profiles of each.

To estimate the next state (15) a set of points (sigmapoints) were sampled according to (14), described in more detail in Azzalini et al. (2022). This process also updates the current covariance P of the process.

Since the precise noise profile Q of the state vector X is not well known for a given neuronal model, and its respective discrepancies to the observation, we employed an adaptive method used in previous work in order to update Q as well as the estimate of the observation noise R over time (Azzalini et al., 2022). See Supplementary material for further details.

The adaptation of Q and R can be tuned by user defined variables a, b, λ₀, and δ₀, Where λ₀ and δ₀ are the minimum weights used when updating Q and R, the percent shift from the current to the new estimated value. To ease computational load these updates to Q and R were not applied after every sample, instead they were applied when a fault was detected (see Supplementary material). The values of λ₀ and δ₀ were both set to a value of 0.2. The variables a and b are weights of how much Q and R will be updated proportional to the magnitude of fault detected. The values of a and b were set according to explorations done in previous work when RAUKF was applied to a conductance based model (Azzalini et al., 2022). a>>b was found to be effective and here we use a = 10 and b = 1.

When applying the RAUKF to an observation of either SINGLE or FULL, θ was set to estimate the same conductance values that were optimized in the reduced model SINGLE. That is,

To generate an observation, an input train of step currents was injected into either SINGLE or FULL models. The input train consisted of 500 ms of zero current, then a step to 30, 60, 90, 0, –30, –60, and –90 pA, repeating 4 times. The current was also mixed with a Gaussian white noise ~N(4,5) pA. The mean of the noisy stimulation was chosen to match the characterized bias current in FULL and SINGLE. This input train of multiple positive and negative pulses is designed to be able to expose a full dynamic range of spiking, resting and subthreshold activities due to the underlying biophysical currents, and can be used in real-time to estimate parameter values.

The initialization of P and Q was derived from the order of magnitude θ, while x was made constant to 1 × 10⁻⁴ and 1 × 10⁻⁸, respectively. For the P of θ the value was set to 2 × the magnitude of the value of θ used in SINGLE, e.g., 340 → 1 × 10² → 1 × 10⁴, and 4 × the order of magnitude was used to initialize Q. Two initial states were used for θ₀, one where θ₀ = θ_SINGLE, known as the optimized initial estimate, and another where θ0=10⌊log10(θSINGLE)⌋, known as the poor initial estimate.

In vivo recordings of OLM cells exhibit a range of baseline firing frequencies that span ≈3 − 25 Hz (Varga et al., 2012; Katona et al., 2014) (depending on behavioral state of movement, sleep etc.). However, our previous spiking resonance explorations in IVL states using full, detailed multi-compartment models were technically constrained to a small range of baseline firing frequencies, specifically ≈12 − 16 Hz (Guet-McCreight and Skinner, 2021), due to how different IVL states could be generated. As such, we were limited in being able to extract biophysical underpinnings of different spike frequency resonances that include theta frequencies. Here, with SINGLE, our reduced model that maintains biophysical balances for the 5 different currents it contains (see Figure 1), we are able to obtain a wide range of baseline firing rates and thus can greatly expand our explorations in IVL states. This puts us in a position to determine what the biophysical determinants underlying theta frequency spiking resonance might be.

We obtained a wide range of baseline firing rates with SINGLE by altering the parameters for the “noisy” excitatory and inhibitory currents it receives (see Section 2). From our parameter search, 414 different sets of parameters that can generate valid IVL states were found. Within these sets, the excitatory conductance and variance were generally smaller than their inhibitory counterparts, suggesting that inputs received by OLM cells in vivo may be inhibitory dominant. From these 414 sets, we chose 10 that could span the full range of in vivo firing frequencies of OLM cells, and refer to them as 10 representative parameter sets. Each representative set has a minimum and maximum baseline firing rate as given in Table 1 along with its parameter set values.

For each representative set, we generated fifty 10 s long IVL states (via 50 random seeds) to obtain 500 IVL states. In Figure 5, we plot the minimum and maximum firing rates of these 500 IVL states. To illustrate what the OLM cell firings in IVL states look like, one example of an IVL state from each of the 10 chosen representative sets is also shown in Figure 5.

We can generate as many IVL states as we desire simply by using >50 random seeds. This was also the case in our previous work using a full multi-compartment model, but the randomness was in regard to synapse placement and presynaptic spike times from specified presynaptic cell types, resulting in a limited baseline firing frequency range (Guet-McCreight and Skinner, 2021). Here, the randomness is less restrictive because of the single compartment model nature of SINGLE and is able to exhibit the needed full range of baseline firing frequencies.

We are not limited in how many random seeds we can use. We already know that we can capture the full extent of in vivo firing ranges, and we aimed to generate a large data set of spike frequency resonances. With a large data set, we anticipated that we could determine whether there are any particular biophysical determinants underlying theta frequency spiking resonance specifically. We generated 50,000 IVL states from the 10 representative sets using 5,000 random seeds. In Figure 6A, we show the computed spiking resonance (f_r) as a function of the baseline firing frequency (f_B) when we used inhibitory perturbations. The analogous plot for excitatory perturbations is shown in Figure 6B. To ensure that the individual dots are visible, we only used the first 50 seeds (500 IVL states) in creating the plots shown in Figure 6. As can be seen, f_B spans the range, unlike the case with FULL (see Supplementary Figure S6). In Supplementary Figure S7, we plot the mean (Vm¯) and variance (σVm2) of the sub-threshold membrane potential, and the ISICV to show that these aspects of IVL states are not determining factors of the resulting f_r's. That is, for any given f_r, there are a range of Vm¯, σVm2 and ISICV values. Equivalently, resonant frequencies are distributed across the ranges of Vm¯, σVm2 and ISICV values.

Not surprisingly, one can see some dependence of f_r on f_B. That is, with a higher frequency f_B, there is an increase in the density of higher frequency f_r's. This can be seen in Figure 6, as well as in the density of colored dots representing f_B in Supplementary Figure S7. That is, IVL states with a higher f_B are more likely to have higher f_r, even though due to stochasticity of the IVL states, it is possible for some higher f_B IVL states to be resonant at low frequencies. However, this observation holds true on average. In essence, we have f_r's at many different frequencies.

Using FULL, we had previously noted that inhibitory perturbations, and not excitatory perturbations, generated f_r's that include theta frequencies during IVL states (Guet-McCreight and Skinner, 2021). Here, with SINGLE, this observation can be more clearly seen, as shown in Figure 6 where dashed red lines to delineate the theta frequency range (3–12 Hz) are drawn. With our greatly extended dataset, we can easily see that the full theta frequency range of f_r with inhibitory perturbations is exhibited. However, with excitatory perturbations, this is not the case—only minimally and at the upper end of theta frequency range are f_r's obtained. We note that the limited baseline frequencies available with FULL prevented us from being able to fully observe this.

As we would like to understand what underlies the ability of OLM cells to exhibit spiking resonance at theta frequencies, we now focus on just the inhibitory perturbations as it is inhibitory, but not excitatory, that encompass the full range of theta frequency spiking resonances (see Figure 6).

To consider this, we undertake a spike-triggered average (STA) analysis from the perspective of all of the underlying biophysical currents. Specifically, 10 different firing rates are chosen (set # 0–9 as shown in Figure 5). As noted in the Methods, since we have a mathematical model in hand, we are able to examine how the biophysical currents preceding spiking contribute for a range of baseline firing rates and spiking resonant frequencies. By performing STA analyses up to 200 ms prior to spiking, particular contributions by the various biophysical currents that influence theta frequency spiking (i.e., approximately 200 ms period) could be exposed. We examined each current type's contribution to spiking in IVL states for all of the different f_r's (27 in total from 0.5 to 30 Hz). In Figure 7A, we show STA plots for six different f_r's from these 27. For each f_r, we show normalized values in the STA plots for the five different biophysical currents (and leak) for each of the 10 representative sets (see Methods for details).

From STA plots using un-normalized values, the proportions of currents are directly comparable and it is clear, for example, that the relative proportion of Ih (to the total current) is much smaller than that of IM. From a thorough examination of STA plots of the un-normalized current values (proportion of the total current), we find that the larger the f_r, the more negative is the Ih and the less positive is the IM. In addition, STA plots with higher f_r frequencies show obvious fluctuations. These fluctuations would seem to be due to the presence of more perturbations, with the peaks and valleys of the fluctuations matching the timing of the perturbations. However, IM appears the least sensitive to the perturbations, with negligible fluctuations (see Figure 7A). Of particular note is that STA plots of IM and Ih show minimum overlap between the different representative sets relative to the other currents. This is more obviously the case for IM whose overlap is visually separable (see Figure 7A). In general, overlap is likely the result of stochasticity in the OU process which makes the currents variable and intersecting with each other. The minimal overlap in IM and Ih indicates an insensitivity to local changes, which could be attributed to them having slower time constants relative to the others. This is especially true for IM which has the slowest time constant (see equations and plots for the biophysical currents in Supplementary material). This would also explain why it does not show fluctuations for larger f_r whereas Ih does. A general interpretation of this minimal overlap is that regardless of what the in vivo firing rate of the OLM cell might be, IM and Ih can “tailor” their response accordingly. That is, they can perceive “global changes” and so the different firing rates are in the “history” of IM and Ih more than the other current types with faster time constants. Therefore, we now reasonably assume that if there is going to be any frequency specificity in the OLM cell's spiking resonance, Ih and IM would be the currents that could potentially control this. We now focus on Ih and IM to determine whether there are any particularities associated with theta frequency spiking resonance. To do this, we use an analysis that involves slopes and spreads. See Methods for computation details.

Let us first consider our slope analyses. We found that as f_r increased, the slope of IM's STA plot became less negative, and the slope of Ih's became less positive. Linear fits are statistically significant. In considering a combination of IM and Ih slopes, we again obtained a statistically significant linear fit, but now the slope crossed zero within the theta frequency range. This is all shown in Figure 8A where the 3–12 Hz theta frequency range is delineated by dashed green lines.

In our spread analysis, we found that as f_r increased, the spread in both IM and Ih STA plots increased as shown in Figure 8B. However, when we examined a combination of Ih and IM, the spread was more variable as shown in Figure 8C. Parabolic fits to each of IM, Ih and their combination were performed and the curve fits are shown in Figures 8B, C. From these fits, it is clear that the combined IM and Ih exhibits a tendency to have a minimal value in the theta frequency range, unlike IM or Ih alone, although none of these parabolic fits were found to be statistically significant. However, one can obtain a statistically significant linear fit to Ih (R² = 83.1%). While the parabolic fits were not statistically significant, a nonlinear local polynomial regression can yield a statistically significant fit with a minimal value for the combined IM and Ih, but neither IM nor Ih alone. Overall, the findings of our spread and slope analyses imply that to specifically have spiking resonance at theta frequencies, we need to have a balance. That is, our analyses indicate that there is a "balanced sweet spot" during which the combination of Ih and IM minimally change approaching spiking (zero slope and minimal spread), regardless of the baseline firing frequency of the OLM cell in vivo. This suggests that the slow time constants of IM and Ih act in concert to allow a theta frequency spiking resonance to emerge despite variability in OLM baseline firing rates. That is, IM decreases and Ih increases prior to spiking in such a way to favor a response to incoming theta frequency perturbations. The combination of their biophysical characteristics is what allows this, and not either current on its own.