Distinguishing cognitive state with multifractal complexity of hippocampal interspike interval sequences
1Neuroscience Program, Wake Forest School of MedicineWinston-Salem, NC, USA
2Department of Physiology and Pharmacology, Wake Forest School of MedicineWinston-Salem, NC, USA
3Department of Biomedical Engineering, Wake Forest School of MedicineWinston-Salem, NC, USA
4Department of Biomedical Engineering, University of Southern CaliforniaLos Angeles, CA, USA
*Correspondence: Dustin Fetterhoff, Department of Physiology and Pharmacology, Wake Forest School of Medicine, Medical Center Blvd., Winston-Salem, NC 27157, USA dfetterh@wakehealth.eduAbstract
Fractality, represented as self-similar repeating patterns, is ubiquitous in nature and the brain. Dynamic patterns of hippocampal spike trains are known to exhibit multifractal properties during working memory processing; however, it is unclear whether the multifractal properties inherent to hippocampal spike trains reflect active cognitive processing. To examine this possibility, hippocampal neuronal ensembles were recorded from rats before, during and after a spatial working memory task following administration of tetrahydrocannabinol (THC), a memory-impairing component of cannabis. Multifractal detrended fluctuation analysis was performed on hippocampal interspike interval sequences to determine characteristics of monofractal long-range temporal correlations (LRTCs), quantified by the Hurst exponent, and the degree/magnitude of multifractal complexity, quantified by the width of the singularity spectrum. Our results demonstrate that multifractal firing patterns of hippocampal spike trains are a marker of functional memory processing, as they are more complex during the working memory task and significantly reduced following administration of memory impairing THC doses. Conversely, LRTCs are largest during resting state recordings, therefore reflecting different information compared to multifractality. In order to deepen conceptual understanding of multifractal complexity and LRTCs, these measures were compared to classical methods using hippocampal frequency content and firing variability measures. These results showed that LRTCs, multifractality, and theta rhythm represent independent processes, while delta rhythm correlated with multifractality. Taken together, these results provide a novel perspective on memory function by demonstrating that the multifractal nature of spike trains reflects hippocampal microcircuit activity that can be used to detect and quantify cognitive, physiological, and pathological states.
Introduction
By analyzing the mono- and multifractal properties of neural temporal dynamics, we may generate new insights concerning how the brain functions with implications for detection of cognitive, physiological, and pathological states. Such analyses have been used successfully to detect pathological conditions such as heart disease (Ivanov et al., 1999), Alzheimer's disease (Lahmiri and Boukadoum, 2013), Parkinson's disease (Zheng et al., 2005), and epilepsy (Serletis et al., 2012; Dutta et al., 2014). Additionally, multifractal analysis detects clear differences in neural activity between wakefulness and sleep stages using EEG signals (Weiss et al., 2009; Zorick and Mandelkern, 2013). Multifractal complexity of time series is believed to indicate functional connectivity because such complexity would provide the essential temporal dynamics to support information transfer via variable patterns of neural activity. Interactions across brain regions, detected as multifractal complexity, regularly fluctuate between task and rest conditions specifically in regions associated with the task (Ciuciu et al., 2012). In order to assess cognitive state detection abilities, a paradigm was implemented to examine how multifractal complexity is reflected by active (i.e., task-related) hippocampal microcircuit processing.
Temporal coding analyses attempt to derive information about brain function from the timing of action potentials generated by neuronal ensembles or from rhythmic neuronal oscillations, such as theta rhythm (Jones and Wilson, 2005). Some temporal coding analyses assign physiological function to activity within frequency bands (i.e., theta-phase precession or cross frequency coupling), but we propose that multifractal analysis can provide new insights into mechanisms of neurological temporal coding because it quantifies the structure of variability and the self-similar (fractal) nature of physiological systems. Scale-free dynamics are often associated with long-range temporal correlations (Linkenkaer-Hansen et al., 2001; Ciuciu et al., 2012) and quantified by the monofractal Hurst exponent. However, a single exponent does not capture the complexity of many physiological signals, supporting the use of a spectrum of scale invariant exponents that describe the multiple, co-existing fractal patterns (Dixon et al., 2012). Memory is commonly believed to occur through repetitive neuronal sequences (Hampson et al., 2012), and therefore multifractal analysis applied to spike train patterns may quantify a possible basis of memory detected as multifractal complexity.
To assess if active memory processing is reflected by multifractal measures, analyses of monofractal long-range temporal correlations (LRTCs) and multifractal complexity in hippocampal interspike interval (ISI) neural sequences were conducted during a working memory task. In vivo electrophysiological recordings of rat hippocampal CA3 and CA1 subregions were conducted during a resting state condition in a neutral (task-independent) context for 25–30 min both before and after performance of the delayed nonmatch-to-sample (DNMS) task. Between the pre-task recording and the DNMS task, rats were injected with vehicle or tetrahydrocannabinol (THC), a psychoactive component of cannabis known to impair memory function (Hampson and Deadwyler, 2000). Prior results demonstrated that hippocampal neurons with memory-correlated firing rate alterations (functional cell types, FCTs) recorded during the DNMS task were more multifractal than non-memory neurons (non-FCTs) and THC administration impaired memory while reducing multifractality (Fetterhoff et al., 2015). By examining the same neurons before, during and after the DNMS task, alterations in multifractality were assessed in a different context. These experiments and analyses were designed to extend previous findings by testing three new hypotheses and facilitating a stronger intuition concerning multifractal properties of hippocampal microcircuits. First, we hypothesized that LRTCs, as indicated by the Hurst exponent, would decrease during the DNMS task compared to resting (pre/post) recording conditions. Since LRTCs arise when distant activity has a greater influence on future activity patterns, we hypothesized a decrease would occur due to the constantly changing requirements of the DNMS task. Second, we hypothesized that an increase in multifractal complexity reflects active memory processing in populations of hippocampal neurons and therefore, spike trains should be more multifractal during the DNMS task compared to both pre- and post-task recording phases. Third, we hypothesized that THC would decrease both multifractal complexity and Hurst exponents during the task and post-recording phases compared to vehicle control recordings during the same phase. Finally, to enrich conceptual interpretation of multifractal complexity and LRTCs and establish the difference between structure (multifractality) and amount of variability, classical spike train variability measures (coefficient of variation, ISI STD and mean ISI) were compared with the mono- and multifractal variables.
The primary goal of this study was to assess the capacity of multifractal analysis to distinguish between recording phases and drug conditions. Fourier analysis of neuronal signals is one commonly employed method to distinguish between physiological and cognitive states (De Carli et al., 2004; Jones and Wilson, 2005; Nguyen et al., 2008; Palva et al., 2010; Van Someren et al., 2011; Garn et al., 2014), and therefore, the performance of multifractal analysis was compared to the frequency content computed from the same spike trains. The results showed that the monofractal Hurst exponent and magnitude of multifractality could differentiate between more recording/drug conditions compared to frequency content (theta and delta) and further support utility of multifractal analysis for this objective (Weiss et al., 2009; Zorick and Mandelkern, 2013). Multifractal analysis has the potential to generate novel insights into the role of neuronal ensembles by quantifying different temporal features compared to other analyses.
Materials and methods
Rats
Male Long-Evans rats (Harlan) aged 6–10 months (n = 10) were tested under protocols approved by the Wake Forest University Institutional Animal Care and Use Committee, and in accordance with the Association for Assessment and Accreditation of Laboratory Animal Care and the National Institute of Health Guide for the Care and Use of Laboratory Animals (NIH Publication No. 8023). All animals were individually housed and allowed free access to food with water regulation to maintain 85% of ad libitum body weight during testing. Upon termination of the study, all rats were anesthesized with ketamine (100 mg/kg) and brains were perfused with formaldehyde for preservation and subsequent histology to confirm electrode placement.
Apparatus
The behavioral testing apparatus for the delayed nonmatch to sample (DNMS) task as used in other studies (Hampson and Deadwyler, 2000; Hampson et al., 2012) consisted of a 43 × 43 × 50 cm Plexiglas chamber with two retractable levers (left and right) positioned on either side of a water trough on the front panel. A nosepoke device (photocell) was mounted in the center of the wall opposite the levers with a cue light positioned immediately above the nosepoke device. A video camera was mounted on the ceiling and the entire chamber was housed inside a commercially built sound-attenuated cubicle.
Delayed nonmatch-to-sample (DNMS) task
The DNMS task consisted of three main phases: Sample, Delay and Nonmatch (Figure 1B). The Sample phase initiated the trial via presentation of either the left or right lever (50% probability), which required the animal to press and make the Sample Response (SR). The lever was then retracted and the Delay phase of the task initiated, as signaled by the illumination of a cue light over a nosepoke photocell device on the wall opposite to where the lever was presented. At least one nosepoke (NP) was required following the interposed delay interval which varied randomly in duration (1–30 s) on each trial during the session. The Nonmatch phase began when the delay timed out, the photocell cue light turned off, and both the left and right levers on the front panel were extended. Correct responses consisted of pressing the lever in the Nonmatch phase located in the spatial position opposite to the position of the SR; in other words, a Nonmatch response (NR). This produced delivery of a 0.4 ml water reward in a reservoir between the two levers. After the NR the levers were retracted for a 10.0 s intertrial interval (ITI) before the Sample lever for the next trial was presented. A lever press at the same position as the SR (match response) constituted an “error” with no water delivery and turned off chamber house lights for 5.0 s with the next trial presented 5.0 s later. Individual performance was assessed as % NRs (correct responses) with respect to the total number of trials (80–100) per daily (1 h) session.
Task and rest (non-task) recording paradigm
All rats were recorded for 25–30 min in a bare, white 38 × 29 × 30 cm plastic container that was inserted within the DNMS testing chamber both before (pre) and after (post) the DNMS task (Figure 1A). This constituted a neutral environment used to record task-independent neuronal activity that was opaque to prevent rats from seeing task components, such as levers or the nosepoke device.
Drug preparation and administration
Δ9-THC was obtained from the National Institute on Drug Abuse as a 50 mg/ml solution in ethanol. Detergent vehicle was prepared from Pluronic F68 (Sigma, St. Louis, MO), 20 mg/ml in ethanol. Δ9-THC was added to the detergent-ethanol solution (0.04–0.12 ml of THC), and then 2.0 ml of saline (0.9%) was slowly added to the ethanol-drug solution. The solution was stirred rapidly and placed under a steady stream of nitrogen gas to evaporate the ethanol (~10 min). This resulted in a detergent-drug suspension (12.5 mg/ml THC), which was sonicated and then diluted with saline to final injection concentrations (1.0–3.0 mg/ml THC). On testing days, animals were injected intraperitoneally with the detergent vehicle solution or the THC drug-detergent solution (1 ml/kg) immediately after the pre-task recording phase and approximately 5–10 min before the start of the behavioral session. At least one vehicle day was imposed between each drug-testing day. All rats received THC on 5–8 days. All drug solutions were mixed fresh each day.
Hippocampal electrode array surgery
All surgical procedures conformed to National Institutes of Health and Association for Assessment and Accreditation of Laboratory Animal Care guidelines, and were performed in a rodent surgical facility approved by the Wake Forest University Institutional Animal Care and Use Committee. Electrode arrays and recordings were the same as described in several prior publications from this laboratory (Hampson et al., 1999, 2012; Hampson and Deadwyler, 2000). After being trained to criterion performance level in the DNMS task animals were anesthetized with ketamine (100 mg/kg) and xylazine (10 mg/kg) and placed in a stereotaxic frame. Craniotomies (5 mm-diameter) were performed bilaterally over the dorsal hippocampus to provide for implantation of 2 identical array electrodes (Neurolinc, New York, NY), each consisting of two rows of 8 stainless steel wires (diameter: 20 μm) positioned such that the geometric center of each electrode array was centered at co-ordinates 3.4 mm posterior to Bregma and 3.0 mm lateral (right or left) to midline (Paxinos and Watson, 1997). The array was designed such that the distance between two adjacent electrodes within a row was 200 μm and between rows was 400 μm to conform to the locations of the respective CA3 and CA1 cell layers. The longitudinal axis of the array of electrodes was angled 30° to the midline during implantation to conform to the orientation of the longitudinal axis of the hippocampus, with posterior electrode sites more lateral than anterior sites. The electrode array was lowered in 25–100 μm steps to a depth of 3.0–4.0 mm from the cortical surface for the longer electrodes positioned in the CA3 cell layer, leaving the shorter CA1 electrodes 1.2 mm higher with tips in the CA1 layer. After placement of the array the cranium was sealed with bone wax and dental cement and the animals treated with Ketoprofen (3.0–5.0 mg/kg) for pain relief over the next 4–6 h. The scalp wound was treated periodically with Neosporin antibiotic and systemic injection of penicillin G (300,000 U, intramuscular) were given to prevent infection. Animals were allowed to recover from surgery for at least 1 week before continuing behavioral testing (Berger et al., 2011; Hampson et al., 2012).
Electrophysiological monitoring, acquisition and waveform sorting of neuronal data
Animals were connected by cable to the recording apparatus via a 32-channel headstage and harness attached to a 40-channel slip-ring commutator (Crist Instruments, Hagerstown, MD) to allow free movement in the behavioral testing chamber. Single neuron action potentials (spikes) were isolated by time-amplitude window discrimination and computer-identified individual waveform characteristics using a multi-neuron acquisition (MAP) processor (Plexon Inc., Dallas, TX, USA). Single neuron spikes were recorded daily and identified using waveform and firing characteristics within the task (perievent histograms) for each of the DNMS events (SR, LNP, and NR). To maintain waveform shape across days, all recorded data was concatenated into one file (separately for each rat) and offline sorting was performed using principal component analysis, peak-valley, and nonlinear energy algorithms in Offline Sorter (Figure S1; Plexon Inc., Dallas, TX, USA). Hippocampal neuron ensembles used to distinguish recording phases and drug treatment conditions consisted of 10–30 single neurons, each recorded from a separate identified electrode location on either of the bilateral arrays. Only isolated neural spike waveforms exhibiting firing rates of CA1 and CA3 principal cells (i.e., 0.5–8.0 Hz average firing rate) and consistent multifractal properties across sessions were included in analyses. Previous work has shown that hippocampal neurons recorded identified in this manner exhibit consistent mean, baseline, and DNMS task modulated firing rate over multiple task sessions (Deadwyler et al., 1996; Hampson et al., 1999, 2003).
Multifractal detrended fluctuation analysis (MFDFA)
Detailed descriptions of multifractal detrended fluctuation analysis (Kantelhardt et al., 2002; Kantelhardt, 2012) and associated Matlab code (Ihlen, 2012) are available elsewhere. We briefly summarize and illustrate the main components of the MFDFA method. A demonstration of this method (Figure 2) was constructed using recordings of a single neuron (Figure S1) recorded during the DNMS task over two sessions/days: one vehicle (blue) and one THC (green). First, a neuronal spike train is converted to a sequence of interspike intervals (ISIs; Figure 2A), represented as x, and commonly referred to as a noise-like time series (Ihlen, 2012). The ISIs are converted into a “random walk-like” time series Y(i) by subtracting the mean and integrating the ISI signal x:
The random walk Y(i) (Figure 2B inset, blue line) is divided into Ns non-overlapping segments of equal length s (scale). The local quadratic trend yv (Figure 2B inset, black line) is calculated for each of the Ns segments v by a least-square fit of the series to determine the local root mean square variation F for each segment v (Figure 2B):
In this way, F2(v, s) is essentially the mean-square error difference between the quadratic fit yv and the walk-like ISI sequence Y. Visualization of the walk like ISI sequence at multiple q-order statistical moments permits observation of the structure of variability quantified by multifractal analysis (Figure 2B). Periods of low activity (i.e., clusters of short ISIs, faster firing rate) are amplified with negative moments (Figure 2B). Conversely, periods of high activity (i.e., clusters of longer ISIs, slower firing rate) are amplified with positive moments (Figure 2B, bottom). Next, the qth order fluctuation function is determined by averaging over all segments v for each qth power (Figure 2C):
The scaling behavior of the fluctuation function is seen by analyzing log-log plots of Fq(s) vs. s for each qth power (Figure 2C). Standard (monofractal) Detrended Fluctuation Analysis calculates the Hurst exponent from the slope of the power-law regression line between the overall root mean square variation F across multiple scales, s, for a single statistical moment, q = 2. Multifractal analysis performs the same linear regression for a broad range of statistical moments q. Fq(s) increases as a power-law for large values of s if the signal x contains long-range temporal correlations (LRTCs):
The Hurst exponent (Hurst, 1951) is the slope of Fq(s) on this log-log plot for H(q = 2). A greater slope yields a larger Hurst exponent and indicates increased LRTCs which defines how fast the overall root mean square variation grows with increasing segment size (scale = s). Hurst exponents greater than 0.5 indicate the time series contains positive correlations (i.e., persistent structure), Hurst exponents ranging from 0 to 0.5 indicate negative correlations and uncorrelated Gaussian noise has a Hurst exponent equal to 0.5. In multifractal signals, changes in LRTCs occur at different q-order statistical moments. These variations are visualized by comparing the q-order Hurst exponents (slopes of regression lines) on log-log plots of Fq(s) vs. scale s, where each line is computed from a different (integer) qth power ranging from −3 to 3 (Figure 2C). Our example shows more variable slopes in the vehicle condition (Figure 2C, blue) compared to the similar slopes in THC condition (Figure 2C, green), indicating greater multifractality during the vehicle condition. The qth order fluctuation function Fq(s) is one way to visualize multifractal properties of variability, but generally the multifractal singularity spectrum (Figure 2D) is constructed to illustrate the distinction between (mono)-fractal and multifractal analyses. The q-order Hurst exponent H(q) is converted to the q-order mass exponent τ(q):
Finally, a Legendre transform relates τ(q) to the fractal dimension D(h) and Hölder exponent h:
The singularity spectrum is like a histogram of the q-order Hurst exponents (slopes of regression lines in Figure 2C). The magnitude of multifractality is determined by the width of the singularity spectrum and consequently the range of Hölder exponents h covered by the ISI signal (Figure 2D). The Hurst exponent is closely approximated by the Hölder exponent at the apex of the singularity spectrum (where D(h) = 1). In our example, the singularity spectrum for the THC condition is narrower, and thus less multifractal than the singularity spectrum computed for the control condition (Figure 2D). The parameters used for all analyses were determined by viewing log-log plots of Fq(s) vs. s for each qth power and multifractal singularity spectra that resembled those found in the literature (Kantelhardt et al., 2002; Ihlen, 2012). MFDFA was performed by fitting a second-order polynomial, scales used ranged from 16 to 256, and singularity spectrum width was determined by computing h(q = 3)−h(q = −3). All analyses were performed in Matlab using publicly available MFDFA code (Ihlen, 2012).
Fourier transform
A fast Fourier transform of the spike train data in 1 ms binary bins where 1 = spike and 0 = no spike was performed for every neuron recorded during every recording phase. Delta and theta power were measured by taking the ratio of signal power in the delta (0.5–4.0 Hz) or theta range (4–8 Hz) over the total power of the normalized signal, bandpass filtered from 0.5 to 12 Hz. All analyses were performed using Matlab.
Statistical analyses
All calculations were performed in Matlab, results were recorded spreadsheets and imported into Statistical Analysis Systems (SAS) software (SAS Institue, Cary, NC) to perform repeated measures ANOVAs and correlations. A total of seven repeated measures ANOVAs were performed using these dependent variables: the Hurst exponent, singularity spectrum width, coefficient of variation, ISI standard deviation, mean ISI, delta power, and theta power. Neurons were the subject identifier, each session as a within-subjects effect, and no group identifier was used. The covariance structure used was compound symmetry. Main effects of drug condition or recording phase were only discussed when the interaction between drug condition and recording phase was non-significant. When significant ANOVA effects were found, post-hoc tests were performed with 10,000 Monte Carlo permutations to determine two-tailed p-values. Significance levels were set to p < 0.0083 to control for multiple comparisons (Bonferoni correction). Correlations (Spearman's rho) were used to examine the associations between mono- and multifractal variables (Hurst exponent and singularity spectrum width) and coefficient of variation, ISI standard deviation, mean ISI, delta and theta power. 55 total correlations were performed to assess the relationships for monofractal and multifractal variables during all three recording phases (pre-task, task, post-task).
Results
To investigate dynamical interspike interval (ISI) patterns associated with different hippocampal microcircuit processing states, a comparision was made between those generated during the DNMS task (Deadwyler et al., 1996; Hampson et al., 1999) and those occuring during a resting state after either vehicle or THC administration (Hampson and Deadwyler, 2000). Hippocampal principal (pyramidal) cells were identified based on mean firing rate (0.5–8.0 Hz), and neurons exhibiting rates outside this range were excluded. Identified neurons from selected wires were “tracked” from day to day by waveform and multifractal properties (i.e., Hurst exponent and width). A total of 197 hippocampal principal (pyramidal) cells were recorded from 10 different rats and each neuron was recorded from the same electrode over multiple days (10–16 recording days per rat). 117 CA1 and 80 CA3 neurons were analyzed for this study. Every neuron included in the analyses was recorded during at least 4 days (2 vehicle and 2 THC).
Delayed nonmatch-to-sample task
Hippocampal spike trains were recorded during a resting state condition in a neutral (i.e., task-independent) environment both before (pre) and after (post) the DNMS task (Figure 1A) to assess the influence that active memory processing (during the task) exerts on the structure of spike train variability, as indicated by multifractal analysis. This approach was designed to assess electrophysiologcal distinctions between three different recording phases (pre-task, task, post-task) and in two drug conditions. Since drug injections (pluronic vehicle or THC) were given immediately after the pre-task/pre-drug resting state recording phase, all computed measures for the pre-task phase are equal across the two drug conditions. THC, the main active ingredient in cannabis (Gaoni and Mechoulam, 1964), was chosen because it impairs memory encoding during the DNMS task (Hampson and Deadwyler, 1999, 2000), reduces LRTCs and multifractal complexity of task-related neuronal spike trains (Fetterhoff et al., 2015) and impairs theta frequency-related working memory performance in both rats (Robbe et al., 2006) and humans (Ilan et al., 2004; Böcker et al., 2010). THC doses were chosen to maximally impair DNMS performance in order to examine effects on associated multifractal spike train characteristics using previously established dose-response relationships (Hampson and Deadwyler, 2000). Working memory was assessed in 10 rats after vehicle or THC administration using the DNMS task (Figure 1B). A within subjects design was used to assess behavioral performance and hippocampal electrophysiology for 5–8 days per drug condition (vehicle or THC) per rat. The DNMS performance was inversely correlated with delay length, as all animals performed worse at longer delays (Figure 3). THC (green line) impaired performance compared to vehicle (Figure 3, blue line).
Frequency and singularity spectra of hippocampal spike trains
Interactions between frequency bands are believed to represent coordination and exact timing relationships within neuronal ensembles (Buzsáki and Moser, 2013). Multifractal dynamics are believed to arise from interactions across spatiotemporal scales (Ihlen and Vereijken, 2010; Kelty-Stephen et al., 2013) and therefore may represent phenomena related to frequency-specific activity. Additionally, both analyses are commonly used to assess differences in sleep/arousal state, and increasing interest is focused on multifractal analysis as an automatic sleep stage detection method (Weiss et al., 2009; Zorick and Mandelkern, 2013). The performance of multifractal and Fourier analyses was evaluated here in terms of their ability to distinguish between recording phases and drug treatment conditions (Figure 1) based on analyses of hippocampal spike trains. Fourier transforms were used to measure both delta (0.5–4 Hz) and theta (4–8 Hz) power in binary spike train representations. These frequency bands were chosen due to their prominence during the DNMS task (Hyman et al., 2010, 2011) and pervasive presence in working memory literature (Sato and Yamaguchi, 2003; Mormann et al., 2008; Assenza et al., 2013; Clemens et al., 2013; Hasselmo and Stern, 2014). MFDFA quantifies the structure of variability of an ISI sequence and permits estimation of the multifractal singularity spectrum (Ihlen, 2012) that describes a given spike train. The singularity spectrum presents two important pieces of information: the Hurst exponent and magnitude of multifractality. The Hurst exponent H is a self-similarity parameter that quantifies long-range temporal correlations (LRTCs) and takes values nearly equivalent to the Hölder exponent h-value at the center apex of the singularity spectrum (Figure 2D) where D(h) = 1 (Ciuciu et al., 2012; Ihlen, 2012). The Hurst exponent quantifies the (mono)-fractal structure of the ISI time series, however, many physiological signals exhibt a wider range of dynamical activity that is better described as multifractal. The magnitude of multifractality quantifies the structure of variability in an interspike interval (ISI) sequence and is directly proportional to the width of the singularity spectrum. The width is defined by the range of local Hölder exponents h (Figure 2D) which quantify the dynamical profile of a time signal.
The singularity and frequency spectra for two different example neurons are shown in Figure 4 while the population singularity spectra are shown in Figure 5. Repeated measures ANOVA results from population analyses are briefly mentioned here and presented fully in the subsequent sections. The first example permits comparision of all three recording phases taken from vehicle treatment conditions (Figures 4A,B, 5B). One example neuron exhibits increased multifractal complexity during the DNMS compared to either resting state recordings (Figure 4A). The frequency spectra for this same neuron exhibits both delta and theta power in all recording phases (Figure 4B). This neuron illustrates the same effect found in the population (Figure 5B): multifractal complexity (width) increases from post-task to pre-task to task (Figure 7F) and LRTCs (Hurst exponent) are larger during the resting states (pre- and post-task) compared to the task (Figure 6F). Although the singularity spectra are discernable across task phases for this neuron, the frequency spectra were not (Figure 4B). However, the population analyses revealed increased theta power during vehicle resting state recordings (pre- and post-task) compared to vehicle task recordings (Figure 8E).
Another example neuron was chosen to illustrate the finding that THC reduced multifractal complexity (width) during the DNMS task compared to vehicle (Figure 4C). This neuron was recorded over 16 total days, 8 vehicle days and 8 THC days. Individual and average singularity spectra are shown (Figure 4C). The frequency spectrum for this neuron shows delta rhythm but not theta (Figure 4D). This neuron is representative of the entire population: during the DNMS task, THC reduces multifractality (width; Figures 5C, 7D,F), but has no significant effect on the Hurst exponent (Figures 6D,F). Additionally, THC did not affect frequency content during the DNMS task compared to vehicle (Figures 8D,E).
Temporal correlations and multifractal complexity of hippocampal population
Two repeated measures ANOVAs were performed to assess the capacity of monofractal and multifractal variables, the Hurst exponent and singularity spectrum width, to distinguish between recording phases and drug treatment conditions. Hurst and width were used as dependent variables to examine main effects of drug condition (vehicle vs. THC) and recording phase (pre-task/pre-drug, DNMS task, post-task) and an interaction between the two. Only interactions (not main effects) are discussed when significant. We primarily wanted to evaluate the ability of computed variables to establish 6 ad hoc chosen differences: vehicle pre-task vs. task, vehicle task vs. post-task, vehicle pre-task vs. post-task, vehicle task vs. THC task, vehicle post-task vs. THC post-task, and THC task vs. THC post-task. Since the distributions were long-tailed, post-hoc Monte Carlo permutation tests were performed with signifiance levels set at p < 0.0083 (Bonferoni correction). All analyses were performed on a total of 5143 individual interspike interval (ISI) sequences: 829 (495 CA1, 334 CA3) individual ISI sequences from vehicle pre-task recordings, 1004 (592 CA1, 412 CA3) from vehicle task recordings, 848 (510 CA1, 338 CA3) from vehicle post-task recordings, 771 (443 CA1, 328 CA3) from THC pre-task/pre-drug recordings, 923 (550 CA1, 373 CA3) from THC task recordings, and 768 (453 CA1, 315 CA3) from THC post-task recordings. Population statistics were performed for all hippocampal neurons (Figures 6, 7) and differences by hippocampal location were shown in the Supplementary Material (Figure S2).
The Hurst exponent is a monofractal self-similarity parameter that quantifies LRTCs in an ISI sequence. A significant main effect of drug condition [F(1, 196) = 97.61, p < 0.0001], a significant main effect of recording phase [F(2, 379) = 182.7, p < 0.0001], and a significant interaction [F(2, 369) = 80.77, p < 0.0001] were found when using the Hurst exponent as the dependent variable. The significant drug treatment by recording phase interaction revealed three out of six important group differences (Figure 6F): both pre- and post-task vehicle recordings contain greater LRTCs than task recordings (Figures 6A,B) and THC reduces LRTCs compared to control during the post-task recording (Figure 6E). No significant difference was found between the Hurst exponent from vehicle pre-task recordings compared to vehicle post-task recordings (Figure 6C), nor was a difference established between the Hurst exponent computed from vehicle or THC task recordings (Figure 6D). After THC treatment, the Hurst exponent of the task condition was similar to those from the post-task (Figure 6F).
Multifractal complexity reflects energy flow through all scales of a dynamical system (Dixon et al., 2012; Kelty-Stephen et al., 2013) and is greater in memory processing neurons compared to randomly spiking ones (Fetterhoff et al., 2015). Therefore, multifractality may selectively arise in hippocampal microcircuit processing when task-specific input-output transformations occur (Berger et al., 2012; Hampson et al., 2012). A repeated measures ANOVA using multifractal width as the dependent variable yielded a significant main effect of drug condition [F(1, 196) = 11.09, p = 0.001], a significant main effect of recording phase [F(2, 379) = 30.8, p < 0.0001], and a significant interaction [F(2, 369) = 13.87, p < 0.0001]. The significant drug condition by recording phase interaction established four out of six significant differences (Figure 7F): in the vehicle condition, all three recording phases are different from each other (Figures 7A–C; task > pre > post). Multifractal complexity was greatest when hippocampal ensembles were processing task-relevant information. THC reduced multifractal complexity (width) during the task (Figure 7D) but had no effect on multifractal complexity during post-task recordings (Figure 7E). Therefore, the effect of THC to reduce multifractality of hippocampal neurons only occurred when memory was impaired during the task. Additionally, multifractal complexity (width) of neurons after THC administration was similar during DNMS task and post-task phases (Figure 7F). The lack of multifractality adjustments between task and post-task recordings after THC administration suggests that THC impairs functional transistions in hippocampal microcircuit activity detected as multifractal complexity.
Standard variability measures are marginally task-specific
Multifractal analysis quantifies the structure of variability (Ihlen, 2012); therefore, to establish the differences between structure and amount of variability, we compared multifractal indices to standard variability measures of the same hippocampal multi-neuron data. To determine if the amount of variability of hippocampal spike trains can account for distinctions between recording phases and drug conditions, the coefficient of variation (CV) was computed for all hippocampal ISI sequences by dividing ISI standard deviation by mean ISI. A CV greater than 1 is an indicator of neuronal bursting and a CV equal to 1 indicates a poisson process. Repeated measures ANOVAs were performed for all three variables: CV, ISI STD and mean ISI. For CV, there was no significant main effect of THC administration [F(1, 196) = 0.68, p = 0.4096], nor was there a significant interaction between drug condition and recording phase [F(2, 369) = 1.61, p = 0.2008]. A significant main effect of recording phase was found [F(2, 379) = 3.76, p = 0.0241] and post-hoc assessment showed that CVs were higher during the pre-task recording phase compared to the task phase (Figure 8A). For ISI standard deviation, neither main effect was significant [drug condition: F(1, 196) = 2.72, p = 0.1004; recording phase: F(2, 379) = 0.1, p = 0.9017], but the interaction between drug condition and recording phase was significant [F(2, 369) = 3.72, p = 0.0252]. Post-hoc assessment of ISI STD revealed only two significant differences: amount of variability is greater during the task compared to either before or after in the vehicle condition (Figure 8B). For mean ISI, neither main effect was significant [drug condition: F(1, 196) = 3.84, p = 0.0514; recording phase: F(2, 379) = 1.82, p = 0.1633], but there was a significant interaction [F(2, 369) = 7.2, p = 0.0009]. Post-hoc assessment for mean ISI revealed three significant differences: neurons fire slower after THC compared to vehicle during post-task recordings, and neurons fire slower during the task compared to pre- and post-task after vehicle administration (Figure 8C). Both the Hurst exponent and multifractal complexity (width) appear to be better indicators of hippocampal microcircuit processing than these standard variability measures since they were able to distinguish task phase and drug condition combinations.
Frequency content varies in a task-specific manner
Independent delta and theta rhythms were found in the human medial temporal lobe and are hypothesized to perform separate cognitive functions (Mormann et al., 2008). Delta power is a classical marker used to assess physiological/sleep state (Harmony et al., 1996; Pereda et al., 1999; Clemens et al., 2013) and is negatively correlated with human working memory performance (Axmacher et al., 2010). Additionally, theta frequency activity is known to coordinate neuronal interactions in working memory networks (Jones and Wilson, 2005), has been correlated with performance during the DNMS task (Hyman et al., 2010, 2011) and is reduced by cannabinoid administration (Ilan et al., 2004; Robbe et al., 2006; Böcker et al., 2010; Kucewicz et al., 2011). Therefore, we assessed the performance of these frequency components in distinguishing between recording phases and drug conditions. Delta and theta power in all hippocampal neurons were assessed by calculating the total power in the delta (0.5–4 Hz) and theta bands (4–8 Hz) using a fast Fourier transform (Figures 4B,D). These two activity bands were chosen because of their prominence in hippocampal spike train recordings. Two repeated measures ANOVAs were performed using delta and theta power as the dependent variable to determine if alterations in frequency content can account for the same differences detected with multifractal analysis.
When using delta power as the dependent variable, a significant main effect of different recording phases was found [F(2, 379) = 5.23, p = 0.0057]. Post-hoc analysis of the recording phase effect showed that delta power is greater during task and post-task compared to pre-task (Figure 8D). There was no significant main effect of drug condition [F(1, 196) = 2.19, p = 0.141] or interaction between drug condition and recording phase [F(2, 369)= 2.52, p = 0.082].
Theta power was assessed in the same manner. The overall effect of drug condition was not significant [F(1, 196) = 3.1, p = 0.0796]. However, there was a significant main effect of recording phase [F(2, 379) = 41.91, p < 0.0001] and a significant drug condition by recording phase interaction [F(2, 369) = 13.28, p < 0.0001]. Post-hoc analysis of the interaction revealed that theta power was significantly higher during both the pre- and post-task recording sessions compared to the task (Figure 8E).
Correlation analysis
The aforementioned results suggest a relationship between spike train temporal coding properties quantified by fractality (Hurst exponent/LRTCs/self-similarity), multifractality (singularity spectrum width), variability measures (mean ISI, ISI STD, CV), and frequency content variables (delta and theta). Correlation analyses were performed in order to improve theoretical understandings of multifractal variables and provide insight for developing computational models of the nervous system that reproduce LRTCs and multifractal complexity of neuronal spike trains. Spearman's rho values were computed for all recording phases and drug condition combinations separately. Correlations between the Hurst exponent and singularity spectrum width yielded small positive rho values ranging from 0.10 to 0.28 (data not shown). This indicates that the Hurst exponent and multifractal width quantify different spike train properties with respect to hippocampal microcircuit processing. Fifty additional correlation analyses were performed to determine relationships of LRTCs (the Hurst exponent) and magnitude of multifractality (width) with CV, ISI STD, mean ISI, delta, and theta power for each recording phase and drug condition. Spearman's rho values greater than |0.5| were considered to indicate “strong” correlations. Pre-task THC values were not reported because they were similar to pre-task vehicle values in all cases.
We tested whether LRTCs were associated with changes in variability or frequency content using correlation. A strong relationship was found between the Hurst exponent and mean ISI during task recordings (Figure 9A), indicating that LRTCs are more prominent with shorter average ISIs (i.e., when neurons fire more frequently). The Hurst exponent was negatively correlated with ISI STD during the task (Figure 9A), indicating that LRTCs and self-similarity occur more frequently when neurons fire with less variable ISIs. Both of these correlations were weaker during the resting phase recordings (pre- and post-task), suggesting a task-specific relationship of fractal vs. firing rate and variability patterns. None of the other three relationships were strong, which indicate that LRTCs and self-similarity quantified by the Hurst exponent are independent from CV, theta and delta power.
To develop a deeper conceptual understanding of multifractal complexity, the singularity spectrum width was correlated with standard variability measures and frequency content. Multifractal complexity (width) was strongly correlated with relative variability, expressed as the coefficient of variation (CV; Figure 9B), during all recording phases and drug conditions. Despite this strong correlation, multifractal complexity was able to distinguish between recording phases and drug conditions better than CV (Figure 7F vs. Figure 8B). This indicates that multifractal complexity is a more sensitive measure for detecting memory-specific alterations in hippocampal microcircuit processing. A strong positive correlation was found between multifractal complexity and delta power during the task that decreased during the resting state recordings (Figure 9B), suggesting that low frequency activity may support multifractal complexity during working memory. Multifractal complexity strongly correlated with “Poison-ness” (CV) and delta power but not with mean ISI, ISI STD, or theta power. Taken together, these results suggest that fractality, multifractality, and theta frequency activity quantify different properties of spike train temporal coding structure.
Discussion
Multifractal analysis revealed something never shown before: the nature of the temporal structure in hippocampal spike trains is altered by performing a working memory task. Monofractality (Hurst exponent) decreased but multifractality (width) increased from pre-task to task and these trends reversed from task to post-task. Prior indications that the neuronal encoding during the task is essential for successful performance (Hampson et al., 2012) support the hypothesis that spike trains represent information with multifractal temporal coding properties (Fetterhoff et al., 2015). Previous results showed that microcircuit activity patterns occurring when encoding information during the sample phase are correlated with correct performance (Berger et al., 2011; Opris et al., 2012) and support the notion that memory functions through repetitive neuronal ensemble codes (Deadwyler and Hampson, 1997; Berger et al., 2012). Therefore, multifractal, self-similar spiking patterns detected during this study might constitute a substrate for memory information transmission.
The main goal of this study was to elucidate the multifractal properties of active hippocampal microcircuit processing by assessing differences between resting state and working memory conditions. The presented results show that multifractal complexity (indicated by singularity spectrum width) permits distinction between all recording phases and support the use of multifractal analysis in extracting variables related to cognitive state better than other commonly applied methods. Our first hypothesis was verified by the finding that long range temporal correlations (LRTCs), indicated by the Hurst exponent, were decreased during task performance compared to the resting state (Figures 6A–C,F). We verified our second hypothesis that multifractal complexity (width) of hippocampal ISI sequences is greater during active memory processing (task recording phase) compared to the resting state (pre- and post-task; Figures 7A–C,F). Unexpectedly, we found that pre-task was more multifractal than post-task (Figure 7C). We believe this could be due to a “priming” effect because the rats are more motivated for water and their hippocampal ensembles are preparing to perform the DNMS task by some form of mental rehearsal. Due to the rats' experience with the habitual daily testing procedure, they may learn to associate pre-task recording conditions with subsequent task recordings when reward is available. The rats are also more satiated during the post-task recording after receiving water during the task and this may reduce the hippocampal inputs responsible for generating multifractality during the task. This finding highlights the advantage of using multifractal analysis for detecting even subtle differences in cognitive states. These results are consistent with other findings showing that task performance elicits multifractal fluctuations restricted to task-related regions while LRTCs decrease in all analyzed brain regions, independent of task involvment (Ciuciu et al., 2012). The structure of variability, detected as multifractal complexity, becomes more intricate when hippocampal microcircuits exhibit a wider range of dynamics required for memory processing (during the task) and this alteration may elicit breakdowns of LRTCs (Figures 6F, 7F). These results suggest that enhancement of the theoretical and conceptual foundations of neuroscience is possible by applying multifractal analysis in order to achieve the larger goal of producing a dynamic portrait of the functioning brain.
THC impairs hippocampal information transmission by disrupting the multifractal spike train patterns that may constitute a basis of memory processing. THC administration reduced multifractal complexity during the DNMS task and decreased LRTCs during post-task recording phase (Figures 6F, 7F). Under control conditions, multifractal complexity (width) and LRTCs became lower and higher, respectively, when rats transitioned from the task to the post-task, but THC administration inhibited these changes. Therefore, it is possible that THC effectively prevented endogenous microcircuit dynamics that facilitate neurophysiological adjustments from memory processing to a resting state. These findings are consistent with others demonstrating that THC promotes default mode network activity at inappropriate times (Bossong et al., 2013), and they support the application of multifractal complexity as a marker of network involvment if the reduction of multifractal complexity is due to decreased memory-related interactions that occur when DNMS performance is impaired by THC.
Computational models attempting to reproduce the temporal coding properties of neuronal spike trains (Goris et al., 2014) can likely be improved by replicating endogenously occuring LRTCs (Fürstenau, 2010) and multifractal complexity. Interpretation of the structure of variability assessed using multifractal analysis was enhanced by comparing performance and relationships between monofractal and multifractal variables (Hurst and width) and standard measures of variability. The correlation analysis showed that the Hurst exponent is negatively correlated with mean ISI and ISI STD during the DNMS task under both vehicle and THC conditions (Figure 9A). This task-specific correlation demonstrates that the increase of LRTCs during resting phase (pre/post) recordings (Figure 6F) is not due to independent changes in mean ISI and ISI STD but may be related to alterations in their ratio when expressed as CV (Figure 9A). The strongest relationship was found between multifractal complexity (width) and the coefficient of variation (Figure 9B), and this finding would suggest that multifractal complexity reflects information directly related to how ISI variability scales with ISI mean. However, the superiority of singularity spectrum width in distinguishing between examined cognitive states is clear in the population ANOVA results: multifractality (width) confirmed four out of six ad hoc distinctions from the interaction between drug condition and recording phase while CV failed to yield a statistically significant interaction and could only distinguish between the pre-task and task recording phases (Figure 7F vs. Figure 8A).
Examination of brain oscillations via frequency content is a commonly used method to study temporal information processing in the nervous system (Lisman, 2005; Axmacher et al., 2006; Dragoi and Buzsáki, 2006). Therefore, it is important to understand how these methods compare to new ways of quantifying and describing temporal coding properties, such as multifractal analysis. It was found that the combination of delta and theta power quantification distinguished between recording phases to a lesser extent than multifractal complexity (width) and LRTCs (Hurst). We found that theta power determination permitted distinction between two out of six ad hoc differences (Figure 8E), and interestingly, both were identical to changes detected by the Hurst exponent (Figure 6F). Further inspection using correlation analyses promoted discovery that the Hurst exponent was not correlated with oscillatory activity in the tested frequency bands and therefore suggests that they are independent dynamical processes. Theta and delta power were negatively and positively, respectively, correlated with multifractal complexity (width). This supports the finding that both delta and theta rhythms exist independently in the human hippocampus (Mormann et al., 2008) and suggests that temporal coding properties detected as multifractality may preferentially arise from the delta frequency activity. It is possible that the precise action potential timing in the delta frequency range, from 1 s to 250 ms, conveys essential information utilized to support cognitive function and reflected as multifractality. It is shown here that multifractal analysis can be used to detect changes in hippocampal processing better than standard Fourier spectrum analyses. However, it will be necessary to compare both methods in order to integrate multifractal concepts with prior findings describing frequency content and to ultimately synthesize new hypotheses about how the brain functions.
Learning and memory require information to be carried from the past into the future, and multifractal complexity of hippocampal neurons may fluctuate depending on how strongly information is received, processed and sent by these neurons. To put our results into a more general perspective, hypothetical singularity spectra were constructed to match our qualitative findings with respect to memory processing (Figure 10). We hypothesize that physiological states are characterized by specific monofractal and multifractal features. All different states must fall within a range of possible dynamics (Figure 10, gray spectrum) and thus a range of possible multifractality. Active memory processing (Figure 10, blue) recruits this system to a stronger degree than resting (Figure 10, orange) and therefore neurons recorded during the task exhibit stronger multifractality. This multifractality may facilitate memory processing by offering a larger range of spike train variability and greater processing capacity. When THC or other memory impairing agents are administered during the task (Figure 10, green), the normal level of multifractal complexity exhibited is reduced and memory performance suffers. Alterations in multifractal complexity may reflect the degree of presently embedded information and therefore would provide information relevant for detection of physiological state.
The application of multifractal analysis revealed that repeating spatiotemporal activity patterns detected in hippocampal spike trains may form a previously undiscovered contribution to the temporal coding of memory information. This demonstrated how analysis of the multifractal structure of temporal dynamics can enable new insights about how the brain functions and can facilitate methodological improvements for detection of alterations in cognitive, physiological, and pathological states (Suckling et al., 2008; Wink et al., 2008). Multifractal analysis is being utilized in pathology diagnosis and has already been successfully applied to Parkinson's Disease (Zheng et al., 2005), seizure detection (Dutta et al., 2014), Alzheimer's Disease (Vysata et al., 2013), and many others (Slezin et al., 2007; Di Ieva et al., 2015). As for physiological state, multifractal analysis was proposed as a method for automatic detection of sleep stages (Weiss et al., 2009; Zorick and Mandelkern, 2013), and an fMRI study showed changes in multifractality selectively occur in task-related brain regions (Ciuciu et al., 2012). The presented assessment of cognitive state can be integrated with these results concerning physiological state: deep sleep contains more LRTCs and less multifractality compared to REM sleep (Weiss et al., 2009), REM sleep contains larger LRTCs than waking (Zorick and Mandelkern, 2013), and as shown here, hippocampal neurons during a task-independent (resting) state exhibit greater LRTCs but less multifractality compared to when recorded during the working memory DNMS task. It is therefore possible that consciousness, described in terms of cognitive/physiological state, occurs on a spectrum that can be quantified and understood via computation of the multifractal singularity spectrum. Multifractal analysis quantifies the scale invariant, self-similar structure that is pervasive throughout biological processes and application of this method to neurophysiological data will improve our understanding of how the nervous system processes cognitive information.
Acknowledgements
The authors thank Dr. Miaohua Jiang for insightful discussion and clear explanations of multifractal analysis. The authors also thank Dr. Andrew Sweatt for his editorial work. This work was supported by contract and grants to SD from the DARPA REMIND Program (N66001-C-09-2080, COL Geoff Ling, Program Manager) and NIDA (DA07625 and DA006634).
Supplementary material
The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fnsys.2015.00130