Local field potentials (LFPs) measured in the extracellular matrix of the brain are postulated to arise from the integration of synaptic ionic currents and spread by volume conduction. However, there is a lack of consensus on whether these spatiotemporal voltage gradients are just an epiphenomenon of spiking or if the LFPs play a functional role in information processing. To examine a potential functional role of LFPs in information processing, we developed a microfluidic device that allows neurons from the hippocampal formation to self-wire through microfluidic channels, effectively isolating the activity of single axons between subregions of the network. We recorded spontaneous theta-band activity (4-10 Hz) in these axons whose power spectra were independent of simultaneous spiking activity. The highest theta amplitudes above noise were seen intermittently in a sparse set of axons from the CA3 into the CA1. Source neurons for the axonal theta were identified through cross correlation. Functionally, sparse axonal theta phase and amplitude correlated with target subregional spiking and more strongly with burst length. These results suggest that theta voltage oscillations in axons may contribute to activation of slow voltage-gated calcium channels to drive stronger synaptic release of transmitter to coordinate hippocampal activity between subregions. We propose that theta oscillations are controlled by specific ion channels distinct from those that generate spikes, a multiplex coding mechanism for inter-regional communication with implications for routing, executive control, EEG, disease states, and artificial neural networks.

Four Chamber, Five Connection Tunnel Device

In brief, a microfluidic device was molded in PDMS. The device was comprised of four chambers for separate culture of subregions of the hippocampus. Each chamber was 30 mm². The subregional chambers were connected by 3 x 5 x 400 μm (HxWxL) microfluidic tunnels. In each array, 67 tunnels connected each subregion with four of those tunnels monitored by electrodes, excluding the EC-CA3 cross connection that had 11 tunnels with five tunnels monitored. These dimensions favor mostly one, two and three axons per tunnel (Narula 2017). Resistance within the tunnels was measured as R=ρL/(W×H), where ρ is the resistivity of the Neurobasal medium, 90 Ωcm. This yields a resistance for the tunnel of 24 MΩ, far higher than the subregional electrode resistance 0.2 MΩ. Therefore, these axonal tunnels produce high amplitude signals, well isolated from neighboring tunnels and the bulk subregion from which they originate. The device was affixed to a Multichannel Systems MEA120 1100 (Multichannel Systems, Reutlingen, Germany) for amplified recordings of spontaneous activity on the 120-electrode microarray.

In vitro hippocampal neuronal network culture

Neurons were isolated and cultured from postnatal day 4 Sprague Dawley rat pups under anesthesia as approved by the UC Irvine Institutional Animal Care and Use Committee (IACUC). Briefly, entorhinal cortex (EC), dentate gyrus and hilus (DG), CA3, and CA1 including the subiculum were microdissected and cultured for 21 days in NbActiv4 medium; Transnetyx BrainBits, Springfield, IL).

Acquisition and processing of electrophysiological activity, spike sorting and direction of axonal information transfer

Spontaneous activity was analyzed from six networks (n=6 arrays) using MC_Rack (Multichannel Systems) and custom MATLAB software at a sampling rate of 25 kHz at 37 °C for 5 min. over the range of 0.5-50 kHz. Arrays with less than 80% active tunnels or with poor growth in one of the subregions were rejected for recording. Methods for spike sorting and determining the direction of axonal communication were described in our previous publications. In brief, axonal spikes were sorted using Wave-Clus with sensitivity set to 5 S.D. of the signal for the tunnels and subregions. Feedback and feed-forward direction of spike propagation was determined by the temporal delay in propagation of action potential over two electrodes spaced at 200 μm using our normalized matching index (NMI) algorithm. Electrodes were considered for evaluation if the spike rate was at least 0.4 Hz averaged over the whole 300 second recording. We examined theta activity in axons via spectrograms generated by MATLAB’s continuous wavelet transform toolbox. After down sampling to 2500 Hz, the Morlet wavelet filter bank was set for the frequency ranges of 3-300 Hz at 10 voices per octave.

Spike contribution to theta power

We wanted to evaluate whether theta oscillations were simply an artifact of repeated spiking. Spike cutouts from recorded signals were repeated at theta burst frequencies and zero-padded to 2 seconds for spectrographic analysis. We compared the power in the spectrogram for a single spike, three spikes pasted at theta timing of 100 ms and three trains of five spikes spaced at 200 Hz with trains repeated at 5 Hz. The spectrogram of each signal was computed for 3-300 Hz and summed over 4-10 Hz for total theta power. The summed magnitudes of the signal were normalized for frequency to reduce the over representation of higher frequency signals. The Hilbert transform was used on the theta filtered signal to visually represent the power of each signal.

Prevalence of Axonal theta amplitudes

Theta amplitude was determined by filtering raw data collected from microelectrodes using a zero-phase 8^th order Butterworth filter for 4-10 Hz and down sampled to 1000 Hz.  A 100 log-binned histogram was used to count all Hilbert-transformed theta amplitudes in each single axon. Next, a two-component univariate Gaussian Mixed Model (GMM) was employed to segregate a noise class (µ₁) and high amplitude theta class (µ₂) in each distribution. The decision boundaries were calculated via Bayesian statistics. To normalize the Gaussian probability distributions, The joint probabilities are found by multiplying the likelihood parameter by the weight values of each class from the GMM and dividing by the sum of joint probabilities.  Each axonal high amplitude class mean (µ₂) was compared to the noise mean (µ₁) of all axons in that subregion via a Dunnett multiple comparisons test. The average class of µ₁ and significantly different µ₂ are reported for each subregion. ANOVA with HSD multiple comparisons adjustment were used to compare the average high amplitude theta class in each subregion.

Identifying Source-Axon-Target Relationships

The goal was to determine the neuronal source of the axon recorded in the tunnel and the corresponding subregional target neuron. To find the source, we iteratively determined which, if any, of the 18 electrodes in the upstream subregion recorded neuronal theta activity that correlated with axonal theta by MATLAB XCORR with a 200 ms window (one theta cycle). Significance of the correlation was calculated using MATLAB CORRCOEF. To find the target neuron(s), we determined which axonal theta amplitudes or phase angles significantly modulated spikes in the target during axonal theta amplitudes above 5 μV. The influence of axonal theta on the target was computed as a modulation index.

Significance of Axonal Phase Angle and Amplitude with Subregional Spiking

Subregional spikes were counted within thresholded periods when axonal theta amplitudes were above background noise (5 μV) and longer than 200 ms (2x the shortest theta cycle). Intervals were connected if they were less than 3 cycles apart. To establish a relationship between the angle of axonal theta with subregional spiking, the Tort Modulation Index (MI) was used to determine if distributions of instantaneous axonal angles at subregional spike times derived from the Hilbert transform deviated significantly from a uniform distribution. In short, a binned probability distribution of spikes over a range of angles was created and normalized for probability. The binning was determined via the Freedman-Diaconis rule (n=20) where the angle was on a linear scale from –180 degrees to 180 degrees. The Kullback-Leibler distance derived from the Shannon Entropy was divided by the log of the number of bins to get the Modulation Index. Nonparametric statistical testing was performed through shuffling high amplitude axonal angles 1000x via randomizing the imaginary component of the FFT against the subregional spikes in those same time intervals and a p-value was obtained. To establish a relationship between the amplitude of axonal theta with subregional spiking, a log-log binned probability distribution of subregional spikes over a range of axonal amplitudes was created. The binning was set equal to the number of angle bins (n=20). To obtain the significance of the amplitude measurements, a log-log linear regression was fit over this distribution for an R², slope, and p-value of the slope significantly different from zero. The optimal range for the regression in each case was found via a grid search for the highest R² value with the start of the regression anchored at the lowest amplitude (5 μV).

Mutual information between axonal voltage oscillations

For specificity among axons, we determined whether multiple axons carried the same theta signals. Mutual information was used to determine how much information was shared in axons communicating in the same direction. Binning was determined by the Freedman-Diaconis rule (n=20). The amplitudes of the axonal Hilbert envelopes of theta waves as the basis for comparison. The range spanned from the minimum amplitude of an axon to the maximum amplitude and was logarithmically spaced. This ensured that coupling of an axon to an electrode did not affect its ability to be compared to another axon, in effect normalizing the amplitudes between two axons. Distributions of binned amplitudes were compared to see how much information was in the same bin between every axon pair in a subregion, computed in bits.

General Linear Models for predicting subregional spiking

General linear models (GLM) were used to predict if the amplitude and phase features of theta in axons had an effect on subregional target neuron spiking. The MATLAB R2024 fitglm function was used to calculate a GLM for predicting subregional spiking during intervals of high amplitude theta (selected using thresholding method explained above). The GLM was computed on signals down sampled to 2500 Hz. The response variable was target neuron spiking. Spiking was coded as indices of logical 1’s. The logit link function was used to determine if the independent variables of theta amplitude and phase predicted target neuron spiking. Theta phase was transformed to the cosine of theta during high theta episodes. Cosine was chosen because it linearized the cyclical relationship of angle and had a stronger prediction of spiking than either angle in degrees or the sine of the angle alone. All other indices below the threshold were not considered. Additionally, we computed GLMs for bursting activity in target neurons. All indices in the burst bounds were defined as “bursting” vs “non-bursting” in the indices outside the burst bounds. We defined the GLM using the formula in Wilkinson Notation:

Wellspikes or WellBursts ~ axonalThetaAmp * (CosThetaAngle)

ThetaAmp was the z-scored amplitude of theta at each index above threshold. We defined bursts as a minimum of four spikes with a maximum of 50 ms between spikes.

Statistics

Significance was determined as p<0.05. Data in bar graphs are means +/- S.E. with n independent observations. Bar averages were compared by ANOVA (unless this is specified above).  To determine whether the peak angles in different subregions were significantly different, the parametric circ_ktest was applied from the MATLAB Circular Statistics Toolbox.