Bursting Neurons in the Hippocampal Formation Encode Features of LFP Rhythms

Burst spike patterns are common in regions of the hippocampal formation such as the subiculum and medial entorhinal cortex (MEC). Neurons in these areas are immersed in extracellular electrical potential fluctuations often recorded as the local field potential (LFP). LFP rhythms within different frequency bands are linked to different behavioral states. For example, delta rhythms are often associated with slow-wave sleep, inactivity and anesthesia; whereas theta rhythms are prominent during awake exploratory behavior and REM sleep. Recent evidence suggests that bursting neurons in the hippocampal formation can encode LFP features. We explored this hypothesis using a two-compartment model of a bursting pyramidal neuron driven by time-varying input signals containing spectral peaks at either delta or theta rhythms. The model predicted a neural code in which bursts represented the instantaneous value, phase, slope and amplitude of the driving signal both in their timing and size (spike number). To verify whether this code is employed in vivo, we examined electrophysiological recordings from the subiculum of anesthetized rats and the MEC of a behaving rat containing prevalent delta or theta rhythms, respectively. In both areas, we found bursting cells that encoded information about the instantaneous voltage, phase, slope and/or amplitude of the dominant LFP rhythm with essentially the same neural code as the simulated neurons. A fraction of the cells encoded part of the information in burst size, in agreement with model predictions. These results provide in-vivo evidence that the output of bursting neurons in the mammalian brain is tuned to features of the LFP.


Bursting neuron model
The equations and parameters of the bursting neuron model have been published before in Constantinou et al. (2015) and are included here for completion. The model consisted of two compartments: dendrites and soma. The input signal I(t) was injected in the dendritic compartment and bursting activity was recorded from the somatic compartment as described in Supplementary Equations 1 and 2, respectively.
The equations of the ionic currents and the model parameters are listed in Supplementary Tables 1 and 2, respectively. The 4 th order Runge-Kutta method with 0.01 ms time step was used for numerical integration of the model.

Discretization of LFP feature signals
LFP features (voltage, slope, phase and amplitude) vary continuously. However, estimating the probability of continuous LFP features from a finite sample would result in an enormous bias. Hence, the continuous signal was discretized into a finite number of bins M chosen by optimizing the trade-off between being large enough to preserve most of the information and, at the same time, small enough to reduce the bias, as in Elijah et al. (2015). In short, we varied M and, for each value, obtained the feature set X, from which we computed the mutual information I(X; N ) = I(X) − I(X|N ), the bias estimate I s (X; N ) = I(X) − I(X|N s ) (where ... indicates average over 100 repetitions), and the bias-corrected information I c (X; N ) = I(X; N ) − I s (X; N ). We computed these measures for the burst distinction code with no time lag, so that N = {1, 2, 3} and X comprised the features occurring at the time of burst onset.
The results were similar for all bursting neurons (example in Supplementary Figure 4). When M < 2 2 , the information measures were underestimated; whereas when M > 2 3 , there was considerable bias in the information estimate I(X; N ) as depicted by the diverging lines of I(X; N ) and the bias-corrected information estimate I c (X; N ). The bias correction method we used was sufficient to correct for this as indicated by the plateau of I c (X; N ) (Supplementary Figure 4). Since there was no considerable gain in information by using M > 2 2 , LFP features were discretized with M = 4.
In the full burst code: N = {0, 1, 2, 3}, in the burst rate code: N = {0, burst}, and in the burst distinction code: N = {1, 2, 3}. We define the variables Since the three variables arise from the same neural response, Y is a deterministic function of W , When Y = burst, the variable Z is also a deterministic function of W , The chain rule states that for any three variables X, A, B, I(X; A, B) = I(X; A) + I(X; B|A).
If we take A = W and B = Y , the chain rule becomes If we now take A = Y and B = W , the chain rule becomes In addition, x∈X w∈W P (x, w|y = 0) log P (x, w|y = 0) P (x|y = 0)P (w|y = 0) + (10) + P (y = burst) x∈X w∈W P (x, w|y = burst) log P (x, w|y = burst) P (x|y = burst)P (w|y = burst) Supplementary Equation 4 implies that when y = 0, there is no alternative but to have w = 0. Therefore, in the first term of Supplementary Equation 10, only the term with w = 0 appears, for all others, P (x, w|y = 0) vanishes. Moreover, P (x, w = 0|y = 0) = P (x|y = 0), so the first term vanishes. Similarly, Supplementary Equation 5 implies that when y = burst, w is equal to z, and P (x, w|y = burst) = P (x, z). Therefore, Supplementary Equation 9 becomes I(X; Y, W ) = I(X; Y ) + P (Y = burst)I(X; Z).

SUPPLEMENTARY RESULTS
Theta rhythms can be separated in two types based on their sensitivity to atropine (Kramis et al., 1975). Urethane preserves the atropine-sensitive theta (3-7 Hz) but eliminates the atropine-resistant theta (7-12 Hz) (Kramis et al., 1975;Clement et al., 2008). We observed that, under urethane anesthesia, the LFP of three rats exhibited shifts in which either delta (∼1 Hz) or theta rhythms (∼4 Hz) were dominant (example in Supplementary Figure 2). We identified 13 bursting units during the epochs containing dominant theta rhythms under anesthesia. These units were also bursting during epochs of dominant delta rhythms. We analyzed these data separately for the theta-dominant epochs and present the results here.
We investigated whether these cells encoded features of theta rhythms in their bursting output. Eleven cells showed evidence of encoding the instantaneous voltage, slope and phase of theta rhythms by the full burst code and the burst rate code (examples in Supplementary Figures 8A-B,D-E and 9A-D); and five of these cells also encoded the instantaneous amplitude. For the latter cells, the information encoded for voltage, slope and phase was twice to ten times higher than of amplitude. Four of the encoding cells also showed evidence of feature encoding by the burst distinction code (example in Supplementary Figures  8C and 9E-H). One cell encoded most information about the instantaneous amplitude of theta rhythms, and less about the voltage and slope, by the full burst code and burst rate code. These results suggest that bursting neurons can encode information conveyed by atropine-sensitive theta rhythms during anesthesia. Supplementary Table 1. Equations describing the ionic currents of the two-compartment model. The last row shows the kinetics equation of the gating variables for the Na, K and slow K currents.

SUPPLEMENTARY TABLES AND FIGURES
Leak currents