Skip to main content

The flow of axonal information among hippocampal subregions: 2. Patterned stimulation sharpens routing of information transmission

Cite this dataset

Lassers, Samuel; Vakilna, Yash; Tang, William; Brewer, Gregory (2023). The flow of axonal information among hippocampal subregions: 2. Patterned stimulation sharpens routing of information transmission [Dataset]. Dryad.


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.

README: The Flow of Axonal Information Among Hippocampal Subregions: 2. Subregionally-Specific Network Responses to Patterned Stimulation Sharpens Routing of Information Transmission

Contained are spike data from MultiChannel Systems MEA 120 electrodes in cultured mammalian hippocampi. Directionally specific axonal transmission was measured by the delay between two electrodes spanning a microfluidic tunnel. Tunnels are narrow enough whereas the majority of tunnels contain only one axon. Two stimulation patterns were developed to induce long-term potentiation. High-frequency stimulation (HFS) trains were delivered at high and low theta intervals. The first set of HFS patterns has five cycles delivered at high theta intervals (HFS 5). The second set of HFS patterns has forty cycles delivered at low theta intervals (HFS 40). There are nine unstimulated recordings and six stimulated recordings per stimulation pattern in the same arrays. Axonal tunnels are as follows:

  • A6-A7
  • B6-B7
  • C6-C7
  • D6-D7
  • E6-E7
  • F1-G1
  • F2-G2
  • F3-G3
  • F4-G4
  • F5-G5
  • H6-H7
  • J6-J7
  • K6-K7
  • L6-L7
  • M6-M7
  • F8-G8
  • F9-G9
  • F10-G10
  • F11-G11
  • F12-G12

More information on the data layout can be found here:

Arrays are cultured to include subregions of the hippocampus from the trisynaptic loop. Clockwise arrays have the EC in the upper left quadrant between the tunnels, DG is in the upper right, CA3 is in the lower right, CA1 is in the lower left. Counterclockwise arrays have the EC in the upper left, DG in the lower left, CA3 in the lower right, and CA1 in the upper right.

We found that after stimulation in communicating axons, interspike intervals and inter burst intervals increased while intraspike intervals decreased and burst duration increased. This shows that information is being delivered at a slower rate, but the packets of neural activity contained more information. Spiking and bursting activity balanced with more feedback activity than feed-forward activity. Additionally, graphic analysis shows that when stimulation is applied, the number of functional connections between wells and tunnels decreases but the reliability of their spike rate proportions increases. This is a network computation strategy that we have called “sharpening.”

Description of the data and file structure

Tunnel and well spike data are 300 seconds per recording sampled at 25 kHz. The zip files contain the Wave_Clus times files with unsorted spikes for axons in tunnels. Our method of axonal identification required one cluster, so Wave_Clus was used to eliminate noise in the times file output. Wave_Clus times files are a structure that contains the spike outlines, the cluster class, and the parameters Wave_Clus used for superparamagnetic clustering. The cluster class contains two columns and is the only variable from this structure that is used, the first contains integers denoting what cluster number the spike belongs to. The second column contains the time in ms that the spike occurred. The spike shape is under the structure name “spikes” where the rows are an individual spike, and the columns are voltage readings of each spike. The forced variable of the times structure contains a logical array that denotes which spike was forced into a cluster using slightly wider error margins. The parameters contain the parameters used for superparamagnetic clustering. We have empirically changed the default parameters of Wave_Clus to better accommodate our data. The MATLAB parameters file for Wave_Clus can be found in the uploaded files under set_parameters.

The zip files are labeled with a specific file ID (FID) and a stimulation pattern (no stimulation, HFS 5, and HFS 40) for the unsorted tunnels with Wave_Clus times files. The axon-sorted data are structured in a MATLAB cell array with sequential FIDs for each stimulation type and will be labeled (stimulation pattern)SortedAxons.mat. Each cell within this data structure contains a table with all tunnels and associated spike trains. If multiple axons were found, multiple [1x7500000] logical arrays would be assigned to the same cell in the table. If no axons were found the cell would be empty. Feed-forward and feedback spike trains are stored for each axon for each tunnel electrode. The "upstream" channel is the channel that comes earlier in the trisynaptic loop (EC>DG>CA3>CA1>EC) and corresponds to the first electrode named in the hyphenated electrode pair. For example, the upstream channel in the tunnel "E6-E7" is E6. The other electrode is the downstream channel. This nomenclature persists for both feed-forward or feedback axons, E6 would still be the upstream channel for either a feed-forward or feedback axon. The conduction time for each axon is also identified for each axon and direction is specified in the table cells. The well spike data is stored in a MATLAB cell array and can be found under (stimulation pattern)WellSpikes.mat. Each cell contains a table that corresponds with the sequential FIDs used. Each table contains four columns with the subregion, electrode name, region id (equal to subregion), and the spike train in a [1x7500000] logical array.

Once these structures are created, we calculate all associated spiking dynamics of the file which are stored in MATLAB tables. These are available in this repository under files with "SpikeDynamics" in the title with the associated stimulation condition. Each row in the table corresponds to an electrode of a file ID where the interspike intervals (ISI), interburst intervals (IBI), spikes per burst (SPNB), burst duration (BD), intraburst spike rate (IBSR), and the bounds of the burst are stored in [1xn] cell arrays where the timing is in seconds. In the case of axons, directionality is also specified with a 1 for feed-forward and a 0 for feedback and only the upstream electrodes are considered in a tunnel. This dataset can be used to calculate all complimentary cumulative histograms, probability density functions, and graphical analysis found in our paper.


MATLAB R2022a was used for analysis. Please see the github repository for more information. where our preprocessing pipelines for both tunnel and well analysis are available and the associated analysis files.


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 mm2 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/mm2 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% CO2 and 9% O2 (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% CO2, 9% O2 (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.


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 R2) 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 R2 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 R2 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.


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.