The subregions of the hippocampal formation are essential for episodic learning and memory formation, yet the spike dynamics of each region contributing to this function are poorly understood, in part because of a lack of access to the inter-regional communicating axons. Here we reconstructed hippocampal networks confined to four subcompartments in 2D cultures on a multi-electrode array that monitors individual communicating axons. In our novel device, somal and axonal activity were measured simultaneously with the ability to ascertain the direction and speed of information transmission. Each subregion and inter-regional axons had unique power-law spiking dynamics, indicating differences in computational functions. After stimulation, spiking and burst rates decreased in all subregions, spikes per burst generally decreased, intraburst spike rates increased, and burst duration decreased, which were specific for each subregion. These changes in spiking dynamics post-stimulation were found to occupy a narrow range, consistent with the maintenance of the network at a critical state. Functional connections between the subregion neurons and communicating axons in our device revealed homeostatic network routing strategies post-stimulation in which spontaneous feedback activity was selectively decreased and balanced by decreased feed-forward activity. Post-stimulation, the number of functional connections per array decreased, but the reliability of those connections increased. The networks maintained a balance in spiking and bursting dynamics in response to stimulation and sharpened network routing. These plastic characteristics of the network revealed the dynamic architecture of hippocampal computations in response to stimulation by selective routing on a spatiotemporal scale in single axons.

Four Chamber In Vitro Hippocampal Neuronal Network Culture

Details were provided previously in Vakilna (Vakilna et al., 2021). Briefly, MEA120 glass multielectrode arrays (MEA) with 120 30 μm diameter electrodes spaced 0.2 mm apart were used as substrates for the culture of neuronal networks (Multichannel Systems, Reutlingen, Germany; ALA Scientific, Farmingdale, NY, USA). The MEA was divided by a custom polydimethylsiloxane (PDMS) device into four subregions each of 9.7 mm² by 1 mm high wells and 51 microfluidic tunnels 3 μm high x 10 μm wide x 400 μm long spaced 50 μm apart (Figure 2). Five of the 51 interregional tunnels were monitored with pairs of electrodes to determine the direction of action potentials in single axons. This dedicated 19 electrodes into each subregional well. After oxygen plasma treatment of the substrate, the wells were coated with poly-D-lysine for cell adhesion. Hippocampal subregions were micro-dissected from postnatal day 4 Sprague-Dawley rat pups under anesthesia as approved by the UC Irvine Institutional Animal Care and Use Committee. Brain cells were dissociated and plated at 1,000 cells/mm² for DG (including the hilus), 330 for CA3, 410 for CA1(including subiculum), and 330 for EC, mimicking the ratio of neural densities in vivo: EC-DG 1:3, DG-CA3 3:1, CA3-CA 11:1.25, and CA1-EC 1.25:1 (Braitenberg, 1981). The cells in 10 μL of NbActiv4 medium (Transnetyx BrainBits, Springfield, IL; (Brewer et al., 2009b) were plated into the wells sequentially to allow for adhesion before the dish was filled with medium after 30 min in the incubator. The cultures were capped with Teflon membrane covers (ALA Scientific, Farmingdale NY) and incubated for 21–26 days in humidified 5% CO₂ and 9% O₂ (Brewer & Cotman, 1989). Half of the medium was changed every 7 days. Activity was recorded 2–5 days after a medium change when the networks reached maturity. 

Multielectrode Array Recording and Patterned Stimulation

Spontaneous five-minute recordings and stimulations were produced from the 120-electrode microarray with a Multichannel Systems MEA 120 1100x amplifier (Multichannel Systems, Reutlingen, Germany) at a sampling rate of 25 kHz at 37°C in humidified 5% CO₂, 9% O₂ (custom Airgas USA, Santa Ana, CA). Recordings were initiated several minutes after removal of networks from the culture incubator, shortly after stable activity was seen in 80% of the tunnels. Arrays with less than 80% active tunnels or that had poor growth in one of the subcompartments were rejected for analysis. The Multichannel Systems software MCRack was used to capture and initially analyze the recordings and stimulate the culture. Stimuli were biphasic at 200 μs, each with current amplitudes first at -10 μA and then +10 μA. Two kinds of stimuli were delivered at three sites of spontaneous activity. They were patterned within a train containing high-frequency pulses spaced at 10 ms, repeated 5 times (100 Hz), and trains spaced 50 ms apart with four train repeats (termed 5 HFS) (Figure 2). The second pattern was also 100 Hz pulses spaced at 10 ms, but repeated 40 times, and trains spaced 200 ms apart with five train repeats (termed 40 HFS). Stimuli were applied at three sites in each subregion in the array with the next train starting as the previous train ended. Spontaneous recordings 67 +/- 20 min (SD) after stimulation are reported from nine separate cultures (9 arrays) listed in Supplemental Table 1; we collected the effects of stimulation for six of these. Clockwise direction of plating: (CW) was EC-DG-CA3-CA1. Counterclockwise (CCW) was CA1-CA3-DG-EC to control for plating order bias. None was detected.

Spike Detection, Sorting, and Axonal Propagation Direction Pipeline

We made automation improvements to our axonal spike directionality algorithm (Vakilna et al., 2021). Raw tunnel data sampled at 25 kHz was filtered through Wave_Clus (Chaure et al., 2018) which uses a second-order elliptic bandpass filter at 300 Hz to 10,000 Hz for spike detection and a fourth-order elliptic bandpass filter at 300 Hz to 3,000 Hz for clustering. Spike detection and clustering was computed above two different thresholds: 5 to 50 S.D. noise and 50.1 to 500 S.D. to ensure that large spikes in the axons were accurately counted and clustered. Through Wave_Clus, the negative peak of spikes was detected with a minimum of 1 ms before the peak and 2 ms after the peak. A refractory period of 1.5 ms was specified. Spike shapes differing by at most three standard deviations from the mean spike shape were included in a single cluster. Clustering was primarily used to discard complex spikes from two or more axons in a microtunnel causing overlapping spike waveforms. All spikes collected from a single electrode in high and low clustering were consolidated into a single cluster. Since previous analyses of these tunnel devices showed that ~63% of tunnels only had one axon (Narula et al., 2017), we could identify single axons by their uniform conduction velocities, i.e., by spike timing delays, which distinguished multiple axons per tunnel. We empirically checked waveforms after “clustering” through timing delays to ensure singular waveforms. Axonal conduction time delay was used to generate a normalized matching indexing (NMI). With two electrodes spanning the same distance in each tunnel, a timing comparison was made, and an NMI algorithm was computed for every tunnel. Tunnels were considered valid for NMI>0.2 (20% of spikes matched between two clusters) which was sufficient for eliminating spurious spike pair correlations during high spike rates. Paired spikes were identified as spikes with a delay between 0.2 and 1 ms over the 0.2 mm distance, corresponding to the physiological bounds of axonal propagation velocities of hippocampal action potentials (0.4–2 mm/s) (Colombe & Ulinski, 1999). A histogram of conduction times was created with empirically determined thresholding and peak prominence values using the MATLAB findpeaks function. The threshold used was 1.9 standard deviations and the peak prominence parameter was at least 12% of the highest bin. Peaks in the histograms were identified as belonging to different axons +/-0.16 ms, validated with spike shapes. Delay times up to 0.12 ms were considered if there was a peak at sufficiently fast conduction times. Positive delay times indicate feed-forward axons while negative conduction times indicate feedback axons. Example histograms of four types of measured spike timing delays for a given tunnel, and the associated waveforms are provided in Supplementary Figure 1. This is the basis for our determination of feed-forward and feedback directionality of axonal communication. In contrast to the tunnels, in the wells, raw data was filtered through Wave_Clus and spikes were detected at +/- 5 S.D. noise. In the wells, neuronal somata density and electrode placement was sparse enough to detect a single neuron in >90% of cases, so spike clustering was not needed.

Measured Spike Dynamics and Probability Distributions

As in Vakilna et al. (2021), the distribution of inter-spike intervals (ISI), inter-burst intervals (IBI), and spikes per burst (SPB) follow log-log distributions and were visualized as normalized complementary cumulative probability distributions (CCDs) with logarithmically spaced bins (Newman 2005). A log-transformed linear model was used to fit the CCD after log transformation.

 ???10(?)=?×???10(?)+?

Where P is the cumulative probability, ? the slope, t the interspike interval in ms and c the intercept. The best fits were found through performing a grid search to find the local maximum (highest R²) with time limits varied up to 50% with a step size of 5%. For ISI, a single fit was found for all CCDs considering probabilities from 1 to 0.1 and intervals from 0.01 to 1.0 s except EC-DG feed-forward axons which were fitted for intervals between 0.1 to 0.2 seconds to account for non-linearities in the distribution. For inter-burst intervals, two linear fits were calculated for all CCDs piecewise to account for the “up states” and “down states” which refer to fast and slow bursting, respectively (Vakilna et al., 2021). The minimum time of the up states was used as the maximum time for the down states. The distributions of intra-burst spike rates and burst durations better followed a log-normal distribution and were visualized with a probability distribution (Vakilna et al., 2021). The median was calculated by fitting a normal distribution to the probability distribution in natural log space. The mean in natural log space (m) was extracted from the fit. A one-way ANOVA was calculated from the fits.

Edge Plot Graph Analysis

In graph analysis, edges are the connections between nodes (neurons and tunnel axons). Edge detection was used to correlate firing rates between well and tunnel electrodes. NMI was designed to be very exact and discriminatory for spike identification and required robust signals from both electrodes in a tunnel. However, examining raster plots, some channels exhibited signs of lower axonal coupling to the electrode. Therefore, the more active electrode was used as the basis for edge detection. The average firing rate for each 100 ms bin was calculated from all combinations of each axon in a tunnel and each adjacent well electrode with a 10 s wide sliding window using 100 ms steps. If the window contained at least four bins with a non-zero firing rate in a combination of well and tunnel, a linear regression correlating the wells and tunnels was calculated. This was repeated for all possible connections. For plotting, all linear regressions with a positive slope greater than 0.1 and an R² greater than 0.2 were used. A slope less than 1 indicated a decreased firing rate through that connection and a slope greater than 1 indicated an increase in firing rate. The reliability of this connection was described through the R² values (Pearson linear regression). These values were averaged over all the time windows for each array to derive weights and graphically depicted to generalize the active connections in the network. MATLAB graph analysis functions were employed to calculate the number of edges (connections from a well node to a tunnel axon node), degrees (number of edges per node), and centrality (degrees weighted for slope or R2) of each node for feed-forward and feedback connections for well-to-tunnel edges.

Statistics

Data were analyzed with custom MATLAB 2020a scripts. Slopes were compared for a statistically significant difference (α<0.05) using analysis of covariance (ANCOVA) and Tukey’s honestly significant difference (HSD) test. The significance within a stimulation set between subregions and between stimulation sets of the same subregion were analyzed by analysis of variance (ANOVA) and Tukey’s HSD test with the null hypothesis rejected for p<0.05. Data were analyzed for nine separately plated networks, six of which were subjected to patterned stimulation. Details are summarized in Supplemental Table 1.