Theta/delta coupling across cortical laminae contributes to semantic cognition

Theta/ delta coupling across cortical laminae contributes to semantic cognition. J Neurophysiol 121: 1150–1161, 2019. activity in populations of neurons is associated with cognitive and motor function. Our understanding of the neuronal mechanisms underlying these core brain functions has beneﬁtted from demonstrations of cellular, synaptic, and network phenomena, leading to the generation of discrete rhythms at the local network level. However, discrete frequencies of rhythmic activity rarely occur alone. Despite this, little is known about why multiple rhythms are generated together or what mechanisms underlie their interaction to promote brain function. One overarching theory is that different temporal scales of rhythmic activity correspond to communication between brain regions separated by different spatial scales. To test this, we quantiﬁed the cross-frequency interactions between two dominant rhythms—theta and delta activity—manifested during magnetoencephalography recordings of subjects performing a word-pair semantic decision task. Semantic processing has been suggested to involve the formation of functional links between anatomically disparate neuronal populations over a range of spatial scales, and a distributed network was manifest in the proﬁle of theta-delta coupling seen. Furthermore, differences in the pattern of theta-delta coupling signiﬁcantly correlated with semantic outcome. Using an established experimental model of concurrent delta and theta rhythms in neocortex, we show that these outcome-dependent dynamics could be reproduced in a manner determined by the strength of cholinergic neuromodulation. Theta-delta coupling correlated with discrete neuronal activity motifs segregated by the cortical layer, neuronal intrinsic properties, and long-range axonal targets. Thus, the model suggested that local, interlaminar neocortical theta-delta cou- pling may serve to coordinate both cortico-cortical and cortico-subcortical computations during distributed network activity. NEW & NOTEWORTHY Here, we show, for the ﬁrst time, that a network of spatially distributed brain regions can be revealed by cross-frequency coupling between delta and theta frequencies in subjects using magnetoencephalography recording during a semantic decision task. A biological model of this cross-frequency coupling suggested an interlaminar, cell-speciﬁc division of labor within the neocortex may serve to route the ﬂow of cortico-cortical and cortico-subcortical information to promote such spatially distributed, func- tional networks. (2–5 M (cid:6) and were bandpass ﬁltered at 0.1 Hz to 300 to provide a local reference for the neocortically generated delta rhythm. Intracellular recordings were performed with 2 M potassium acetate-ﬁlled glass micropipettes (50–150 M (cid:6) ) and were low-pass ﬁltered to 2.5 kHz. Neuron subtypes were quantiﬁed by response to a depolarizing current step (0.2 nA, 200 ms) from resting membrane potential. Statistical analysis of excitatory postsynaptic potentials (EPSPs) and action potential (AP) outputs was performed with a repeated-measures ANOVA to examine the effects of experimental condition (different concentrations of carbachol) over delta phase. Signiﬁcant [using the false discovery rate (FDR) test] coupling between the ﬁeld delta rhythm and EPSP/AP proﬁles is displayed in polar plots (Fig.


INTRODUCTION
Rhythmic electrical activity in discrete frequency bands accompanies a broad range of motor, affective, and cognitive processes in the brain. Evidence for a direct, causal role in cognitive processing has been postulated for some time (e.g., Başar et al. 2001), with different frequencies involved to different extents in different tasks. For example, gamma rhythms organize primary sensory information to facilitate higher-order processing (Fries 2015); beta rhythm generation correlates with task performance, requiring short-term memory and prediction (Arnal and Giraud 2012); alpha rhythms control access to stored memories (Klimesch 2012); theta rhythms appear to be required for sequential processing of sensory information (Remondes and Wilson 2013); delta rhythms appear vital for semantic processing (Harmony 2013). Where the information is available, these rhythms appear to have an origin in subsets of neurons in local cortical and thalamic circuits (Roopun et al. 2008), but the temporal organization they impart onto local circuit outputs is vital for control of information flow within the wider cortical mantle (e.g., Akam and Kullmann 2010;Li et al. 2017).
Having a large library of rhythms available to local cortical circuits makes for a complex temporal landscape. However, the situation is further complicated by observations showing that discrete rhythms are rarely generated alone. Coexistence of multiple frequencies of activity, each temporally interacting with one another, is a common feature of brain dynamics (Lakatos et al. 2005), but we understand little about the mechanisms that facilitate these interactions, nor the computational advantages they may impart: Do they just represent a simple additive process, with one brain region involved simultaneously in multiple cortical processes, or is there synergy at work, with the presence of multiple, interacting rhythms superadditive for cortical function? There is increasing evidence for the latter, leading to the current working hypothesis that different frequencies of rhythm chaperone cortical communication on different spatial scales (Canolty and Knight 2010; Kopell et al. 2000).
To address this issue, we use a semantic cognition task to quantify changes in the outcome-dependent pattern of interaction between the two cardinal rhythms involved (theta and delta frequency activity). Semantic cognition assigns meaning to the storm of sensory inputs that we experience during wakefulness (Corbett et al. 2009) and has been shown to involve interaction between many brain regions on multiple spatial scales (Binder et al. 2009). Functional neuroimaging studies suggest a "hub and spoke"-like structure to semantic representation networks (Patterson et al. 2007). Activity locally seen within hub regions (e.g., Mollo et al. 2017) provides amodal rendering of semantic information (e.g., anterior temporal lobe, ATL) and comparators to prior schema leading to appropriate behavioral responses [e.g., posterior middle temporal gyrus (pMTG), intraparietal sulcus (IPS)]. This local activity needs to be distributed over broader spatial scales for appropriate semantic cognition to take place. Interaction between the above hub regions and language centers [left inferior frontal gyrus (LIFG); Sharot et al. 2012], executive control regions (Duncan 2010), and parietal areas that are involved in polymodal sensory integration, attention, and spatial cognition is also required [angular gyrus (AG); Binder et al. 2009].
Flow of information between different regions involved in the semantic network most likely depends upon theta and delta frequency's temporal patterns. Evidence for a role for theta-alpha rhythms (5-12 Hz) in semantic processing has been seen, in general, in both the cortex (Klimesch et al. 1994) and, specifically, in the ATL (van Ackeren and Rueschemeyer 2014). Indirect evidence also suggests an involvement of activity at the delta frequency. The N400 component of event-related potentials-thought to represent a phase-resetting of ongoing delta rhythms (Van Petten and Luka 2006)-is a marker for semantic cognition (Koelsch et al. 2004) and communication between frontal and parietal regions of semantic relevance and is mediated by delta rhythms (1-4 Hz) during decision tasks (Nácher et al. 2013). The decision-making component of semantic processing also suggests importance for theta and delta rhythms. Activity across the theta-delta bands (2-9 Hz) relates to models of activity in which neurons coding for a particular outcome receive increasing levels of excitation over time (i.e., a "ramp" of synaptic input) (Hunt et al. 2012;Purcell et al. 2010;Wang 2002). In addition, iterative switching between discrete neuronal activity states at approximately theta frequency correlates with decision making in frontal regions (Rich and Wallis 2016).
Given these proposed roles for theta and delta frequency activity in semantic cognition, we then use a local circuit, biological model of these rhythms (Carracedo et al. 2013) to identify possible "cellular" origins for their pattern of interaction. The profile of theta/delta interaction changes during the semantic task could be accurately reproduced in this local circuit model. The known connectivity properties of the cellular substrates identified by the model support the "different rhythms-different spatial scales" hypothesis (Canolty and Knight 2010) and suggest further refinement: Different frequencies of rhythm interact to coordinate cortico-cortical and cortico-subcortical functional connectivity.

Participants
Seventeen healthy, female, right-handed native English speakers (mean age: 23.3 yr, age range: 20 -35 yr) with normal or correctedto-normal vision participated in this study under written informed consent and following ethical approval (York Neuroimaging Centre, University of York).

Task
The task is shown schematically in Fig. 2A. Stimuli were adapted from Badre et al. (2005) and involved sequential visual presentation of word pairs with varying levels of association. Subjects were requested to consider whether the word pairs were semantically related. Subjects were not aware of this aspect of the task structure. Related word combinations were selected using the Edinburgh Associative Thesaurus to identify words that were produced frequently together by 22% of participants. Unrelated word combinations were created by shuffling these data sets and removing coincidentally related pairs. The words were all nouns with a concreteness rating Ͼ500 as specified in the Medical Research Council psycholinguistic database (Wilson 1988). Each word appeared once in both the related and unrelated conditions. Red fixation lines were always present on the screen, with words shown in light gray on a dark gray background (as in Fig. 2A). The screen was positioned 75 cm from the viewer, and the image was projected in a way that words subtended 1°vertically and 5°horizontally at the retina. First word presentation occurred 800 ms into the task for a duration of 200 ms. The second word presentation occurred at 1,150 ms, again for 200 ms. The trial length was 3,550 ms with a jitter of between 0 and 1,000 ms between adjacent trials. Ten percent of trials were catch trials, in which participants were cued via a question mark on the screen to press a button to indicate whether the words in the pair were related. These trials were used to monitor reaction time (RT) performance and the behavioral response of interest: semantic decision ("related" versus "unrelated") for each subject/word pair. The catch trial and following trial were excluded from further analyses to avoid contamination from motor artifacts.

Data Acquisition
High-resolution T1-weighted MR structural scans (method adapted from Kozinska et al. 2001) were acquired on a GE 3.0T Signa Excite HDX with an image resolution of 1.13 mm and used for magnetoencephalography (MEG) coregistration. MEG data were acquired using a 4D Neuroimaging Magnes 3600 Whole Head 248 Channel SQUID magnetometer-based system. Coregistration was achieved with a Polhemus Fastrak (Colchester, VT), and five fiducial coils allowed for movement tracking. Data were acquired in a dark, magnetically isolated room and at a sample rate of 678.17 Hz (low-pass filtered to 200 Hz). Raw MEG data are available on request from E. Jefferies (e-mail: beth.jefferies@york. ac.uk).

Virtual Electrodes
"Virtual electrodes" were constructed using a linearly constrained minimum variance beamformer using the technique described by Van Veen (1997). The spatial filters were constructed using a covariance matrix calculated from all of the data and a multiple-spheres forward model (Huang et al. 1999). Following coregistration with each subject's structural MRI scan, coordinates were specified in standard Montreal Neurological Institute (or MNI) space, with individual coordinate spaces transformed to standard space. The first principal component analysis component of the three orientations was used to create time series data. The extraction of beamformed data, and the analyses below, was performed using the York Neuroimaging Analysis Framework (publicly available at https://vcs.ynic.york.ac.uk/ docs/naf/index.html) and standard MATLAB signal processing routines.

Selection of Regions of Interest
A 3D lattice of points was constructed across the whole brain with 5-mm spacing, and beamformers were used to compute the total power at each point using the Neural Activity Index (NAI) (Van Veen 1997). This approach was used to assess any changes in magnitude of the two frequencies of interest (theta and delta rhythms; Fig. 3B). However, owing to combinatorial issues and computational demand, more detailed spatiotemporal analyses were performed on a much smaller number of sites chosen for their precedented involvement-in terms of task-and state-dependent significant functional MRI (fMRI) activations-in visual, primary sensory (Bankó et al. 2011), and semantic processing (e.g., Visser et al. 2010). Regions involved in control, default, memory, and network switching tasks (e.g., Christoff et al. 2009) were also included for comparison.
Chosen coordinates for left (L) and right (R) hemispheres were as follows (numbered in bold, according to the matrices shown in

Analysis
Preprocessing. For each trial/subject, beamformed channels containing movement artifacts, eye blinks, and external noise were excluded, leaving 70 -130 trials for each condition (related or unrelated word pairs) for each subject. Data epochs for each trial (length 3,350 ms, see "task" above) were then filtered from 0.5 to 45 Hz using a zero-phase, finite impulse response filter (MATLAB) to preserve both phase information and the original sample rate (678.17 Hz) for analysis of temporal structure. For analyses focusing on theta and delta rhythms, we first identified the subject mean modal peak frequencies in the time series data epochs from 1st stimulus presentation to RT. These were 5.2 Hz (used as the theta frequency reference) and 2.5 Hz (delta frequency) (see Fig. 3A). For theta/delta phase-amplitude coupling analyses, band-pass-filtered time series were generated from 4.5 to 7.5 Hz and from 1.0 to 4.0 Hz.
Region activation. Changes in event-related potentials (ERPs), particularly, the N400 component, have been linked to semantic processing previously (Koelsch et al. 2004). Therefore, we quantified stimulus-induced deviations from baseline (the first 800 ms of each trial before presentation of the first stimulus), as magnitude changes above threshold in the virtual electrode data. Threshold was defined by analyzing the distribution for all virtual electrode activity for each subject and set at two standard deviations above the baseline mean (i.e., above the 5% significance level, Fig. 1A). This method captured both the initial ERP and any subsequent, significant deviations leading up to the decision time (i.e., 800 ms; RT for each trial/subject. Regions shown to be active during the task, in general (at any time from presentation of 1st stimulus to RT), were then subdivided according to each subject's semantic decision ("related" or "unrelated") on a trial-by-trial basis.
Temporal cross-covariances. Two measures of pairwise temporal relationships between virtual electrode activity were used on the data epoch between first stimulus presentation to RT. First, synchrony is recognized as a key mechanism used by the cortex to code for properties of sensory stimuli (Gray et al. 1989). We calculated pairwise synchrony as the mean magnitude of the cross-covariogram around 0 ms (Ϫ25 ms to ϩ25 ms) lag for each pair of electrodes ( Fig.  1B) using the signal processing toolbox in MATLAB. Data were Hamming windowed with window length of 300 ms and overlap of 285 ms. Second, covariance between electrode pairs at theta frequency-particularly, with half a theta cycle temporal separation-has been demonstrated as a signature of functional connectivity between brain regions (Mizuseki et al. 2009). Therefore, we took the magnitudes of the pairwise cross-covariogram around ϩ96 ms and Ϫ96 ms (192 ms ϭ 5.2 Hz, the mean theta frequency for this subject cohort). These were meaned between 25 ms on either side of these phase values. Pairwise synchrony and theta covariance thresholds were again set at two SDs over the mean cross-covariance data distribution.
Phase-amplitude coupling. More complex temporal signatures, particularly cross-frequency coupling, are known to play an important role in long-range network function during cognition (Canolty and Knight 2010). Therefore, we focused on interactions between the two dominant frequencies expressed between stimulus presentation and RT. Delta-theta coupling was quantified using the method described by Kramer and Eden (2013), using meaned traces from each region node for each subject. Electrode data were band-pass filtered to extract activity in the delta band and theta band based on the cohort mean delta and theta frequencies (e.g., see above and Fig. 3A). The absolute component of the Hilbert transform of the theta channel was plotted relative to the phase of each delta period between stimulus presentation and response time (2-5 delta periods per subject/electrode). This method exposed any changes in theta rhythm magnitude relative to the concurrent delta rhythm phase (Fig. 1C). Analysis data are presented on a linear plot to demonstrate any changes in the phase-amplitude coupling (PAC) profile for each task outcome. A false-detection rate analysis (Storey and Tibshirani 2003) was then performed to extract any statistical differences from shuffled theta channel data (P Ͼ 0.05), and the significant PAC scores were displayed on polar plots (e.g., Fig. 4 and Fig. 7). Raw data for this analysis were selected to meet the criteria suggested by Aru et al. (2015) for cross-frequency analysis in all aspects except "input-related nonstationarities"-phase reset could not be discounted as it is a fundamental property of the delta and theta rhythms and vital for their role in cognition (Calderone et al. 2014, see discussion; Cobb et al. 1995).
In vitro model. Cognitively relevant phenomena in human EEG, MEG, and fMRI data sets are often modeled to suggest underlying network mechanisms. Neural, computational, and statistical models have been applied to semantic cognition (Ralph et al. 2017) but rarely inform with sufficient objectivity and biological detail to suggest celland local network-level mechanisms. This is particularly pertinent when considering the origin of brain dynamic signatures of cognitive relevance for two reasons: First, the minimum network substrates underlying known mechanisms of brain rhythms are all, to date, contained in local circuits anatomically smaller than the maximum spatial and/or temporal resolution of noninvasive techniques. Second, the different neuronal subtypes involved in generation of different rhythms have different anatomical projection profiles, which limitand thus point toward-their possible involvement in communication across multiple, distributed regions (see DISCUSSION). Coupled delta and theta rhythms have been modeled at the local circuit level previously in rodent association cortical regions anatomically corresponding to areas known to be involved in human semantic cognition (Carracedo et al. 2013). In this model, the two coexistent rhythms were generated locally by subnetworks involving neurons with very different long-range projection profiles: regular spiking neurons (RS, theta) projecting purely cortico-cortically and intrinsic bursting neurons (IB, delta) projecting subcortically (Groh et al. 2010;Kim et al. 2015). Therefore, we used this model to identify possible neuronal substrates for semantic decision-related changes in theta/delta coupling in local circuits. We then related model data to the "different frequencies-different spatial scales" hypothesis.
Briefly, horizontal and thalamocortically oriented slices (450 m thick) containing parietal cortex were taken from male Wistar rats (150 -300 g), according to the UK Home Office Animals (Scientific Procedures) Act 1986. Slices were transferred to an interface chamber and perfused with oxygenated artificial cerebrospinal fluid (aCSF) (in mM: 126 NaCl, 3 KCl, 1.25 NaH 2 PO 4 , 0.6 mM MgSO 4 , 1.2 CaCl 2 ,24 NaHCO 3 , and 10 glucose) at 32°C. Delta rhythms occurred spontaneously and persistently in the presence of 2-6 M carbachol-an analog of ACh. Extracellular field recordings were performed with aCSF-filled glass micropipettes (2-5 M⍀) and were bandpass filtered at 0.1 Hz to 300 Hz to provide a local reference for the neocortically generated delta rhythm. Intracellular recordings were performed with 2 M potassium acetate-filled glass micropipettes (50 -150 M⍀) and were low-pass filtered to 2.5 kHz. Neuron subtypes were quantified by response to a depolarizing current step (0.2 nA, 200 ms) from resting membrane potential. Statistical analysis of excitatory postsynaptic potentials (EPSPs) and action potential (AP) outputs was performed with a repeated-measures ANOVA to examine the effects of experimental condition (different concentrations of carbachol) over delta phase. Significant [using the false discovery rate (FDR) test] coupling between the field delta rhythm and EPSP/AP profiles is displayed in polar plots (Fig. 7).

Basic ERP Electrophysiology Revealed Sensory Task-Dependent But not Semantic Interpretation-Dependent Brain Regions
No significant differences in behavioral performance were seen between each semantic interpretation (P Ͻ 0.05, Fig. 2,  A and B). Similarly, the magnitude of the suprathreshold ERP changes also showed no correlation with behavioral performance. A spatiotemporal sequence of regional activations was seen for each subject distributed around the time of the second stimulus presentation (Fig. 2C, i). However, although no set of regions active at any time between stimulus presentation and RT was common to all subjects when analyzed by semantic outcome (word-pairs seen as related or unrelated), a set of regions showed activation (Ͼ2 SDs above baseline) in response to task, in general, in all subjects (Fig. 2C, ii). These regions were all posterior and consisted of bilateral primary visual, posterior cingulate, and parietal areas, along with both hippocampi, although measurement from this deep-lying structure with noninvasive methods is notoriously unreliable (Fig. 2C, iii). No involvement of the frontal or temporal semantic regions used for analysis was exposed using this basic measure of local activity. Thus, region activation alone failed to capture the semantic decision-making process.

Theta/Delta Spectral Content Dominated the Epoch between Stimulus Presentation and Response but Did not Associate with Semantic Interpretation
The above data demonstrated that different semantic decision outcomes did not associate with magnitude changes in ERP in the individual regions analyzed here. However, data also showed dominant, spectrally discrete delta and theta frequency peaks during the time from stimulus presentation to behavioral response (Fig. 3A, i, ii). Band-pass filtering of these delta (1-4 Hz) and theta (4.5-7.5 Hz) components revealed no significant global differences in power in the time period from presentation of the stimuli to a "related" and "unrelated" behavioral response when averaging across all subjects (Fig.  3B, P Ͼ 0.05). The presence of rhythmic activity is a known substrate for communication between brain regions vital for cognition (Fries 2015). While rhythm power may be relevant in some cases, the key property of such rhythms is to provide a temporal framework with which to coordinate activity over cortical distance. Therefore, we next examined basic measures of this temporal framework associated with delta and theta rhythms.

Pairwise Regional Covariance Metrics also Revealed Taskbut not Semantic Interpretation-Dependent Networks
We measured the degree of synchrony from broadband signals as a measure of pairwise regional communication between regions (Gray et al. 1989). Significant (P Ͻ 0.05) pairwise regional synchrony values for each correlation mean showed sequences of coactivation from 1st word presentation to behavioral response for each subject (Fig.  3C, i). A small number of local region pairs in frontal and parietal areas correlated with task for Ͼ15/17 subjects (94%, Fig. 3C, i and iii). However, no long-range correlations were apparent (Fig. 3C, iii) nor, again, did synchrony between any set of region pairs correlate with either "related" or "unrelated" semantic outcome for all subjects (Fig.  3C, i).
The above observations demonstrated that both the strongest ERP responses and the most broadband-synchronous region pairs did not correspond with either the core semantic network structure revealed by fMRI studies (see INTRODUC-TION), or the nature of the semantic decision made on presentation of the word pairs. Therefore, we next examined the precedented pattern of theta-mediated functional connectivity-the half theta period phase separation of activity between regions (Mizuseki et al. 2009). Again, no common set of region pairs correlated with semantic interpretation (Fig. 3D, i). However, compared with the 0-ms synchrony metric (Fig. 3C), a greater number of region pairs were shared by Ͼ15/17 subjects (94%) when considering semantic task alone. Connected regions formed a distributed network involving medial prefrontal cortex, medial frontal gyrus, hippocampus, visual cortex, and temporal areas. Unlike the region pairs revealed by absolute synchrony, these regions have been shown previously to be vital for semantic processing: ATL-aSTG, pMTG, IFG, dmPFC, IPS, and AG (Whitney et al. 2011;Jefferies 2013), Fig. 3D, iii); for executive control: vmPFC and MFG; memory: both left and right hippocampi; and task-associated primary sensory processing: left and right V1 (Fig. 3D, iii). The resulting network was far more global than that seen for zero lag/lead synchrony alone (Fig. 3C, iii versus Fig. 3D, iii) and consisted of region pairs with highly significant, greater anatomical separation (P Ͻ 0.01, n ϭ 6, 15).

PAC Between Theta and Delta Rhythms Predicted Semantic Interpretation
The above data established that temporal interactions associated with theta rhythms exposed a network of regions previously linked to semantic cognition. To investigate this further, we next examined any relationship between the theta rhythmic activity and the other dominant component of the spectra revealed in Fig. 3A-the delta rhythm. Delta and theta rhythms show cross-frequency coupling during cognition (Lakatos et al. 2005), and cross-frequency coupling, in general, is thought to provide a temporal framework for coordinating local and distal cortical dynamics (Canolty and Knight 2010). Therefore, we analyzed the profile of theta coupling to the delta rhythm for each region. We used the mean, by region for each subject, delta-rhythm phase from the left hemisphere only as a reference, as we wanted to prevent possible interference from interhemispheric phase Fig. 2. Raw response magnitude demonstrates common regions related to subjective task performance but does not predict outcome. A: pictorial examples from the task. B: percentage errors and reaction times after second stimulus for all subjects for the two conditions-related (black) and unrelated (red). C, i: an example of region activity changes with time. Data from one subject, averaged for "related" trials is displayed as colormap where regions form each row (arranged by time to maximal response). Trial time is along the x-axis and magnetoencephalography (MEG) region activations (Ͼ2 SD) are on the color axis. Times for stimulus onset are shown as dashed lines. Note region numbers (y-axis) were reorganized by time of initial activation and do not represent the numbers given in methods. C, ii: plot shows the number of active regions (Ͼ2 SDs from the mean signal from 0.8 to 2.0 s) across subject. Note no regions were common to more than 7/17 subjects when considering outcome ("related"black line, "unrelated"-red line, "related" and "unrelated"-blue line). C, iii: map shows the location of each of the 10 regions found common to Ͼ15 subjects (94%) overlaid on a horizontal brain representation viewed from below.
lags. The left hemisphere was chosen because previous studies have suggested a dominance of the left hemisphere in semantic processing (Vigneau et al. 2006;Whitney et al. 2011). Plotting the absolute (magnitude) or the theta rhythm against this mean delta phase metric showed a switch from single to bimodal peaks per delta period before a "related" and "unrelated" interpretation respectively. These data suggested a difference in local delta-theta phase amplitude coupling for different semantic interpretations (P Ͻ 0.01, Kologmorov-Smirnov, n ϭ 17; Fig. 4). Fig. 3. Synchrony and theta covariance differentially demonstrate common region pairs related to subjective task performance but do not predict outcome. A: average waveform (left) and spectrograms (right) for four example regions [left posterior cingulate cortex (LPCC), left angular gyrus (LAG), right anterior insular cortex (RAIC), and right posterior middle temporal gyrus (RpMTG)] averaged over all trials for one example in one subject for the related (black) and unrelated (red) conditions. Scale bar: 200 ms, 2 nA·m. B: mean (by subject) difference between "related" and "unrelated" delta and theta power shown as colormap of percent difference calculated for each point of grid of beamformed nodes with 5-mm spacing throughout the brain. Note that none of these regional differences was significant when corrected for multiple comparisons [family-wise error (FWE): P Ͻ 0.05 in specialized proresolving mediators (SPM)]. C, i: pairwise, normalized cross-covariance between regions demonstrated trajectories of synchrony above threshold (Ͼ2 SDs) from stimulus presentation (inset). Note region numbers (y-axis) were reorganized by time of initial activation and do not represent the numbers given in METHODS. Main graph shows number of region pairs with correlations Ͼ2 SDs above mean for each behavioral outcome. No synchronous region pairs were common to Ͼ8 subjects when separated into outcome ("related"-black line, "unrelated"-red line). 0-ms synchrony showed only five region pairs interacted in Ͼ15 subjects when considering task (blue line). C, ii: cross-covariance matrix of all region pairs, color-mapped onto commonality across the subject pool. C, iii: location of each of the five region pairs common to Ͼ15 (94%) of subjects. D, i: theta covariance between regions demonstrated trajectories from stimulus presentation (inset). Main graph shows number of region pairs with theta-lagged correlations Ͼ2 SDs above mean for each behavioral outcome. Theta covariance showed 14 region pairs interacted in Ͼ15 subjects when considering task (blue line). D, ii: cross-covariance matrix of all region pairs, color-mapped onto commonality across subjects. Dark colors compared with light colors indicate commonality in greater numbers of subjects. D, iii: location of each of the 14 region pairs common to Ͼ15 (94%) or subjects. Inset shows the distribution of interregional distances for interacting regions in each case. Note the significant (P Ͻ 0.05) shift from short-to long-distance functional connections when considering theta covariance (orange) versus synchrony (purple).

A Biological Model Predicted the Changes in PAC Profile Correlated with Altered Neocortical Neuronal Subtype Outputs
It is notoriously difficult to assign a genuine biological process to cross-frequency coupling phenomena using analytical models (Aru et al. 2015). To establish a possible neuronal network mechanism to the above modification of delta-theta PAC (Fig. 4), we used an in vitro biological model shown to capture delta-theta PAC previously. The model chosen generates local, concurrent neocortical theta and delta rhythms in isolated parietal cortex (Carracedo et al. 2013). The delta rhythm responses seen in the MEG data sets used here appeared as a consequence of stimulus-locked averaging of ongoing delta rhythms in individual trials ("phase resetting," Fig. 5, A and B), as seen previously for visual stimuli (Klein et al. 2016). It should perhaps be noted that delta rhythms are a prominent feature of both slow-wave sleep and waking (Hall et al. 2014), with amplitude differences relative to higher frequencies responsible for the subtle changes in the spectral power law behavior (Wen and Liu 2016). Therefore, we first established that the model could also demonstrate this phase resetting. Using thalamocortical slices, time-discrete (50 s) thalamic stimulation (bipolar electrode, 10 -100 A) was extremely effective at phase resetting the persistent neocortical delta rhythm (Fig. 5, C and D).
However, while phase reset is a fundamental response of ongoing cortical rhythms to input, it is also a potential source of spurious cross-frequency coupling measures (Aru et al. 2015). Therefore, we designed a set of model conditions that attempted to capture potential differences in cortical dynamics associated with semantic interpretation in the absence of discrete inputs. The in vitro model is dependent on low, but non-zero, cholinergic neuromodulation. ACh levels in the cortex vary on a subsecond timescale and carry information about the subjective nature of sensory stimuli. This form of neuromodulation is positively correlated with stimulus presentation properties, the degree of novelty, saliency, and the subject's degree of attention to the stimulus (Acquas et al. 1996). Thus, we used different levels of cholinergic neuromodulation to explore the effects on the profile of concurrent theta-delta rhythm interactions in the three main principal neuron subtypes in deep and superficial neocortex (Fig. 6A).
Lower and higher levels of pharmacologically induced cholinergic excitation (2 M and 6 M carbachol, respectively) generated significantly different profiles of local delta rhythmnested neuronal excitation (EPSPs) and response (AP generation) in all three neuronal subtypes tested (Fig. 6). IB neurons in layer 5 (the primary generators of the local neocortical delta rhythm) generated trains of action potentials overlying ramped excitatory postsynaptic potentials (EPSPs; Fig. 6D, Fig. 8). A small, but significant, decrease in outputs from these neurons was seen at higher cholinergic excitation (P Ͻ 0.05).
In contrast, for RS neurons in layers 5 and 2 and 3, both EPSP inputs and AP output probabilities became bimodal under higher cholinergic excitation, with dual EPSPs leading to dual, discrete periods of output generation (Fig. 6, B and C). These dual peaks in input and output were separated, on average, by one theta period within each accompanying delta period (L5 RS, 175 Ϯ 12 ms, L2/3 RS, 205 Ϯ 25 ms)-an observation similar to the peak separation in the MEG crossfrequency coupling measures (Fig. 4). MEG is thought to measure mean local dendritic current flow in populations of neurons throughout the cortical mantle. Therefore, we compared the model data with the subject data by pooling the mean EPSP profiles (Fig. 7B) and spike probabilities (Fig. 7C) from all three neuron subtypes and compared with the global mean delta-theta PAC profile (Fig. 7A). The transition from single to dual peaks was significant in each case (FDR, P Ͻ 0.05, Fig.  7A, ii-C, ii). In addition, no significant difference was seen between the corresponding model and human data cross-frequency coupling profiles. Those associated with "related" or "unrelated" semantic interpretations in subjects compared with model data in low-and high-cholinergic conditions, respectively (P Ͼ 0.05 for pooled mean EPSP and AP probability measures). This suggested that reaching an "unrelated" interpretation for paired words involved cholinergically mediated, dual-state changes at theta frequency in RS neuron inputs/ outputs temporally aligned to activity in delta-generating IB neurons (see DISCUSSION).

DISCUSSION
The data presented here demonstrate that theta and delta frequency dynamics dominate the spectral profile of cortical activity during semantic processing. Temporal coordination at theta frequency exposed the core regions of the semantic cognition network revealed by fMRI studies (see INTRODUC-TION). Additional temporally coupled regions exposed were the task-relevant primary sensory regions (bilateral V1), the hippocampi (but with the caveat that signal from deep structures can be unreliable), and key control regions involved in decision making (MFG, vmPFC). The profile of delta-theta PAC within this network provided a predictor of semantic decision outcome in the word-pair task used here, with PAC profile in the period from stimulus presentation to RT showing dual peaks for unrelated decisions and single peaks for related decisions. Dual peaks in delta-theta PAC were reproduced in a biological model in a manner dependent on cholinergic neuromodulation, Fig. 4. The profile of delta/theta phase amplitude coupling predicted semantic interpretation. Phase-amplitude coupling of region-averaged data, filtered for theta magnitude (4.5-7.5 Hz), with reference to the mean delta rhythm phase. The cross-frequency coupling profiles were significantly different (P Ͻ 0.05, n ϭ 17, Kologmorov-Smirnov) when comparing activity before a "related" (black) versus an "unrelated" (red) interpretation. Inset: significant phaseamplitude coupling (PAC) epochs polar plotted by delta phase.
suggesting multiple activations of cholinergically excited theta-generating neurons were required to reach an "unrelated" compared with a "related" semantic decision. How could such a temporal pattern arise, and how does it fit with current theories of the role of cross-frequency coupling in cognition?
Cross-frequency coordination of brain rhythms is a potent substrate for a broad range of cognitively relevant processes (Hyafil et al. 2015), and the pattern of coordination may represent a range of different neuronal network architectures if it is assumed the phenomenon is biologically generated. Unfortunately, apparent interactions can be generated nonbiologically in many different situations related to data type, quality, analysis methods, and external inputs (Aru et al. 2015). However, several observations do suggest the present observations were a consequence of a biological substrate: 1. The phase reset that generated the theta and delta average responses seen occurred on presentation of the first word not the second. This may explain why ERP measures alone exposed only the regions primarily receiving the sensory information (Fig. 2C) rather than those involved in the semantic task-insufficient information being presented at the time of phase-reset to perform any semantic computation. 2. The second word presentation was not time-locked to the intrinsic frequency of the delta rhythms or its nested theta rhythm (-350 ms, 2.9 Hz) stimulus separation (refer to 2.5 Hz and 5.2 Hz, respectively). 3. The results presented were from all the delta/theta periods occurring from stimulus presentation to decision point. Pooling multiple periods of potentially cross-frequency-coupled rhythms over time minimizes any effects of time-discrete inputs (Aru et al. 2015). 4. Modeling the observed pattern of cross-frequency coupling reproduced the results from human data with high fidelity in the absence of any time-discrete external inputs (Fig. 7). In addition to using the model to help ablate any possible artifacts from sensory inputs contributing to the results, it also allowed us to speculate on the biological origins of the crossfrequency interactions seen. Despite the species and temperature differences, the model reproduced the same frequency range of theta seen in MEG data and slightly slower delta frequency. Pharmacological generation of theta/delta rhythms in this manner has been shown previously to work equally well in human and rodent tissue in vitro (Carracedo et al. 2013), and the data presented here suggest that the persistent rhythm in the rat can be phase-reset in a similar manner to that seen in the human MEG data (Fig. 5). The consequent stimulus-locked delta rhythm survived for a longer time in the model than in subjects, but this may be a consequence of the absence of Stimuli were presented at the times indicated by the arrows. Middle: single example of raw data (black) and the corresponding filtered data (blue) to illustrate the continual presence of delta activity pretrial and posttrial. Bottom: phase for each of the 92 signals (black) and the average phase (red). B: spectrograms were constructed from A in two different ways: Upper spectrogram shows the power in the broadband signal derived from the average of each single trial, while the lower spectrogram shows the power from the averaged signal. Note the stimulus-induced delta power increase is only apparent from the average signal, the actual delta activity is persistent. C: in vitro model of persistent delta rhythms reproduced the stimulus-induced phase reset of persistent delta rhythms. Top: signal from parietal cortex aligned to electrical stimulus is shown, with stimulation of thalamus occurring in thalamocortical slices. Individual signals are shown in black, whereas the stimulus-aligned average is shown in red. Lower figure shows the corresponding phase of the recorded oscillation. D: phase reset plot for thalamic input to parietal cortex in the model. Note the highly linear, "near-absolute" nature of the phase reset.
further stimuli related to reaching the decision point that may be present in the MEG trial data.
The overall change in PAC profile was remarkably similar in the model and human task-dependent data. In this respect, the model predicted a highly discrete division of labor across the two, coupled spectral components dominating the MEG signals. Theta rhythmicity, both in terms of excitatory input and action potential output, was the sole property of regularly spiking (RS) neocortical neurons. These neurons are purely cortico-cortical in terms of their long-range projections; that is to say, they communicate exclusively with their immediate neighbors and distal cortical regions directly. In contrast, delta rhythmicity arises purely from layer 5 intrinsic bursting neurons (Carracedo et al. 2013), and this principal neuron subtype projects to immediate neighbors and exclusively to distal subcortical regions (Groh et al. 2010;Kim et al. 2015). In addition, the local mechanisms of generation of these two rhythms are very different; the delta rhythm is generated in the thalamus and neocortex by synaptically coupled bursting neurons operating in networks controlled by GABA B receptormediated synaptic inhibition (Carracedo et al. 2013;Hughes et al. 1998), whereas the theta rhythm (5-8 Hz) is generated by complex interactions between fast and slow GABA A receptormediated inhibitory synaptic inputs to principal cells (Gillies et al. 2002).
This division of labor has at least two consequences for the dynamical differences shown here before reaching a semantic decision. First, the observed bimodal nature of deep and superficial layer RS neuronal output probabilities, separated by one theta period, would suggest a role in reinforcing interaction between local cortical layers. Reciprocal interactions between deep and superficial layers have been proposed to represent a biological substrate for unsupervised learning and recall in layered neural networks (Carracedo et al. 2013;Dayan et al. 1995). Longer-range cortico-cortical interactions at theta frequency are also suggested to underlie error detection with reference to memory (Jacobs et al. 2006;Jensen and Tesche 2002), and in prefrontal, medial frontal, and parietal areas (Luu et al. 2004;Trujillo and Allen 2007;van Driel et al. 2012). These regions formed part of the task-related circuit revealed by theta-related region coupling (Fig. 3D). This suggested a commonality between error signaling and the "actual versus expected" nature of the comparison between 1st and 2nd presented words in the task. The overt "state-change"-like nature of the output probabilities in RS neurons resembled activity patterns shown in orbitofrontal neurons during subjective decision making previously (Rich and Wallis 2016) (Fig. 8).
Second, the theta-frequency state changes seen in RS neurons were strongly phase-related to the delta rhythm manifest in IB neurons. Outputs from this cell subtype target subcortical structures, such as colliculus, striatum, and thalamus, all of which have been shown to be involved in decision making. The ramped inputs to IB neurons during delta in the model presented here closely resembled those inferred from existing models of neuronal dynamics during decision making (see INTRODUCTION, Hunt et al. 2012;Fig. 8). However, in these models, the ascending ramp of magnitude of excitation to "decision" arises predominantly from the sensory input containing the "values" used for choice. How can this be replicated in a model that has no sensory input (the isolated parietal cortex preparation)? The nature of the model indicated that "ramps" in excitatory synaptic input were generated locally within neocortex via feedforward interactions between slow (predominantly NMDA receptor-mediated) recurrent connections between adjacent IB neurons. This fits the known profile of excitatory connections to neocortical neurons. Only~5% of inputs to neocortical cells arise from ascending cortical input, even in primary sensory areas (Peters and Payne 1993), and it is considered extremely rare to find a neuron that actually Fig. 6. Biological model of concurrent theta/delta activity reproduces magnetoencephalography (MEG) outcome-specific dynamics via altered cholinergic excitation. A: schematic of cells targeted in the in vitro parietal cortex slice. Regularly spiking (RS) neurons were recorded in layers 2 and 3 and layer 5. Intrinsically bursting (IB) neurons were recorded in layer 5. Spontaneous delta activity was recorded in the presence of cholinergic agonist (carbachol) at two different concentrations (2 M, black lines, 6 M, red lines). B: means Ϯ SE probability of action potential generation in layers 2 and 3 RS neurons binned according to phase of layer 5 field potential delta rhythm. *Significantly larger action potential probabilities (P Ͻ 0.05; n ϭ 10 replicates per n ϭ 5 slices/neurons) were seen for the higher level of cholinergic excitation from 30°to 90°. Example traces show activity at resting membrane potential (spikes, upper examples) and at Ϫ70 mV [excitatory postsynaptic potentials (EPSPs), lower examples]. C: similar profile of action potential output changes was seen when comparing outputs from layer 5 RS neurons. D: action potential bursts dominated layer 5 IB neuron activity. Elevated cholinergic excitation significantly reduced these burst outputs (*P Ͻ 0.05; n ϭ 10 replicates per n ϭ 5 slices/neurons). Fig. 7. Neuronal EPSP (input) and action potential generation (output) profiles match theta/delta crossfrequency coupling profile differences for "related" versus "unrelated" semantic interpretations. Decision-dependent differences in magnetoencephalography (MEG) theta-delta cross-frequency coupling (A) were not significantly different from the mean excitatory postsynaptic potential (EPSP) profiles (B) and action potential (C) profiles pooled from all three neuron subtypes (P Ͼ 0.05, Kologmorov-Smirnov, MEG "related" versus model lower cholinergic excitation (black), and MEG "unrelated" versus model higher cholinergic excitation (red). Lower graphs show polar plots for values in the upper-quartile range of the distribution. Note the presence of an additional peak later in the delta period when comparing the two outcomes and model conditions. Fig. 8. Biological model predicts theta rhythm involvement in cortico-cortical communication and delta rhythm involvement in cortical-subcortical communication. Top: layer 5 intrinsically bursting (IB) neuronal behavior in the experimental model of coupled theta and delta activity displayed overt ramplike excitatory postsynaptic potentials (EPSPs), indicative of increasing synaptic drive during the active part of one delta rhythm duty cycle. This behavior resembled that of accumulator neurons used to model decision making. These neurons are all subcortically projecting in the neocortex (Groh et al. 2010;Kim et al. 2015). Accumulator model figure was adapted from Hunt et al. 2012. Bottom: layer 5 and layers 2 and 3 regularly spiking (RS) neuronal behavior in the experimental model showed state-like changes in action potential output probabilities driven by multiple EPSP inputs in a manner dependent on subsequent decision. Similar changes in the probability of neuronal activity were reported for frontal cortical neurons using linear discriminant analysis (LDA) of unit activity. LDA probability figure adapted from Rich and Wallis (2016) with permission from Springer Nature [Nat Neurosci, Decoding subjective decisions from orbitofrontal cortex, Rich EL and Wallis JD, 2016]. Regularly spiking neurons are exclusively cortico-cortical projecting neurons. responds to a single sensory stimulus or stimulus property (Abeles 1988). The vast majority of inputs arise from local and distal cortical efferents (Young 2000). Thus, the "ramp" of inputs is best considered an intrinsic property of cortical networks, with the cortical regional profile of sensory input dictating the onset and gradient of the "ramp"-the degree and spatial extent of phase reset of the ongoing delta rhythm (e.g., Fig. 3A).
The strong correlation between semantic decision outcome and the profile of theta-delta cross-frequency coupling suggested that the above two phenomena need to interact to perform this cognitive task. The model used here suggested that cholinergic neuromodulation was critical in shaping this interaction. While cholinergic responses are the rule rather than the exception in all neocortical neuron subtypes, the converse-neocortical neurons influencing cholinergic neuronal excitability-is very rare. Interestingly, the one neocortical region shown to have a relatively strong influence on the magnitude of cholinergic neuronal activity is the orbitofrontal cortex (Do et al. 2016). This suggests a feedforward relationship between "state-change" activity [i.e., periodic, theta frequency alterations in output probability (Rich and Wallis 2016)] in this area and the extent of cholinergic drive mediating theta frequency state changes in the wider neocortex. Thus, modulation of the model local circuit via cholinergic excitation may not be solely related to attention and stimulus novelty (Acquas et al. 1996), but also to the ongoing activity associated with semantic processing within the neocortex.
We suggest that state change (Rich and Wallis 2016) and accumulator models (Wang 2002) may, therefore, reflect two temporally coherent, exquisitely self-governing dynamic signatures of the same local neocortical circuits, synergistically modifying both local and brain-wide neuronal activity before semantic decision. In the present data, these dynamic signatures were manifested as theta and delta frequency rhythms arising from different neuronal subtypes. The projection profile of these neurons suggested that cross-frequency coupling may not be just facilitating simultaneous interactions over different spatial scales, but rather coordinating cortico-cortical and cortico-subcortical interactions. These rhythms represented the majority of the spectral content of the MEG data between stimulus presentation and decision time point. However, they also occur coupled to faster frequencies (e.g., Lakatos et al. 2005) that have also been linked to semantic processing in key brain areas (Mollo et al. 2017). More comprehensive studies, over a broader dynamic range will, therefore, likely reveal much richer detail of the semantic decision process and the role of cross-frequency coupling in brain-wide communication, in general.

GRANTS
This article was supported by the Wellcome Trust, IBM, NIH/NINDS R01NS044133, and the ERC (SEMBIND, FLEXSEM).

DISCLOSURES
No conflicts of interest, financial or otherwise, are declared by the authors.