DOI: 10.5061/dryad.7h44j1013
title: The flow of axonal information among hippocampal subregions: 2. Patterned stimulation sharpens routing of information transmission

description: 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.
description: 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<sup>2</sup> 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<sup>2</sup> 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<sub>2</sub> and 9%
 O<sub>2</sub> (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<sub>2</sub>, 9% O<sub>2</sub> (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&gt;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 &gt;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 <em>P</em> is the cumulative probability, ? the
 slope, <em>t</em> the interspike interval in ms and
 <em>c</em> the intercept. The best fits were found through
 performing a grid search to find the local maximum (highest
 R<sup>2</sup>) 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<sup>2</sup> 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<sup>2</sup> 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 (α&lt;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&lt;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.
description: # 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&gt;DG&gt;CA3&gt;CA1&gt;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. ## Code/Software 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.