ORIGINAL RESEARCH
published: 12 April 2019
doi: 10.3389/fnins.2019.00316

Theta Oscillations Alternate
With High Amplitude Neocortical
Population Within
Synchronized States
Erin Munro Krull 1,2* , Shuzo Sakata 3 and Taro Toyoizumi 1
1
RIKEN Center for Brain Science, Tokyo, Japan, 2 Beloit College, Beloit, WI, United States, 3 Strathclyde Instituteof Pharmacy and Biomedical Sciences, University of Strathclyde, Glasgow, United Kingdom

Edited by:
Ritchie Edward Brown,
VA Boston Healthcare System,
United States
Reviewed by:
Kenji Mizuseki,
Osaka City University, Japan
Adrien Peyrache,
McGill University, Canada*Correspondence:
Erin Munro Krull
munroe@beloit.edu
Specialty section:
This article was submitted toSleep and Circadian Rhythms,
a section of the journal
Frontiers in Neuroscience
Received: 23 December 2018
Accepted: 20 March 2019
Published: 12 April 2019
Citation:
Munro Krull E, Sakata S andToyoizumi T (2019) Theta Oscillations
Alternate With High Amplitude
Neocortical Population Within
Synchronized States.
Front. Neurosci. 13:316.
doi: 10.3389/fnins.2019.00316

Synchronized states are marked by large-amplitude low-frequency oscillations in thecortex. These states can be seen during quiet waking or slow-wave sleep. Within
synchronized states, previous studies have noted a plethora of different types ofactivity, including delta oscillations (0.5–4 Hz) and slow oscillations (<1 Hz) in the
neocortex and large- and small- irregular activity in the hippocampus. However, it
is not still fully characterized how neural populations contribute to the synchronizedstate. Here we apply independent component analysis to parse which populations
are involved in different kinds of neocortical activity, and find two populations thatalternate throughout synchronized states. One population broadly affects neocortical
deep layers, and is associated with larger amplitude slower neocortical oscillations. The
other population exhibits theta-frequency oscillations that are not easily observed in
raw field potential recordings. These theta oscillations apparently come from below the
neocortex, suggesting hippocampal origin, and are associated with smaller amplitude
faster neocortical oscillations. Relative involvement of these two alternating populations
may indicate different modes of operation within synchronized states.
Keywords: synchronized state, independent component analysis, theta oscillation, non-REM, slow oscillation,delta oscillation, cortical state

INTRODUCTION
Cortical state can be described along a spectrum from synchronized to desynchronized (Harrisand Thiele, 2011). Desynchronized states are characterized by low amplitude high frequency
activity where local neuronal activity is uncorrelated. Desynchronized activity is seen during
active waking (engaging in a task such as navigating a maze) and rapid-eye movement sleep(REM). In rats, while the neocortex is desynchronized, the hippocampus exhibits theta oscillations,sinusoidal oscillations between 4 and 8 Hz (Diekelmann and Born, 2010). Synchronized states arecharacterized by high-amplitude, low-frequency oscillations where neuronal populations fluctuate
between UP states marked by frequent neuronal firing, and DOWN states marked by neuronalsilence. Synchronized states are seen during quiet waking (resting or engaging in a routine task suchas eating) as well as non-REM sleep. In rats, while the cortex is synchronized, the hippocampus
exhibits irregular activity featuring sharp waves, quickly depolarizing waves generated in CA3
accompanied by ripple oscillations (>100 Hz) (Diekelmann and Born, 2010). Within non-REM

Frontiers in Neuroscience | www.frontiersin.org

1

April 2019 | Volume 13 | Article 316

Munro Krull et al.

Two Modes Alternate Within Synchronized States

Furthermore, we found that during synchronized states when
BL5 is active neocortical UP/DOWN states tend to be slowerand have a larger amplitude, while faster small-amplitude
UP/DOWN states are seen when SUB is active. These results
suggest two modes of operation within the brain during corticalsynchronized states.

sleep, neocortical activity progresses from K-complexes (large
amplitude biphasic waves) to slow-wave activity (SWA) (<1 Hz)
mixed with delta oscillations (1–4 Hz) (Crunelli et al., 2006).
At the same time, the hippocampus can display either lowor high-amplitude irregular activity (Miyawaki et al., 2017).
While activity during synchronized states has been widely
studied, neural ensemble dynamics within synchronized states
are not still fully characterized (Diekelmann and Born, 2010;de Andrés et al., 2011; Crunelli et al., 2015; Neske, 2016;
Mizuseki and Miyawaki, 2017).
Previous studies on cortical dynamics primarily analyzed
activity by looking at either spiking activity or the local fieldpotential (LFP) (Einevoll et al., 2013). Single- or multi-unit
spiking activity gives detailed information about individual
neural interactions. However, because the number of neurons we
can detect is limited, we must infer how cell populations behave asa whole. On the other hand, LFP recordings reflect activity frommany populations simultaneously, including neurons and nonspiking glia (Buzsáki et al., 2012). However, signals from different
populations tend to overlap at each recording site, as illustrated
in Figure 1. The fact that signals overlap makes these recordings
difficult to parse, i.e., it isn’t clear which populations are involved
in different kinds of activity (Kajikawa and Schroeder, 2011;
Buzsáki et al., 2012; Kajikawa and Schroeder, 2015). Current
source density (CSD) analysis helps reduce signal spread byfocusing on current sinks and sources, but it may still bedifficult to parse which cell populations are involved, particularlyin the neocortex where neuronal populations are very dense(Gulyás et al., 2016).
To address the issue of overlapping population signals, weturned to independent component analysis (ICA) (Brown et al.,2001). Based on the assumption that the LFP is made up ofsummed population signals, and that the relative amplitudeof a single population across the electrode remains fixed overtime, ICA decomposes recordings so that signals from different
populations are separated into different components. ICA works
by finding a matrix that decouples the signals so that they are
as statistically independent as possible. Since ICA decouples
signals by matrix multiplication, the number of extracted
components is the same as the number of recording channels.
If there are more populations than recording channels, thenICA groups populations with similar spatial profiles together inthe same component. Even with a limited number of channels,
ICA has proven to be very effective for parsing signals fromdifferent types of neural recordings (Jung et al., 2001; Makarov
et al., 2010; Herreras et al., 2015). Moreover, unlike principal
component analysis which focuses on high amplitude signals,
ICA can easily parse low amplitude signals that may be missedfrom the raw LFP.
In this study, we applied ICA to neocortical laminarrecordings from unanesthetized and urethane anesthetized rats,
and consistently observed two ICA components that alternate
throughout synchronized states. One component has a highamplitude signal that broadly covers much of the lower layers(broad layer 5, BL5); the other component shows intermittent
theta oscillations and has an anatomical profile indicating
it is below the neocortex (THETA) (Herreras et al., 2015).

Frontiers in Neuroscience | www.frontiersin.org

MATERIALS AND METHODS
Animals
We used twenty-three adult Sprague-Dawley rats (male, 270–
413 g, n = 18 for anesthetized experiments, n = 5 forunanesthetized experiments). All procedures for anesthetizedexperiments were performed in accordance with the UK Animals(Scientific Procedures) Act of 1986 Home Office regulations andapproved by the Home Office (PPL60/4217 and 70/8883). Allprocedures for unanesthetized experiments (Sakata and Harris,
2009, 2012) were approved by the Institutional Animal Care and
Use Committee of Rutgers University.

Electrophysiological Experiments
Detailed experimental procedures are described in previous
studies (Sakata and Harris, 2009, 2012; Sakata, 2016). Briefly,
for anesthetized experiments, we anesthetized animals with 1.5–
1.6 g/kg urethane. We also administered lidocaine (2%, 0.1–
0.3 mg) subcutaneously at the site of incision. After attaching ahead-post in the frontal region with bone screws, one of which
was used as an electrode for cortical electroencephalograms(EEGs), we placed the animal in a custom head restraint thatleft the ears free and clear. Two additional bone screws were
implanted in the cerebellum as ground. Body temperature wasmaintained at 37◦ C with a feedback temperature controller (4090-8C, FHC). After reflecting the left temporalis muscle, weremoved the bone over the left auditory cortex (AC) and carefully
performed a small duratomy for each site.
For unanesthetized experiments, in initial surgery, weanesthetized animals with ketamine (100 mg/kg) and xylazine(10 mg/kg), and placed them in a stereotaxic apparatus (David
Kopf Instruments). We attached a head-post (Thorlabs, Inc.) withdental cement (3M ESPE, RelyX Luting Cement), removed theleft temporal muscle, and covered the exposed bone over the left
AC with biocompatible glue and dental cement. After a recovery
period, animals were lightly water-deprived, and handling (5–10 min/day) and head-fixation training began. We trained
the rats for at least five sessions, during which the durationof restraint was gradually extended. Ten percent sucrose was
frequently given during training and water was freely available
for at least 1 h after daily training. On the day of recording, we
carefully performed a craniotomy and duratomy under isoflurane
anesthesia (5% for induction and 0.8% for maintenance). Neither
skin nor muscle was cut during this surgery. After a short
recovery period (>1 h), recording began.
During recording for both experiments, we covered the brain
with 1% agar/0.1 M phosphate buffered saline (PBS) to keep
the cortical surface moist. We inserted a 32 channel siliconprobe (A1x32-10mm-50-177-A32, NeuroNexus Technologies)

2

April 2019 | Volume 13 | Article 316

Munro Krull et al.

Two Modes Alternate Within Synchronized States

FIGURE 1 | Example of ICA applied to neural data. Data generated from a single anesthetized experiment used as an example throughout this study. (A) Different
populations generate field potentials across the electrode, resulting in overlapping signals recorded at each channel. For example, the signal from a single cell (1)may have the highest amplitude in a single channel, while other channels have much lower amplitude. The signal from a larger population (2) may have highamplitude across multiple channels with a single dipole. A population below the electrode (3) may generate the highest amplitude signal in the channel closest to the
population, which decays with channels further from the population and doesn’t display a dipole (Herreras et al., 2015). (B) Example CSD from 32-channel electrodealong with ICA decomposition. The ICA components are listed in order from largest to smallest amplitude over the entire experiment. In the CSD, the auxiliary
channel 0 indicates a tone is played at time 0, and we see an immediate reaction to the stimulus in channels 11 through 17. This immediate reaction to the tone issummarized in ICA component 4. Each ICA component is associated with an anatomical map showing the amplitude of that component across the electrode. The
anatomical map for ICA component 4 is consistent with the large population (2) in (A). Similarly, high frequency activity seen in the CSD in channels 20 and 21 isreflected in ICA component 9, whose anatomical map is consistent with the single cell (1) in (A). Finally, there is a 4 Hz sinusoidal oscillation in ICA component 7 thatis not apparent in the original CSD. The anatomical map gradually increases toward the lowest channel and there is no dipole, indicating that the oscillation was
generated from below the electrode as in population (3) in (A).

where channels are 50 µm apart, either manually or slowly(2 µm/sec or slower) with a motorized manipulator (DMA1511, Narishige) into the AC (1400–1770 µm from the surface).
All electrophysiological experiments were performed in a singlewalled soundproof box (MAC-3, IAC Acoustics) with the interior
covered by 3 inches of acoustic absorption foam. Broadband
signals (0.07–8 kHz) from the silicon probes were amplified (1000times) (Plexon, HST/32V-G20 and PBX3), digitized at 20 kHz
and stored for offline analysis (PXI, National Instruments).
A typical recording schedule was as follows: after insertion of
the probe and an additional waiting period (at least 30 min), we
started recording with a silent period (at least 5 min), followedby sound presentations, and ended with another silent period(at least 5 min).
Acoustic stimuli were generated digitally (sampling rate
97.7 kHz, TDT3, Tucker-Davis Technologies) and delivered infree-field through a calibrated electrostatic loudspeaker (ES1)located approximately 10 cm in front of the animal. We calibrated
the tone presentations using a pressure microphone (PS9200KIT1/4, ACO Pacific, Inc.) close to the animal’s right ear. Acoustic

Frontiers in Neuroscience | www.frontiersin.org

stimuli for this study consisted of short pure tones (50 ms longwith 5 ms cosine ramps, 1/6 or 1/8 octave steps, 3–48 kHz,
10 dB steps, 0–80 dB SPL), long pure tones (300 ms with 10 ms
cosine ramps), and unanesthetized experiments also contained
brief click trains.

Histology
After electrophysiological experiments, we perfused ratstranscardially with physiological saline followed by 4%paraformaldehyde/0.1 M phosphate buffer, pH 7.4. Afterovernight post-fixation in the same fixative, we incubatedbrains in 30% sucrose solution for cryoprotection, cut into
100 µm coronal sections with a sliding microtome (SM2010R,
Leica), and the sections were collected and placed in 0.1 M
PBS. For verification of silicon probe placement, the freefloating sections were counterstained with NeuroTrace (1/500,
N-21480, Life Technologies) in PBS with 0.1% Triton X100 for 20 min at room temperature. The sections were
mounted on gelatin-coated slides and cover-slipped with
antifade solutions.

3

April 2019 | Volume 13 | Article 316

Munro Krull et al.

Two Modes Alternate Within Synchronized States

for all experiments was set to the median LFP-Power overunanesthetized experiments.

Data Analysis Preprocessing
We performed all analysis using Matlab (R2017a, Mathworks,
Waltham, MA, United States). Since the sample rate of theoriginal recordings was 20 kHz, we first filtered the data below
650 Hz using the fast Fourier transform (FFT, fft in Matlab)
with a bump function pass band, where the bump functionsmoothly transitioned from 0 to 1 over 0.25 Hz. We thenremoved 50 Hz line noise by linearly interpolating the Fouriertransform. Finally, we down-sampled the data to 2 kHz, andhigh-pass filtered the data over 0.1 Hz using FFT with a bumpfunction pass band.
We noted three kinds of artifacts to remove from analysis:line drift, electrode pops, and epochs with low firing rates.
To mark epochs with line drift (large amplitude deviationsin the LFP for a single channel), we calculated the standard
deviation of each channel. We then marked sections thatexceeded 5 standard deviations for that channel. For electrodepops (large jumps in the LFP), we calculated histograms ofvoltage slopes over the entire experiment for each channel.
We marked points where the slope was >2,500 V/s for atleast one channel, or >1,000 V/s for at least 2/3 of thechannels within 1 ms. To determine overall firing rate, wesummed the multi-unit activity (MUA) over all channels, andthen computed the average firing rate of the summed MUAover 6 s windows shifted by 600 ms. We marked points thatwere below 40 Hz for anesthetized experiments and 100 Hz
for unanesthetized experiments, where these thresholds wereset by the distribution of firing rates for anesthetized andunanesthetized experiments, respectively.

ICA Application
We first spatially filtered the data using a Gaussian kernel
with width 50 µm. We then calculated the spline inverse
CSD over the LFP (Pettersen et al., 2006). Since the CSD
partially separates sources, previous work suggests that ICA
separates the CSD more cleanly than the LFP (Lȩski et al.,2010). For the ICA training data, we randomly selected
100s worth of data points from synchronized data only. We
used the extended infomax algorithm to minimize effects
of kurtosis in the data (Lee et al., 1999) implemented in
EEGLAB (Schwartz Center for Computational Neuroscience,
La Jolla, CA, United States). We applied ICA 10 times,
each on different training data, to yield 10 different sets of
ICA components. We then chose the ICA component set
that had minimal mutual information between components(Kraskov et al., 2004). To estimate mutual information, we
took samples of 100 to 1,000 s of data points. We thencalculated a mutual information factor (MIF) by estimating
entropy of individual components, using bins with approximately
10 data points per bin. The MIF is the summed entropy of
all individual components, minus ln| W| where W is the ICA
unmixing matrix.
Out of the 32 components generated by ICA, weselected the THETA component first and then the BL5component. ICA lists components according to amplitude.
We chose the first component with a narrow band oscillation
between 2–5 Hz in anesthetized experiments and 5–9 Hz
in unanesthetized experiments, reflecting theta ranges inurethane anesthetized and unanesthetized conditions. The BL5component is the highest amplitude component, not including
the SUB component.
Independent component analysis component amplitude canvary based on experiment as well as populations included in
the component. Therefore, we scaled BL5 and SUB according
to their activity levels. We considered BL5 active when its
signal was correlated with the sMUA. We calculated thecorrelation over 6 s windows shifted by 600 ms to match
spectrogram data. Since the ICA algorithm guesses the signof a component based on the spatial profile, we switched thesign of BL5 if necessary so that the correlation was negative
overall. Time points with correlation ≥0 were considered
inactive. If there were fewer than 30 points that met thiscriterion, then points with correlation ≥−0.1 were considered
inactive. We then centered the signal based on the medianand interquartile range of inactive points. To measure theactivity level of the THETA component, we first determined
the peak frequency between 2–5 Hz for anesthetized data and
5–9 Hz for unanesthetized data. We then smoothed the peakfrequencies using a Gaussian with standard deviation 10 s,
weighting each point according to the amplitude of the peak.
We then found the total amplitude around the peak (width
1 Hz for anesthetized, 2 Hz for unanesthetized), along with theamplitude of the surrounding frequencies which consisted ofa half width window above and below the peak window. The

UP State Detection
We used MUA to determine UP and DOWN states. First
we calculated the smoothed MUA (sMUA) by summing the
MUA over all channels. We then smoothed the summed MUA
using a Gaussian filter with standard deviation of 6.25 ms foranesthetized data and 2.5 ms for unanesthetized data (0.25
divided by minimum firing rate). We then computed UP and
DOWN state transitions based on Sakata and Harris (2009). The
threshold for transitioning from a DOWN to an UP state was
the geometric mean of the sMUA over points where sMUA > 0,
while the threshold for transitioning from an UP to a DOWN
state was 1/5 the UP state threshold. The minimum UP/DOWN
state length was 50 ms for anesthetized data and 25 ms forunanesthetized data.

Spectrogram Calculation
All spectrograms were calculated over 6 s windows shifted by
600 ms using Welch’s power spectral density estimate (pwelch
in Matlab). Spectrograms were used to analyze frequency
content of the LFP as well as ICA components. To estimate
synchronized states, we computed the spectrogram of the
LFP for each channel and then took the median over all
channel spectrograms, keeping time and frequency fixed.
The low-frequency LFP power (LFP-Power) is the summed
power from 1 to 5 Hz (Harris and Thiele, 2011) of themedian spectrogram. The threshold for synchronized states

Frontiers in Neuroscience | www.frontiersin.org

4

April 2019 | Volume 13 | Article 316

Munro Krull et al.

Two Modes Alternate Within Synchronized States

activity level of THETA was defined as the normalized peak
amplitude, where amplitude was centered on the median andinterquartile range of the surrounding frequencies over the entire
experiment for anesthetized data. Because surrounding frequency
amplitude varied a great deal over the course of unanesthetizedexperiments, we decided to normalize the peak amplitude based
and median and interquartile range over 6s instead of thewhole experiment.

RESULTS
We used ICA to investigate which cell populations are active
during synchronized states in 18 urethane anesthetized and 5
unanesthetized head-fixed rats. Anesthetized rats exhibited sleeplike states (Clement et al., 2008; Pagliardini et al., 2013), while
unanesthetized rats could be waking or sleeping. We recorded
from all layers of the AC using 32-channel linear probes (Sakata
and Harris, 2009, 2012). From these recordings, we calculated
the CSD before applying ICA in order to get clearer components(Lȩski et al., 2010). We applied ICA to the CSD 10 times,
each time with different randomly selected sets of artifact-freetraining data, and chose the ICA run with the least estimatedmutual information between components (Kraskov et al., 2004)
to promote as much separation as possible.

Sleep Parameter Threshold Estimation
We estimated sleep stage parameters on the LFP medianspectrogram for unanesthetized experiments. We first calculated
the LFP delta (0.5–4 Hz), theta (6–9 Hz), and sigma (10–14 Hz)
power (Benington et al., 1994; Gross et al., 2009). We thensmoothed the power using a Gaussian kernel with width 2.5 s.
As sigma × theta power indexes waking vs. sleeping, we selected
the sigma × theta threshold first based on largest population
density (Benington et al., 1994). We calculated a histogram with
100 bins over the lowest 90% of data points. Then we smoothed
the histogram using a Gaussian with width over 1 bin. We
calculated the max height of the histogram hc and mean heightof the histogram hm . The threshold is the first point above hc
where the histogram falls below (hc +hm )/2 as sigma × theta
increases. The delta/theta threshold, which indicates the border
between REM and non-REM sleep, was based on the non-REM
to REM indicator value (NIV) (Benington et al., 1994). We
calculated candidate transition to REM (NRT) segments wheredelta power dropped [delta-DP < 1 in Benington et al. (1994)]
and sigma × theta is above threshold. For each NRT segment, we
calculated the NIV and the minimum delta/theta value over theentire segment (DT-min). The threshold was set to the median
DT-min over all segments with NIV values in the top 50%.

Two Populations Consistently Revealed
by ICA
We consistently noticed two ICA components. One component
exhibited sinusoidal oscillations in the theta range (Buzsáki, 2002;
Buzsáki and Draguhn, 2004; Lakatos et al., 2005; Montgomery
et al., 2009) that appeared to originate below the neocortex(THETA). The other was the highest amplitude component
which broadly affected all layers, centered on layer 5 (BL5).
We define the THETA component as the highest-amplitude
component with clock-like sinusoidal oscillations in the theta
range (2–5 Hz anesthetized, 5–9 Hz unanesthetized). We founda THETA component in 16/18 anesthetized experiments and
5/5 unanesthetized experiments. Figure 2A shows that theamplitude of the THETA component did not have a clear
dipole and gradually increased toward the lowest channel
in anesthetized experiments, suggesting that signals volume
conducted from below the neocortex (Sirota et al., 2008; Kajikawa
and Schroeder, 2011; Herreras et al., 2015). In unanesthetizedexperiments, the amplitude also did not have a dipole and waseven over all channels. These anatomical maps were uniquein that they were the only components without a dipole (seeSupplementary Figure S1). At the same time, we noticed thatsignals associated with motion, which affects all channels evenly,
were also separated into the THETA component. (Note that
motion-associated signals may be present outside time intervals
marked as motion artifact, see section “Materials and Methods.”)
Therefore, it is possible that ICA grouped theta oscillations andmotion signals together because of their similar anatomical maps.
We define the BL5 component as the largest amplitude
component, excluding the THETA component. Because BL5
was simply the largest component, each experiment had a BL5component. Anatomical maps of BL5 consistently had large
amplitude across all lower layers, with the largest amplitude in
L5 and a dipole in L4. These anatomical maps are consistent
with a large population generating the signal across multiple
channels as illustrated in Figure 2A. BL5 and THETA wereactive intermittently, as shown in the example anesthetizedexperiment in Figures 2B,C. When BL5 was active, it exhibitedlarge amplitude low frequency oscillations that closely resembled
the CSD. Because of this, we define the BL5 amplitude toLaminar Estimation
To display channel data, we estimated the central L4 channel
based on the maximum CSD amplitude evoked by the preferred
tones. The channel with the highest amplitude was labeled as athalamic recipient layer (L4).

Statistical Analysis for
Neocortical Activity
To compare neocortical activity across states, we used 2-way
ANOVAs based on conglomerated artifact-free data points acrossall anesthetized or unanesthetized experiments with BL5/THETA
index ϕ between −45 and 135 degrees. All 2-way ANOVAs were
performed on the variable of interest as a function of LFP-Powerand ϕ, where LFP-Power was considered continuous and ϕ was
split into BL5 (ϕ ≤ 45) and THETA (ϕ > 45). For UP state
location, we also included the experiment as a random factor.
UP phase amplitude was defined as the largest amplitude
occurrence after Gaussian smoothing over time with a kernelof 2 ms. Peak frequency was defined as the maximum
frequency after Gaussian smoothing with a kernel of
0.25 Hz. Peak frequency width was determined by the Matlab
function findpeaks.

Frontiers in Neuroscience | www.frontiersin.org

5

April 2019 | Volume 13 | Article 316

Munro Krull et al.

Two Modes Alternate Within Synchronized States

FIGURE 2 | Properties of BL5 and THETA components. (A) Spatial maps of both components across 32 channels used to record the LFP, along with diagram of
putative population sources. Results for the example anesthetized experiment shown in (C) and unanesthetized experiment in Supplementary Figure 2 arehighlighted in blue throughout this study. Anatomical maps indicate BL5 is made up of a larger population located within the neocortex, while THETA may begenerated by a population below the neocortex. Amplitude in mA/mm3 , distance in µm. (B) Segments when BL5 is active (strong low-frequency oscillations) andwhen THETA is active (strong sinusoidal oscillation of approximately 4 Hz). Note that the CSD and smoothed multi-unit spiking activity (sMUA) show low-frequency
oscillations in both cases. Scale bar is 85th percentile amplitude of sMUA, CSD in channel 17, and ICA component taken from channel with maximum amplitude inanatomical map. (C) Spectrograms of the LFP, BL5, and THETA components for an example anesthetized experiment. (D) Average frequency spectra when BL5
and THETA are active (scaled amplitude ≥2). Frequencies used to determine activity level are highlighted in pink. (E) Correlation between BL5 normalized amplitudewith LFP-Power and sMUA low frequency power (1–5 Hz), and THETA normalized amplitude with LFP and sMUA relative theta amplitude. Density plot for example
experiments on the left, histograms of correlation from density data on the right. (F) Correlation between BL5/THETA signal and CSD from L5 (channel 350 µm
below L4) and sMUA. Example experiments on left, histograms of median correlations over points where BL5/THETA are active on the right.

be the amplitude from 0.5–4.5 Hz (anesthetized) or 0.5–7 Hz(unanesthetized). When THETA was active, it exhibited a sharp
peak frequency between 2–5 Hz (anesthetized) or 5–9 Hz(unanesthetized). We therefore define the THETA amplitude
based on the peak amplitude within the theta range (seeFrontiers in Neuroscience | www.frontiersin.org

section “Materials and Methods”). Figure 2D shows the averagefrequency spectra of both BL5 and THETA when they are active,
as defined below.
In order to compare activity level across experiments, wenormalized the BL5 and THETA amplitude according to the6

April 2019 | Volume 13 | Article 316

Munro Krull et al.

Two Modes Alternate Within Synchronized States

sMUA (sMUA scaled theta amplitude, T -TEST, ANESTHETIZED:N = 15, P = 0.4623, UNANESTHETIZED : N = 5, P = 0.5536; sMUA,
T - TEST, ANESTHETIZED : N = 15, P = 0.4033, UNANESTHETIZED :N = 5, P = 0.1243) where results for individual experiments can beseen in Figure 2F. Since the THETA component was extracted
from the LFP, we may expect that there is some correlation
between them. However, the fact the THETA is not correlatedwith local spiking activity substantiates that THETA doesn’t haveneocortical origin.

noise level inherent in their signals. We found that BL5 hasa high correlation with the CSD taken from the center of L5
and smoothed Multi-Unit Activity (sMUA) when active, andwas uncorrelated with CSD and sMUA otherwise. Therefore,
we found the median and inter-quartile range of the BL5
amplitude over all time points where BL5 had low correlation
with the sMUA. We then normalized the BL5 component using
this median and inter-quartile range (see section “Materials
and Methods”). Likewise, in the THETA spectrogram inanesthetized experiments, there was a sharp theta peak whichstood out among other frequencies when THETA was active,
and blended in with surrounding frequencies when it was not.
We therefore found the median and interquartile range of thesummed amplitude from frequencies above and below the theta
range over the entire experiment. We then normalized the
THETA amplitude using the surrounding frequency amplitude
median and inter-quartile range (see section “Materials andMethods”). In unanesthetized experiments, while the theta
peak stood out during desynchronized states and intermittently
during synchronized states, the surrounding frequencies changed
according to overall low-frequency LFP power (LFP-Power, 1–
5 Hz power of median spectrogram over all channels) (seeSupplementary Figure 2). In order to not overestimate the
THETA amplitude during synchronized states, we normalized
THETA amplitude according to surrounding frequencies withina 6 s window. Since both BL5 and THETA are normalizedaccording to the median and interquartile range of their inherent
noise level, we consider BL5 and THETA to be active if theiramplitude is ≥2.
Since anatomical maps indicate that BL5 has neocorticalorigin while THETA originates from below the neocortex, weinvestigated the relationship between BL5 and THETA andthe local activity. In order to investigate the relationship, we
looked at both the correlation between the amplitudes of BL5
and THETA with LFP-Power as well as the direct correlation
between the component signals and CSD and sMUA. We did
this because, while the amplitude focuses on a relevant frequency
range, the raw signals may also reveal correlations. We foundthat the BL5 scaled amplitude is correlated with LFP-Power(T -TEST, ANESTHETIZED: N = 17, p < 0.0001, UNANESTHETIZED:N = 5, P < 0.0001) (see Supplementary Table 1 for details
on statistical tests) and sMUA 1–5 Hz amplitude (T -TEST,ANESTHETIZED : N = 17, P < 0.0001, UNANESTHETIZED : N = 5,
P < 0.0062) as illustrated in Figure 2E. When it is active,
BL5 is also directly correlated with the CSD from L5 (mediancorrelation over points with BL5 scaled amplitude > 2, T -TEST,ANESTHETIZED : N = 17, P < 0.0001, UNANESTHETIZED : N = 5,
P < 0.0001) as well as the sMUA (T - TEST, ANESTHETIZED :N = 17, P < 0.0001, UNANESTHETIZED : N = 5, P = 0.0027) asshown in Figure 2F, where the correlation is taken over 6 s
windows in the same manner as the amplitude. Both of thesetests confirm that BL5 is correlated with local activity when BL5
is active.
We found that while THETA was correlated with the LFP(LFP scaled theta amplitude, T -TEST, ANESTHETIZED: N = 15,
P = 0.0260, UNANESTHETIZED : N = 5, P = 0.0532; CSD L5, T - TEST,ANESTHETIZED : N = 17, P < 0.0001, UNANESTHETIZED : N = 5,
P < 0.0001) as seen in Figure 2E, it was not correlated with the
Frontiers in Neuroscience | www.frontiersin.org

THETA Amplitude Alternates With
BL5 Amplitude
Figure 2B appears to show that THETA and BL5 alternatewith each other. To test this possibility, we calculated thecorrelation between the BL5 and THETA scaled amplitudes as
illustrated in Figure 3Ai, and found that they were anticorrelated
over all experiments (Figure 3Aii, T -TEST, ANESTHETIZED:N = 15, P = 0.0002, UNANESTHETIZED : N = 5, P = 0.0054)(see Supplementary Table 2 for details on statistical tests).
Because BL5 and THETA are anti-correlated, we quantified therelative BL5-THETA participation in the CSD by calculating
the angle ϕ and radius r for each BL5-THETA data point asshown in Figure 3Ai. In the conglomerate density plots of
ϕ vs. r over all experiments (Figure 3Aiii), we can see peakdensities where BL5 is dominant (ϕ ≤ 45) and where THETA
is dominant (ϕ > 45). Thus, it makes sense to refer to BL5- andTHETA-dominant states.
We then looked into how long BL5- and THETA-dominantstates last. Figure 3Bi illustrates how ϕ changes over the entire
example anesthetized experiment, where global trends tend
to switch on the order of 100 s of seconds with finer-grainedchanges occurring throughout the experiment. Similar switching
can be seen in unanesthetized experiments, as shown in
Supplementary Figure 3. In Figure 3Bii, we see that althoughBL5-dominant states tend to last longer, both states can laston the order of 10–50 s. The amount of time spent in a BL5dominant vs. THETA-dominant state was similar (Figure 3Biii,
MEAN ± STANDARD DEVIATION (STD) ANESTHETIZED :62.83 ± 16.56%, UNANESTHETIZED: 46.59 ± 8.71%).
Furthermore, tones played during experiments do not appear toaffect the length of states (Figure 3Biv, percent difference during
tone vs. spontaneous activity, MEAN ± STD, ANESTHETIZED:
1.65 ± 18.39%, UNANESTHETIZED: −7.91 ± 14.72%).

BL5 and THETA Amplitudes Alternate
During Synchronized States
In Figure 2B, we see instances where either BL5 or THETA
are active in the presence of noticeable CSD oscillations as seen
during synchronized states. Synchronized states were previously
characterized by high LFP-Power (Harris and Thiele, 2011). Since
we showed in the above section that BL5 and THETA alternate
within experiments, we further investigated whether BL5 andTHETA alternate within synchronized states as well. Figure 4A
shows the density plot of LFP-Power vs. ϕ for the same example
experiments as Figure 2B and Supplementary Figure 2. Within
this plot, ϕ ranges from 0 to 80 degrees for higher values of
LFP-Power in these experiments. Calculating the range (5–95
7April 2019 | Volume 13 | Article 316

Munro Krull et al.

Two Modes Alternate Within Synchronized States

FIGURE 3 | BL5 and THETA alternate within experiments. (A,i) Density plot of BL5 vs. THETA normalized amplitude for an example anesthetized experiment (sameas Figure 1B). Pink and blue marks indicate points that are used for example BL5 and THETA active segments in Figure 1B, respectively. We set ϕ as the angle
from the BL5-axis, while r is the radius. (ii) Correlations of BL5 vs. THETA normalized amplitude for each experiment. (iii) Conglomerate density plots of ϕ vs. radius
for all anesthetized and unanesthetized experiments. (B,i) Time course of ϕ over the same example anesthetized experiment used in Figure 1. (ii) Pie charts
showing the amount of time spent in BL5- vs. THETA-dominant episodes. For example, for anesthetized experiments during THETA-dominant states about 1/3 oftime spent is in an episode of length 10–50 s. (iii) Histogram of the fraction of time spent in BL5- vs. THETA-dominant states. (iv) Histograms of the difference infraction of time dedicated to BL5-dominant states during tone presentations vs. spontaneous activity.

percentile) over the density for all experiments reveals a wide
range of ϕ where both BL5- and THETA-dominant states areincluded in the same echelons of LFP-Power. In contrast, ϕ often
took high values (>80) when the LFP-Power is low.
We then explored how much time is spent in BL5- vs.
THETA-dominant states when LFP-Power is high vs. when
LFP-Power is low. Since many experiments didn’t have a clear
threshold for LFP-Power, we set the threshold as the median
LFP-Power over all unanesthetized experiments (Figure 4A and
Supplementary Figure 4). With this threshold, virtually allanesthetized experiments had high LFP-Power, so we focusedon unanesthetized experiments. The length of low vs. high LFPPower episodes was fairly even, while THETA-dominant states
were longer for low LFP-Power and BL5-dominant states werelonger for high LFP-Power (see Figure 4B). When LFP-Powerwas low most of the time was spent in THETA-dominant states,
while the amount of time spent in BL5- vs. THETA-dominantstates when LFP-Power was high was similar to the ratio foundin anesthetized experiments, with at least 25% of time spent in

Frontiers in Neuroscience | www.frontiersin.org

THETA-dominant states (Figure 4C, MEAN ± STD, PERCENT
HIGH LPF- POWER : 50.33 ± 9.49%, BL5 AND LOW LFP- POWER :19.06 ± 11.58%, BL5 AND HIGH LFP-POWER: 73.73 ± 3.34%).
We further investigated the relationship between THETA
amplitude and LFP-Power directly, since theta oscillations werepreviously noticed predominantly during desynchronized states,
and not synchronized states (Diekelmann and Born, 2010).
Supplementary Figure 5 reveals a similar relationship between
LFP-Power and THETA oscillations as Figure 4A, with a wide
range of amplitudes when LFP-Power is high. Furthermore, thelength of episodes with significant THETA normalized amplitude
during high LFP-Power is approximately the same as episodes
when THETA is inactive (Supplementary Figure 5b). Likewise,
about 50% of the time spent with high LFP-Power had significant
THETA normalized amplitude (Supplementary Figure 5c,
MEAN ± STD, PERCENT LOW SUB THETA , ANESTHETIZED :37.16 ± 20.59%, UNANESTHETIZED: 34.86 ± 9.97%; LOW SUBTHETA AND LOW LFP-P OWER : 23.31 ± 11.66%, LOW SUBTHETA AND HIGH LFP-P OWER : 46.98 ± 6.53%).

8

April 2019 | Volume 13 | Article 316

Munro Krull et al.

Two Modes Alternate Within Synchronized States

FIGURE 4 | BL5 and THETA alternate within synchronized states. (A) Top plots show density of LFP-Power vs. ϕ for single experiments along with 5–95 percentile
ranges. Bottom plots show ranges for all experiments with results for example experiments highlighted in blue. Anesthetized ranges are plotted against centered
LFP-Power (LFP-Power minus LFP-Power with maximum range). Vertical black lines indicate median LFP-Power over unanesthetized experiments, used as thethreshold between low vs. high LFP-Power. Note, by this measure all anesthetized experiments had high LFP-Power. (B) Distribution of episode lengths as apercentage of time spent in episodes of that length. (C) Fraction of time spent in each state. The top panel shows the fraction of time spent with high LFP-Power vs.
low LFP-Power. The lower two panels show the fraction of time spent in BL5- vs. THETA-dominant states. Notably, approximately 30% of time in high LFP-Power is
THETA-dominant. (B,C) Shows results for unanesthetized experiments only.

BL5 and THETA Amplitudes Alternate
Within Sleep Parameters

MEAN ± STD: 5.05 ± 5.35% IN ACTIVE WAKE REGION ,9.76 ± 8.85% IN QUIET WAKE REGION, 30.62 ± 14.47%
IN REM REGION , 68.46 ± 5.40% IN NON -REM REGION ).
Remarkably, the percentage of time spent with THETA-dominantis similar to high LFP-Power states. At the same time,
approximately 50% of time within the estimated non-REM region
has significant THETA amplitude (Supplementary Figure 6,
MEAN ± STD PERCENT LOW SUB THETA : 12.10 ± 10.30%
IN ACTIVE WAKE REGION , 18.32 ± 15.51% IN QUIET WAKEREGION , 26.23 ± 7.34% IN REM REGION , 46.51 ± 6.79% IN
NON -REM REGION ).
One reason for seeing theta during non-REM sleep could
be the transition to REM state, where theta oscillations
become strong just before REM (Benington et al., 1994;
Gottesmann, 1996). To investigate whether THETA oscillations
could represent the transition to REM, we plotted histogramsof BL5- and THETA-dominant points over delta/theta values,
where sigma × theta is above threshold shown in Figure 5C.
We then removed points in the 30 s leading up to thedelta/theta threshold, signifying a transition to REM. The
histograms show that BL5- and THETA-dominated points havea similar distribution above the delta/theta threshold whentransition points are removed. We furthermore checked therelationship of THETA amplitude with these sleep parameters
directly. When we removed transitional points from the nonREM region, we also found similar distributions between high
THETA amplitude vs. low THETA amplitude points. Together,
these figures suggest that THETA-dominant states do not simply
comprise transition-to-REM states.

Since synchronized states may occur during both waking andsleeping, we further characterized the occurrence of BL5- vs.
THETA-dominated states by comparing their amplitudes withtraditional measures of sleep state. In particular, we calculated
the delta (0.5–4 Hz), theta (6–9 Hz), and sigma (10–14 Hz)
power for each unanesthetized experiment. We then compared
ϕ with the measures sigma × theta, used to indicate wake vs.
sleep, and delta/theta, used to indicate non-REM vs. REM sleep(Benington et al., 1994; Gross et al., 2009). Figure 5A showshow ϕ changes for a single unanesthetized experiment withrespect to LFP-Power, delta/theta, and sigma × theta. Although
the first few episodes with high LFP-Power and sigma × theta
appear to be BL5-dominated, the last episode appears to have
many points where THETA is dominant. We then comparedBL5-dominated vs. THETA-dominated points directly on thedelta/theta and sigma × theta plane. In Figure 5B, we can seethat the top right quadrant which indicates non-REM sleep has
clear overlap of both BL5 and THETA points. As expected,
we predominantly see THETA-dominant states when delta/thetaor sigma × theta is low. We divided the delta/theta andsigma × theta plane into regions (see section “Materials andMethods”) in order to estimate the percentage of time dedicatedto BL5 vs. THETA. We found that, while episodes with BL5dominant states tend to be longer within the estimated non-REM
sleep region, a significant amount of time is spent in both BL5and THETA-dominant states, with approximately 70 vs. 30%
of time in each state (Supplementary Figure 6, percent BL5,

Frontiers in Neuroscience | www.frontiersin.org

9April 2019 | Volume 13 | Article 316

Munro Krull et al.

Two Modes Alternate Within Synchronized States

FIGURE 5 | BL5 and THETA dominant states alternate within traditional sleep parameters. (A) LFP-Power, delta/theta power, and sigma × theta power compared toϕ levels for the same unanesthetized experiment as Supplementary Figure 1. (B) THETA-dominant points (ϕ > 45) mapped on top of BL5-dominant points (ϕ≤ 45) on the delta/theta vs. sigma × theta plane used to help distinguish sleep stages. The upper part of the plane indicates sleep, while the lower part indicateswaking. High delta/theta indicates synchronized states such as non-REM sleep, while low delta/theta indicates desynchronized states such as REM sleep. Blackboundaries indicate thresholds used for (C). (C) Histograms of BL5- vs. THETA-dominant points on the left, and high THETA vs. low THETA amplitude on the right.
Delta/theta is scaled by the delta/theta threshold for each experiment. Upper plots include all points above the sigma × theta threshold, set to be just above thepeak sigma × theta density. Middle plots include points at least 30 s away from a transition across the delta/theta threshold. Lower plots include points within 30 s ofa delta/theta threshold crossing. Results for example experiment highlighted in bold.

Neocortical Activity Differs Between
BL5- and THETA-Dominant States

individual UP states. We found that UP state CSD amplitude
significantly decreased with ϕ beyond LFP-Power (TWO-WAY
ANOVA WITH LFP-POWER AND ϕ CONGLOMERATED OVER
EXPERIMENTS , ANESTHETIZED : P < 0.0001, UNANESTHETIZED :P < 0.0001) (see Supplementary Table 3 for details on
statistical tests). Since UP states appeared to be stronger in
L5 during BL5-dominant states, we also tested the UP state
amplitude peak neocortical depth. Similarly, we found that UP
state neocortical depth decreased with respect to ϕ beyond
LFP-Power (TWO-WAY ANOVA BETWEEN LFP-POWER ANDϕWITH EXPERIMENT AS RANDOM VARIATE , ANESTHETIZED :P < 0.0001, UNANESTHETIZED : P = 0.0023).
Finally, we investigated how the frequency content of the
LFP changes with respect to ϕ. Figure 6C shows that the peakfrequency tends to increase for anesthetized experiments, andthe peak in the spectrogram narrows indicating that oscillations
become more regular. We found that this trend is significant over
all experiments independent of LFP-Power (two-way ANOVA
between LFP-Power and ϕ, p = 0 for peak frequency, p = 0.008
for peak frequency width).

Figure 2B shows two versions of slow oscillatory activity inthe L5 CSD and sMUA, one where BL5 is dominant and
one where THETA is dominant. Looking at the CSD andspiking activity across neocortical layers in Figure 6A, wealso notice qualitative differences in the activity between BL5and THETA-dominant states. In particular, we noticed that
BL5 states tended to have larger UP states centered in L5
with a slower frequency. At the same time, THETA states
tended to have smaller UP states focused in L4 with ahigher frequency.
To see if UP state CSD amplitude systematically changed
with ϕ, we compiled the amplitude of all UP states during
synchronized states (LFP-Power greater than median LFPPower over all unanesthetized experiments). We noticed adown-ward trend in amplitude across ϕ in each experiment(see Figure 6B). We then tested whether UP state amplitudedecreased with respect to ϕ independently from LFP-Power,
because LFP-Power may also predict amplitude changes in

Frontiers in Neuroscience | www.frontiersin.org

10

April 2019 | Volume 13 | Article 316

Munro Krull et al.

Two Modes Alternate Within Synchronized States

FIGURE 6 | Neocortical activity differs depending on BL5- vs. THETA-dominant states. (A) CSD (color, mA/mmˆ3) and spiking activity (black asterisks) from example
anesthetized and unanesthetized experiments seen in previous figures. UP state onsets are marked by solid black lines, while DOWN state onsets are marked bydashed lines. (B) Top plot shows density of UP state amplitude and UP state amplitude peak location for individual experiments, along with the mean with respect toϕ. Density is normalized for each ϕ bin, and the mean for a ϕ bin is calculated only if there are ≥ 30 points. We only included synchronized data points. Lower plotsshow the mean for all experiments, with highlighted results for the example experiment above. UP peak location is centered along the mean location for eachexperiment (UP location minus UP location mean). (C) Top left plot shows average spectrogram values over ϕ for the same example experiments as above, alongwith the mean peak frequency over ϕ. Other plots show means for all experiments, with highlighted results for the example experiment.

Frontiers in Neuroscience | www.frontiersin.org

11

April 2019 | Volume 13 | Article 316

Munro Krull et al.

Two Modes Alternate Within Synchronized States

during non-REM sleep can be as strong as theta power
during REM sleep (Gaztelu et al., 1994; Pedemonte et al.,2005), but these oscillations may be more difficult to observe
because of ongoing irregular activity. Finally, theta oscillations
can be seen with increased tiredness during waking states(Vyazovskiy and Tobler, 2005).
Independent component analysis is a form of analysis that
separates overlapping signals in a set of recordings (Bell andSejnowski, 1995; Amari et al., 1996; Hyvärinen and Oja, 1997;
Hyvärinen and Oja, 2000; Brown et al., 2001; Herreras et al.,2015). Previously, ICA has been used a great deal in artifact
removal as well as identifying regions of interest in fMRI andEEG (McKeown et al., 2003; Onton et al., 2006; Eichele et al.,2009; Esposito and Goebel, 2011; Huster et al., 2015; Calhoun andde Lacy, 2017). More recently, ICA was used in LFP recordings
to parse population involvement in the hippocampus (Makarov
et al., 2010; Agarwal et al., 2014; Herreras et al., 2015; Benito et al.,2016). To our knowledge, this is the first study to use ICA to study
synchronized states.

The frequency content of unanesthetized experiments proved
to be richer than anesthetized experiments, with strong low
frequency oscillations throughout synchronized states and anincreasing envelope of delta frequencies. We tested whether theenvelop of delta frequencies increased for all experiments bymeasuring the peak frequency above 2.5 Hz. We saw that thepeak frequency increased significantly with ϕ beyond LFP-Power(TWO-WAY ANOVA WITH LFP-POWER AND ϕ, p < 0.0001).
The peak frequency below 2.5 Hz did not change significantlywith ϕ. However, the peak frequency width appeared to narrow
with ϕ, indicating that oscillations became more regular. We
saw that the peak frequency width decreased significantly over
all unanesthetized experiments both below 2.5 Hz (TWO-WAY
ANOVA WITH LFP-POWER AND ϕ, P < 0.0001) and above 2.5 Hz(TWO-WAY ANOVA WITH LFP-POWER AND ϕ, P < 0.0001),
indicating that both facets of the frequency spectra became
more regular.

DISCUSSIONInterpretation of Results With
ICA and Limitations

We applied ICA to recordings from rat AC and consistently
found two components that alternate within synchronized states.
The BL5 component has a high amplitude and follows the overall
activity. The THETA component generates a theta oscillation
which appears to orginate below the neocortex. Neocortical UP
state activity varies depending on which component is active.
These results indicate that there are two modes of operation
within synchronized states.

One reason why ICA has not been used more on LFP recordings
to date may be because LFP recordings have access to relatively
detailed information, so there are more options available foranalysis than fMRI and EEG. At the same time, ICA is limitedin the fact that the number of components it can return islimited by the number of recordings. Because there are potentiallymore populations detected by neural recordings than can possibly
be parsed, it has been unclear how ICA would behave when
faced with real neural data. Recent simulations in Makarov
et al. (2010) and Gła̧bska et al. (2014) show that ICA tends
to group similar populations together, which is consistent with
our results. With more recordings, ICA would be able to pickout more sources at finer resolution (Agarwal et al., 2014;
Jun et al., 2017).
As ICA groups populations with similar anatomical profiles
together, BL5 and THETA may not stand for single neural
populations. In fact, ICA will parse any signal that affects
recordings, including glia and movement artifacts. ICA
components also may be composed of sets of synapses onto
cells acting synchronously (Herreras et al., 2015). Furthermore,
it is possible that two populations that act independently
at some time points, and synchronously at others, may beseparated into three components: one component for each of
the populations individually and a third for the combination(Herreras et al., 2015).
Although we consistently found a THETA component
which likely exhibited sub-neocortical theta oscillations, thetheta oscillations were more difficult to see in unanesthetizedexperiments compared to anesthetized experiments during
synchronized states. In particular, the theta oscillation wasoften mixed with other low-frequency activity in unanesthetizedexperiments, while the theta peak was isolated in all anesthetizedexperiments. The theta-frequency peak may have been clearer
for anesthetized THETA components because: (1) activity
under anesthesia is simpler than activity without anesthesia

Related Studies on Differences
Within Synchronized States andApplication of ICA
There are several other studies that note different modes
within synchronized states. Recently, Miyawaki et al. (2017)
describe low-amplitude synchronized states (LOW) seen in the
hippocampus, which builds on previously noted small-irregular
activity. LOW states may be related to our THETA-dominantstates since UP state CSD amplitude is lower during THETAdominant states. Notably, LOW states occur simultaneously inthe hippocampus, entorhinal cortex, neocortex, and thalamus.
Therefore, if there is a correlation between LOW states and
THETA-dominant states, then low-amplitude irregular activity
should be seen in the hippocampus during synchronized
THETA-dominant states. Thus, our study hints that there areunderlying theta oscillations during low-amplitude states in
the hippocampus.
Although rat hippocampal theta oscillations have only
been noted during desynchronized states in most studies,
there are a few studies that look at theta oscillations during
synchronized states. Hippocampal theta is seen in non-REM
sleep during transition-to-REM states shortly before REM
starts (Benington et al., 1994; Robert et al., 1999). However,transition-to-REM lasts up to 10 s (Gottesmann, 1996), whereas
THETA-dominant states can last significantly longer and do
not necessarily lead to REM sleep. Several studies havenoted that theta power in the hippocampus and thalamus

Frontiers in Neuroscience | www.frontiersin.org

12

April 2019 | Volume 13 | Article 316

Munro Krull et al.

Two Modes Alternate Within Synchronized States

involvement during SWA (Szabó et al., 2017). Therefore, BL5
may include glia.
THETA represents a population exhibiting sub-neocorticaltheta oscillations. Theta oscillations are clock-like sinusoidal
oscillations previously reported in the hippocampus along withother limbic structures (Buzsáki, 2002). Since the hippocampus
lies directly below the rat auditory neocortex, then THETA
may originate from the hippocampus or possibly deeper
limbic structures (Sirota et al., 2008; Kajikawa and Schroeder,
2011). At the same time, the theta oscillations we observed
had a lower frequency, matching type 2 theta oscillations
previously associated with sensory processing and emotionallysalient stimuli during waking states (Bland and Oddie, 2001;
Tendler and Wagner, 2015). If type 2 theta oscillations areassociated with emotionally relevant memory processing,then they could also be linked to hippocampal-amygdala coactivation (Pelletier and Paré, 2004; Hutchison and Rathore,2015). It has also been hypothesized that theta oscillations
correspond to information coming in to the hippocampus,
while sharp waves indicate information transferred out of
the hippocampus (Hasselmo, 1999). THETA oscillations
during synchronized states may mean that there are lowamplitude theta oscillations within the hippocampus at thesame time as sharp waves, where the theta oscillations standfor information coming in from another structure such as
the amygdala, while information is begin transferred out
via sharp waves.
Further research would involve verification of the rolesof various populations in BL5- and THETA-dominant states.
Being able to focus on one state within synchronized states
may also help to shed light on not only mechanisms foroscillatory activity, but functional roles for activity. For instance,
memory processing has been shown to be hampered fromsleep deprivation for the latter half of the night (Rasch
and Born, 2013). Previously, this was tied to a lack of
REM sleep which plays a larger role later in the night.
However, stage 2 non-REM sleep frequently alternates with
REM in the latter half of the night and further study
may elucidate a unique role for this stage in memory
processing (Diekelmann and Born, 2010; Rasch and Born, 2013;
Hutchison and Rathore, 2015).
We may also use this analysis to find different modes ofoperation, and specifically theta oscillations, in other species. For
example, while rat hippocampal theta is easy to pick out, theta
oscillations tend to be more difficult to identify and study in otherspecies, such as humans and bats (Jacobs, 2013).

and so ICA could better isolate the theta oscillations inanesthetized experiments, or (2) the THETA ICA component inunanesthetized experiments included movement artifact whichallowed neocortical signal to leak into the component. In the firstcase, the THETA component may even isolate theta oscillations
from multiple sources that originate below the cortex. For
instance, the hippocampus, which is directly below the AC andis known to show mixed frequencies during synchronized states(Buzsáki, 2002).
Another limitation of our study is that instead of sleep
staging, we compared ICA components to sleep stage parameters.
Thresholds for these parameters are based on the data fromrat recordings lasting 1–1.5 h, which may not include all
sleep stages. This means that sleep regions may not necessarily
line up with sleep stages seen in other experiments. While
our analysis does not provide a direct comparison with sleep
stages, a comparison to the parameters allows us to see how
BL5 and THETA signals relate to each other regardless of
how sleep stages may be decided. Likewise, while there is no
exact threshold for LFP-Power to determine synchronized vs.
desynchronized states, we are still able to compare LFP-Power
with BL5 and THETA signals.

Possible Implications for SynchronizedStates and Further Research
Our results show that BL5 UP states have relatively highamplitude, while THETA UP states have lower amplitude.
Amplitude during sleep has previously been characterized by
SWA, which is high at sleep onset and slowly subsides throughout
the night (Dijk et al., 1990; Aeschbach and Borbély, 1993;
Vassalli and Dijk, 2009; Rodriguez et al., 2016). Because ofthis change, we hypothesize that BL5-dominant states may beemphasized during early sleep, while THETA-dominant states
may be emphasized later, consistent with LOW states (Miyawaki
et al., 2017). According to the two-process model of sleep, early
sleep is associated with homeostatic processes along with SWA(Borbély, 1982; Achermann et al., 1993; Borbély et al., 2016).
At the same time, sleep is hypothesized to serve two functions:synaptic homeostasis and memory consolidation (Vassalli andDijk, 2009; Diekelmann and Born, 2010; Rasch and Born, 2013).
In particular, late sleep and lower amplitude stage 2 sleep has
been shown to benefit certain types of memory consolidation(Diekelmann and Born, 2010; Rasch and Born, 2013; Hutchison
and Rathore, 2015).
If BL5 and THETA reflect different modes within
synchronized states, then what populations could they represent?
BL5 encompasses much of the cortex, is relatively highamplitude, and mimics the CSD. This means that BL5 could
comprise a large population within the cortex, a high amplitude
population, or both. Early work shows that glial recordings
have relatively high amplitude and follow the LFP (Amzica
and Steriade, 1998). Glia are also known to be involved inhomeostatic processes during sleep (Clinton et al., 2011; deAndrés et al., 2011; Poskanzer and Yuste, 2011; Fellin et al., 2012;
Bjorness et al., 2016). Furthermore, recent work with calcium
imaging shows heavy glial involvement prior to neuronal

Frontiers in Neuroscience | www.frontiersin.org

CONCLUSION
Independent component analysis reveals that theta oscillations
exist within synchronized states, and that they alternate with ahigh-amplitude population in the neocortex. These populations
are associated with different neocortical activity. Therefore,
these data indicate that synchronized states are not uniform,
but may be separated according to at least two distinct
functional roles.

13

April 2019 | Volume 13 | Article 316

Munro Krull et al.

Two Modes Alternate Within Synchronized States(EMK), JSPS KAKENHI Grant No. JP18H05432 (TT), AMED
under Grant No. JP19dm020700 (TT), Medical Research Council(MR/J004448/1 to SS), and Biotechnology and Biological Sciences
Research Council (BB/K016830/1 and BB/M00905X/1 to SS).

DATA AVAILABILITY
The datasets generated for this study can be access at https://doi.org/10.15129/1225d3ad-9728-43a9-8091-5259f091f7a7 or areavailable on request to the corresponding author.

ACKNOWLEDGMENTS
AUTHOR CONTRIBUTIONS
We would like to thank Tansi Khodai for performing theanesthetized experiments and Austin Brockmeier for helpful
discussions with ICA. Unanesthetized data was originallypublished in Sakata and Harris (2009, 2012).

EMK: conceptualization, formal analysis, and writing – originaldraft. EMK and TT: methodology. SS: in vivo experiments. TT
and SS: resources. EMK, TT, and SS: writing – review and editing.
TT: supervision.

SUPPLEMENTARY MATERIALFUNDING
The Supplementary Material for this article can be foundonline at: https://www.frontiersin.org/articles/10.3389/fnins.
2019.00316/full#supplementary-material

This work was supported by the RIKEN Brain Science Institute(EMK and TT), Beloit College Keefer Junior Faculty Funds

REFERENCES

Brown, G. D., Yamada, S., and Sejnowski, T. J. (2001). Independent component
analysis at the neural cocktail party. Trends Neurosci. 24, 54–63. doi: 10.1016/
S0166-2236(00)01683-0
Buzsáki, G. (2002). Theta oscillations in the hippocampus. Neuron 33, 325–340.
doi: 10.1016/S0896-6273(02)00586-X
Buzsáki, G., Anastassiou, C. A., and Koch, C. (2012). The origin of extracellular
fields and currents - EEG, ECoG, LFP and spikes. Nat. Rev. Neurosci. 13,
407–420. doi: 10.1038/nrn3241
Buzsáki, G., and Draguhn, A. (2004). Neuronal oscillations in cortical networks.
Science 304, 1926–1929. doi: 10.1126/science.1099745
Calhoun, V. D., and de Lacy, N. (2017). Ten key observations on the analysisof resting-state functional MR imaging data using independent component
analysis. Neuroimag. Clin. N. Am. 27, 561–579. doi: 10.1016/j.nic.2017.
06.012Clement, E. A., Richard, A., Thwaites, M., Ailon, J., Peters, S., and Dickson,
C. T. (2008). Cyclic and sleep-like spontaneous alternations of brain state
under urethane anaesthesia. PLoS One 3:e2004. doi: 10.1371/journal.pone.00
02004
Clinton, J. M., Davis, C. J., Zielinski, M. R., Jewett, K. A., and Krueger, J. M. (2011).
Biochemical regulation of sleep and sleep biomarkers. J. Clin. Sleep Med. 7(5Suppl.), S38–S42. doi: 10.5664/JCSM.1360
Crunelli, V., Cope, D., and Hughes, S. (2006). Thalamic T-type Ca2+ channels
and NREM sleep. Cell Calcium 40, 175–190. doi: 10.1016/j.ceca.2006.
04.022
Crunelli, V., David, F., Lorincz, M. L., and Hughes, S. W. (2015).
The thalamocortical network as a single slow wave-generating
unit. Curr. Opin. Neurobiol. 31, 72–80. doi: 10.1016/j.conb.2014.
09.001de Andrés, I., Garzón, M., and Reinoso-Suárez, F. (2011). Functional anatomy ofnon-REM sleep. Front. Neurol. 2:70. doi: 10.3389/fneur.2011.00070
Diekelmann, S., and Born, J. (2010). The memory function of sleep. Nat. Rev.
Neurosci. 11, 114–126. doi: 10.1038/nrn2762
Dijk, D.-J., Brunner, D. P., and Borbély, A. A. (1990). Time course of EEG power
density during long sleep in humans. Am. J. Physiol. 258(3 Pt 2), 650–661. doi:
10.1097/WNP.0b013e3181879d56
Eichele, T., Calhoun, V. D., and Debener, S. (2009). Mining EEG-fMRI using
independent component analysis. Int. J. Psychophysiol. 73, 53–61. doi: 10.1016/
j.ijpsycho.2008.12.018
Einevoll, G. T., Kayser, C., Logothetis, N. K., and Panzeri, S. (2013). Modelling andanalysis of local field potentials for studying the function of cortical circuits.
Nat. Rev. Neurosci. 14, 770–785. doi: 10.1038/nrn3599
Esposito, F., and Goebel, R. (2011). Extracting functional networks with spatial
independent component analysis: the role of dimensionality, reliability and

Achermann, P., Dijk, D.-J., Brunner, D. P., and Borbély, A. A. (1993). A modelof human sleep homeostasis based on EEG slow-wave activity: quantitative
comparison of data and simulations. Brain Res. Bull. 31, 97–113. doi: 10.1016/
0361-9230(93)90016-5
Aeschbach, D., and Borbély, A. A. (1993). All-night dynamics of thehuman sleep EEG. J. Sleep Res. 2, 70–81. doi: 10.1111/j.1365-2869.1993.
tb00065.x
Agarwal, G., Stevenson, I. H., Berényi, A., Mizuseki, K., Buzsáki, G.,
and Sommer, F. T. (2014). Spatially distributed local fields in the
hippocampus encode rat position. Science 344, 626–630. doi: 10.1126/science.
1250444
Amari, S., Cichocki, A., and Yang, H. (1996). A new learning algorithm
for blind signal separation. Adv. Neural Inform. Process. Syst. 8,
757–763.
Amzica, F., and Steriade, M. (1998). Electrophysiological correlates of sleep delta
waves. Electroencephalogr. Clin. Neurophysiol. 107, 69–83. doi: 10.1016/S00134694(98)00051-0
Bell, A. J. A. J., and Sejnowski, T. J. T. J. (1995). An information-maximization
approach to blind separation and blind deconvolution. Neural Comput. 7,
1129–1159. doi: 10.1162/neco.1995.7.6.1129
Benington, J. H., Kodali, S. K., and Heller, H. C. (1994). Scoring transitions
to REM sleep in rats based on the EEG phenomena of Pre-REM sleep: an
improved analysis of sleep structure. Sleep 17, 28–36. doi: 10.1093/sleep/
17.1.28
Benito, N., Martin-Vazquez, G., Makarova, J., Makarov, V. A., and Herreras, O.
(2016). The right hippocampus leads the bilateral integration of
gamma-parsed lateralized information. eLife 5, 1–26. doi: 10.7554/eLife.
16658
Bjorness, T. E., Dale, N., Mettlach, G., Sonneborn, A., Sahin, B., Fienberg,
A. A., et al. (2016). An adenosine-mediated glial-neuronal circuit for
homeostatic sleep. J. Neurosci. 36, 3709–3721. doi: 10.1523/JNEUROSCI.390
6-15.2016
Bland, B. H., and Oddie, S. D. (2001). Theta band oscillation and synchrony inthe hippocampal formation and associated structures: the case for its role insensorimotor integration. Behav Brain Res. 127, 119–136. doi: 10.1016/S01664328(01)00358-8
Borbély, A. A. (1982). A two process model of sleep regulation. Hum. Neurobiol. 1,
195–204. doi: 10.1111/jsr.12371
Borbély, A. A., Daan, S., Wirz-Justice, A., and Deboer, T. (2016). The two-processmodel of sleep regulation: a reappraisal. J. Sleep Res. 25, 131–143. doi: 10.1111/
jsr.12371

Frontiers in Neuroscience | www.frontiersin.org

14

April 2019 | Volume 13 | Article 316

Munro Krull et al.

Two Modes Alternate Within Synchronized Statesaggregation scheme. Curr. Opin. Neurol. 24, 378–385. doi: 10.1097/WCO.
0b013e32834897a5
Fellin, T., Ellenbogen, J. M., De Pittà, M., Ben-Jacob, E., and Halassa,
M. M. (2012). Astrocyte regulation of sleep circuits: experimental andmodeling perspectives. Front. Comput. Neurosci. 6:65. doi: 10.3389/fncom.2012.
00065
Gaztelu, J. M., Romero-Vives, M., Abraira, V., and Garcia-Austt, E. (1994).
Hippocampal EEG theta power density is similar during slow-wave sleep
and paradoxical sleep. A long-term study in rats. Neurosci. Lett. 172, 31–34.
doi: 10.1016/0304-3940(94)90655-6
Gła̧bska, H., Potworowski, J., Łȩski, S., and Wójcik, D. K. (2014). Independent
components of neural activity carry information on individual populations.
PLoS One 9:e105071. doi: 10.1371/journal.pone.0105071
Gottesmann, C. (1996). The transition from slow-wave sleep to paradoxical sleep:
evolving facts and concepts of the neurophysiological processes underlying theintermediate stage of sleep. Neurosci. Biobehav. Rev. 20, 367–387. doi: 10.1016/
0149-7634(95)00055-0
Gross, B. A., Walsh, C. M., Turakhia, A. A., Booth, V., Mashour, G. A., andPoe, G. R. (2009). Open-source logic-based automated sleep scoring software
using electrophysiological recordings in rats. J. Neurosci. Methods 184, 10–18.
doi: 10.1016/j.jneumeth.2009.07.009
Gulyás, A. I., Freund, T. F., and Káli, S. (2016). The effects of realistic synaptic
distribution and 3D geometry on signal integration and extracellular field
generation of hippocampal pyramidal cells and inhibitory neurons. Front.
Neural Circ. 10:88. doi: 10.3389/fncir.2016.00088
Harris, K. D., and Thiele, A. (2011). Cortical state and attention. Nat. Rev. Neurosci.
12, 509–523. doi: 10.1038/nrn3084
Hasselmo, M. E. (1999). Neuromodulation: acetylcholine and memory
consolidation. Trends Cogn. Sci. 3, 351–359. doi: 10.1016/S1364-6613(99)01365-0
Herreras, O., Makarova, J., and Makarov, V. A. (2015). New uses of LFPs: Pathwayspecific threads obtained through spatial discrimination. Neuroscience 310,
486–503. doi: 10.1016/j.neuroscience.2015.09.054
Huster, R. J., Plis, S. M., and Calhoun, V. D. (2015). Group-level component
analyses of EEG: validation and evaluation. Front. Neurosci. 9:254. doi: 10.3389/
fnins.2015.00254
Hutchison, I. C., and Rathore, S. (2015). The role of REM sleep theta activity inemotional memory. Front. Psychol. 6:1439. doi: 10.3389/fpsyg.2015.01439
Hyvärinen, A., and Oja, E. (1997). A fast fixed point algorithm for independent
component analysis. Neural Comput. 9, 1483–1492. doi: 10.1162/neco.1997.9.7.
1483
Hyvärinen, A. A., and Oja, E. E. (2000). Independent component analysis:
algorithms and applications. Neural Netw. 13, 411–430. doi: 10.1016/S08936080(00)00026-5
Jacobs, J. (2013). Hippocampal theta oscillations are slower in humans thanin rodents: implications for models of spatial navigation and memory.
Philosoph. Trans. R. Soc. B Biol. Sci. 369:20130304. doi: 10.1098/rstb.2013.
0304
Jun, J. J., Steinmetz, N. A., Siegle, J. H., Denman, D. J., Bauza, M., Barbarits, B.,
et al. (2017). Fully integrated silicon probes for high-density recording of neuralactivity. Nature 551, 232–236. doi: 10.1038/nature24636
Jung, T. P., Makeig, S., Westerfield, M., Townsend, J., Courchesne, E., andSejnowski, T. J. (2001). Analysis and visualization of single-trial eventrelated potentials. Hum. Brain Mapp. 14, 166–185. doi: 10.1002/hbm.
1050
Kajikawa, Y., and Schroeder, C. E. (2011). How local is the local field potential?
Neuron 72, 847–858. doi: 10.1016/j.neuron.2011.09.029
Kajikawa, Y., and Schroeder, C. E. (2015). Generation of field potentials
and modulation of their dynamics through volume integration ofcortical activity. J. Neurophysiol. 113, 339–351. doi: 10.1152/jn.00914.
2013
Kraskov, A., Stögbauer, H., and Grassberger, P. (2004). Estimating mutual
information. Phys. Rev. 69(6 Pt 2):66138. doi: 10.1103/PhysRevE.69.066138
Lakatos, P., Shah, A. S., Knuth, K. H., Ulbert, I., Karmos, G., and Schroeder,
C. E. (2005). An oscillatory hierarchy controlling neuronal excitability andstimulus processing in the auditory cortex. J. Neurophysiol. 94, 1904–1911.
doi: 10.1152/jn.00263.2005

Frontiers in Neuroscience | www.frontiersin.org

Lee, T. W., Girolami, M., and Sejnowski, T. J. (1999). Independent component
analysis using an extended infomax algorithm for mixed subgaussianand supergaussian sources. Neural Comput. 11, 417–441. doi: 10.1162/
089976699300016719
Lȩski, S., Kublik, E., Swiejkowski, D. A., Wróbel, A., and Wójcik, D. K. (2010).
Extracting functional components of neural dynamics with independent
component analysis and inverse current source density. J. Comput. Neurosci.
29, 459–473. doi: 10.1007/s10827-009-0203-1
Makarov, V. A., Makarova, J., and Herreras, O. (2010). Disentanglement of local
field potential sources by independent component analysis. J. Comput. Neurosci.
29, 445–457. doi: 10.1007/s10827-009-0206-y
McKeown, M. J., Hansen, L. K., and Sejnowsk, T. J. (2003). Independent
component analysis of functional MRI: What is signal and what is
noise? Curr. Opin. Neurobiol. 13, 620–629. doi: 10.1016/j.conb.2003.
09.012
Miyawaki, H., Billeh, Y. N., and Diba, K. (2017). Low activity microstates during
sleep. Sleep 40:zsx066. doi: 10.1093/sleep/zsx066
Mizuseki, K., and Miyawaki, H. (2017). Hippocampal information processing
across sleep/wake cycles. Neurosci. Res. 118, 30–47. doi: 10.1016/j.neures.2017.
04.018Montgomery, S. M., Betancur, M. I., and Buzsaki, G. (2009). Behavior-dependent
coordination of multiple theta dipoles in the hippocampus. J. Neurosci. 29,
1381–1394. doi: 10.1523/JNEUROSCI.4339-08.2009
Neske, G. T. (2016). The slow oscillation in cortical and thalamic networks:
mechanisms and functions. Front. Neural Circ. 9:88. doi: 10.3389/fncir.2015.
00088
Onton, J., Westerfield, M., Townsend, J., and Makeig, S. (2006). Imaging human
EEG dynamics using independent component analysis. Neurosci. Biobehav. Rev.
30, 808–822. doi: 10.1016/j.neubiorev.2006.06.007
Pagliardini, S., Gosgnach, S., and Dickson, C. T. (2013). Spontaneous sleeplike brain state alternations and breathing characteristics in urethane
anesthetized mice. PLoS One 8:e70411. doi: 10.1371/journal.pone.007
0411
Pedemonte, M., Gambini, J. P., and Velluti, R. A. (2005). Novelty-induced
correlation between visual neurons and the hippocampal theta rhythm insleep and wakefulness. Brain Res. 1062, 9–15. doi: 10.1016/j.brainres.2005.
07.069Pelletier, J. G., and Paré, D. (2004). Role of amygdala oscillations in the
consolidation of emotional memories. Biol. Psychiatry 55, 559–562. doi: 10.
1016/j.biopsych.2003.08.019
Pettersen, K. H., Devor, A., Ulbert, I., Dale, A. M., and Einevoll, G. T.
(2006). Current-source density estimation based on inversion of electrostatic
forward solution: effects of finite extent of neuronal activity and conductivity
discontinuities. J. Neurosci. Methods 154, 116–133. doi: 10.1016/j.jneumeth.
2005.12.005
Poskanzer, K. E., and Yuste, R. (2011). Astrocytic regulation of cortical UP
states. Proc. Natl. Acad. Sci. U.S.A. 108, 18453–18458. doi: 10.1073/pnas.111
2378108
Rasch, B., and Born, J. (2013). About sleep’s role in memeory. Physiol. Rev. 93,
681–766. doi: 10.1152/Physrev.00032.2012
Robert, C., Guilpin, C., and Limoge, A. (1999). Automated sleep staging
systems in rats. J. Neurosci. Methods 88, 111–122. doi: 10.1016/S0165-0270(99)00027-8
Rodriguez, A. V., Funk, C. M., Vyazovskiy, V. V., Nir, Y., Tononi, G.,
and Cirelli, C. (2016). Why does sleep slow-wave activity increase after
extended wake? assessing the effects of increased cortical firing during
wake and sleep. J. Neurosci. 36, 12436–12447. doi: 10.1523/JNEUROSCI.161
4-16.2016
Sakata, S. (2016). State-dependent and cell type-specific temporal processing
in auditory thalamocortical circuit. Sci. Rep. 6:18873. doi: 10.1038/srep
18873
Sakata, S., and Harris, K. D. (2009). Laminar structure of spontaneous and sensoryevoked population activity in auditory cortex. Neuron 64, 404–418. doi: 10.
1016/j.neuron.2009.09.020
Sakata, S., and Harris, K. D. (2012). Laminar-dependent effects of cortical state on
auditory cortical spontaneous activity. Front. Neural Circ. 6:109. doi: 10.3389/
fncir.2012.00109

15

April 2019 | Volume 13 | Article 316

Munro Krull et al.

Two Modes Alternate Within Synchronized States

Sirota, A., Montgomery, S., Fujisawa, S., Isomura, Y., Zugaro, M., and Buzsáki, G.
(2008). Entrainment of neocortical neurons and gamma oscillations by the
hippocampal theta rhythm. Neuron 60, 683–697. doi: 10.1016/j.neuron.2008.
09.014
Szabó, Z., Héja, L., Szalay, G., Kékesi, O., Füredi, A., Szebényi, K., et al.
(2017). Extensive astrocyte synchronization advances neuronal coupling inslow wave activity in vivo. Sci. Rep. 7, 1–18. doi: 10.1038/s41598-01706073-7
Tendler, A., and Wagner, S. (2015). Different types of theta rhythmicity are induced
by social and fearful stimuli in a network associated with social memory. eLife
2015, 1–22. doi: 10.7554/eLife.03614
Vassalli, A., and Dijk, D.-J. (2009). Sleep function: current questions and new
approaches. Eur. J. Neurosci. 29, 1830–1841. doi: 10.1111/j.1460-9568.2009.
06767.x

Frontiers in Neuroscience | www.frontiersin.org

Vyazovskiy, V. V., and Tobler, I. (2005). Theta activity in the waking EEG is amarker of sleep propensity in the rat. Brain Res. 1050, 64–71. doi: 10.1016/j.
brainres.2005.05.022
Conflict of Interest Statement: The authors declare that the research was
conducted in the absence of any commercial or financial relationships that could
be construed as a potential conflict of interest.
Copyright © 2019 Munro Krull, Sakata and Toyoizumi. This is an open-access article
distributed under the terms of the Creative Commons Attribution License (CC BY).
The use, distribution or reproduction in other forums is permitted, provided theoriginal author(s) and the copyright owner(s) are credited and that the original
publication in this journal is cited, in accordance with accepted academic practice. No
use, distribution or reproduction is permitted which does not comply with these terms.

16

April 2019 | Volume 13 | Article 316