The Impact of Cortical Lesions on Thalamo-Cortical Network Dynamics after Acute Ischaemic Stroke: A Combined Experimental and Theoretical Study

Abstract

The neocortex and thalamus provide a core substrate for perception, cognition, and action, and are interconnected through different direct and indirect pathways that maintain specific dynamics associated with functional states including wakefulness and sleep. It has been shown that a lack of excitation, or enhanced subcortical inhibition, can disrupt this system and drive thalamic nuclei into an attractor state of low-frequency bursting and further entrainment of thalamo-cortical circuits, also called thalamo-cortical dysrhythmia (TCD). The question remains however whether similar TCD-like phenomena can arise with a cortical origin. For instance, in stroke, a cortical lesion could disrupt thalamo-cortical interactions through an attenuation of the excitatory drive onto the thalamus, creating an imbalance between excitation and inhibition that can lead to a state of TCD. Here we tested this hypothesis by comparing the resting-state EEG recordings of acute ischaemic stroke patients (N = 21) with those of healthy, age-matched control-subjects (N = 17). We observed that these patients displayed the hallmarks of TCD: a characteristic downward shift of dominant [alpha] -peaks in the EEG power spectra, together with increased power over the lower frequencies ([delta] and [theta]

-range). Contrary to general observations in TCD, the patients also displayed a broad reduction in [Beta]-band activity. In order to explain the genesis of this stroke-induced TCD, we developed a biologically constrained model of a general thalamo-cortical module, allowing us to identify the specific cellular and network mechanisms involved. Our model showed that a lesion in the cortical component leads to sustained cell membrane hyperpolarization in the corresponding thalamic relay neurons, that in turn leads to the de-inactivation of voltage-gated T-type Ca2+-channels, switching neurons from tonic spiking to a pathological bursting regime. This thalamic bursting synchronises activity on a population level through divergent intrathalamic circuits, and entrains thalamo-cortical pathways by means of propagating low-frequency oscillations beyond the restricted region of the lesion. Hence, pathological stroke-induced thalamo-cortical dynamics can be the source of diaschisis, and account for the dissociation between lesion location and non-specific symptoms of stroke such as neuropathic pain and hemispatial neglect.

Abstract

The neocortex and thalamus provide a core substrate for perception, cognition, and action, and are interconnected through different direct and indirect pathways that maintain specific dynamics associated with functional states including wakefulness and sleep. It has been shown that a lack of excitation, or enhanced subcortical inhibition, can disrupt this system and drive thalamic nuclei into an attractor state of low-frequency bursting and further entrainment of thalamo-cortical circuits, also called thalamo-cortical dysrhythmia (TCD). The question remains however whether similar TCD-like phenomena can arise with a cortical origin. For instance, in stroke, a cortical lesion could disrupt thalamo-cortical interactions through an attenuation of the excitatory drive onto the thalamus, creating an imbalance between excitation and inhibition that can lead to a state of TCD. Here we tested this hypothesis by comparing the resting-state EEG recordings of acute ischaemic stroke patients (N = 21) with those of healthy, age-matched control-subjects (N = 17). We observed that these patients displayed the hallmarks of TCD: a characteristic downward shift of dominant [alpha] -peaks in the EEG power spectra, together with increased power over the lower frequencies ([delta] and [theta]

-range). Contrary to general observations in TCD, the patients also displayed a broad reduction in [Beta]-band activity. In order to explain the genesis of this stroke-induced TCD, we developed a biologically constrained model of a general thalamo-cortical module, allowing us to identify the specific cellular and network mechanisms involved. Our model showed that a lesion in the cortical component leads to sustained cell membrane hyperpolarization in the corresponding thalamic relay neurons, that in turn leads to the de-inactivation of voltage-gated T-type Ca2+-channels, switching neurons from tonic spiking to a pathological bursting regime. This thalamic bursting synchronises activity on a population level through divergent intrathalamic circuits, and entrains thalamo-cortical pathways by means of propagating low-frequency oscillations beyond the restricted region of the lesion. Hence, pathological stroke-induced thalamo-cortical dynamics can be the source of diaschisis, and account for the dissociation between lesion location and non-specific symptoms of stroke such as neuropathic pain and hemispatial neglect.

Author Summary

The thalamus is involved in the relay and processing of most sensory information, and provides an interface between subcortical structures and the neocortex. However, disruptions in the subcortical communication with the thalamus are known to lead to thalamo-cortical dysrhythmia (TCD), which is linked to symptoms in a range of illnesses including Parkinson's disease, neurogenic pain syndrome and tinnitus. Thus far, TCD has solely been interpreted in terms of changes within subcortical pathways, but here we investigate how cortical disturbances (i.e., ischaemic stroke) may affect thalamic function in a similar manner. We do so by analysing the electroencephalogram (EEG) of stroke patients with a cortical lesion, and show that their EEG power spectra display the characteristic features of TCD. We subsequently built a detailed spiking model of thalamo-cortical circuits to identify the local cellular, circuit, and network properties and dynamics that lead to the development of this stroke-induced TCD. Together, our results shed light on less-understood symptoms of stroke such as neuropathic pain and hemispatial neglect, help inform future brain monitoring and diagnostics post-stroke, and suggest potential new treatments for stroke and related neurological conditions.

Citation: van Wijngaarden JBG, Zucca R, Finnigan S, Verschure PFMJ (2016) The Impact of Cortical Lesions on Thalamo-Cortical Network Dynamics after Acute Ischaemic Stroke: A Combined Experimental and Theoretical Study. PLoS Comput Biol 12(8): e1005048. doi:10.1371/journal.pcbi.1005048

Editor: Saad Jbabdi, Oxford University, UNITED KINGDOM

Received: March 24, 2016; Accepted: July 5, 2016; Published: August 10, 2016

Copyright: © 2016 van Wijngaarden et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Data Availability: All model files, including parameter sets and simulation results, are available from http://opendata.specs-lab.com/tcd_model/tcdModel.zip

Funding: This work was supported by the European Research Council under the European Union's Seventh Framework Programme FP7/2007-2013/ERC grant agreement n. 341196 [CDAC], and the Research and Technological Development Programme FP7-ICT-612139 in the WYSIWYD project to PFMJV. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing interests: The authors have declared that no competing interests exist.

Introduction

The brain is a complex network of segregated, functionally specialised, yet densely interconnected regions that exchange and integrate signals with high spatio-temporal precision [1,2]. Pathological perturbations of a focal area within this network can spread to distant regions, and affect their function via maladaptive processes like diaschisis, that leads to dysfunction in a (otherwise seemingly healthy) brain region due to its connectivity to a distant, damaged brain region [3,4]. One of the most common focal perturbations of neuronal tissue are regionally localised lesions due to stroke. With ischaemic stroke (IS), blood supply to the brain is blocked by the occlusion of one or more cerebral arteries, resulting in a cascade of pathophysiological events that eventually lead to neuronal cell death [5]. In the early minutes to hours after an ischaemic episode, most clinical symptoms result from functional impairment within the infarcted (dying) "core" region, and the surrounding "penumbra", an area also affected by ischaemia but potentially salvageable if blood flow is restored in a timely manner. In case penumbral regions are reperfused (e.g., by treatments delivered within hours of stroke onset), the symptomatology becomes more stable and specific, reflective of the loss of function within infarcted brain tissue [6,7]. Such symptoms can include sensory-motor deficits (i.e., sudden unilateral numbness or weakness in face or arm muscles), language disorders (i.e., aphasia), and/or cognitive deficits [8]. However, patients often suffer additional indirect and non-specific symptoms such as post-stroke pain and fatigue, hemispatial neglect, and mood-related disorders, that show an apparent dissociation with the lesion location, each with its own particular temporal dynamics in on- and off-set [9-12]. The mechanistic origins of many of these indirect symptoms of stroke are not well understood, and interventions remain undefined. Here we test whether stroke-induced pathological changes to the thalamo-cortical system (TCS) could serve as an underlying mechanism for generating such symptoms following an ischaemic episode.

The TCS displays a basic recurrent structure, reciprocally connecting thalamic nuclei with neocortical areas [13-15]. Disruptions within these circuits have been linked to symptoms in several neurological illnesses including Parkinson's disease, neurogenic pain syndrome, major depressive disorder, and tinnitus [16-19]. Associated symptoms are believed to develop as the result of an excessive expression of low-frequency bursting in thalamic nuclei and their propagation through the neocortex. This alteration of thalamo-cortical interactions is called thalamo-cortical dysrhythmia (TCD) [20]. TCD is caused by an imbalance of glutamatergic, cholinergic and/or GABAergic inputs to thalamic nuclei, creating hyperpolarised conditions that de-inactivate voltage-gated T-type Ca2+-channels, driving thalamic neurons into a low-frequency bursting mode of low-threshold calcium spikes (LTS) [21,22]. Synchronised bursting on a population level produces abnormal pathological oscillations in the [theta]

-range that propagate through TC efferents, influencing their afferent regions in the neocortex. This low-frequency rhythmic activation of cortical areas creates local asymmetry of lateral cortical inhibition, also known as the "edge-effect", which follows the characteristic shape of a Ricker wavelet, eliciting hyperactivity in neighbouring cortical modules [20,23]. Such hyperactivity, through increased coherence and cross-frequency coupling in the [theta]

and [Beta]-range, has been linked to symptoms including Parkinsonian tremors and sensations of pain [24,25]. This link is further supported by the use of interventions that target thalamic nuclei and their afferent structures, such as deep brain stimulation (DBS), where repolarisation of pathological thalamic neurons can reinstate normal neurophysiological function [26]. Alternatively, surgical removal of the affected thalamic nucleus (central lateral thalamotomy) is associated with 70-95% pain relief for patients with neurogenic pain syndrome [25,27].

Thus far, TCD has solely been interpreted in terms of changes to the subcortical drive onto thalamic nuclei. Our hypothesis is that in ischaemic stroke, cortico-thalamic network dynamics are disturbed; in particular, the neocortical drive onto the thalamus is attenuated, leading to the emergence of TCD. We test this hypothesis by comparing the resting-state electroencephalograms (EEG) of ischaemic stroke patients (N = 21) with those of healthy age-matched controls (N = 17). We show that the EEG of stroke patients displays the characteristic features of TCD relative to controls, and we explain the genesis of this phenomenon with a computational model of the TCS. Using our spiking model, we demonstrate that lesions in cortical circuits resulting from stroke lead to excessive hyperpolarisation of thalamic neurons, switching them to a LTS bursting regime by de-inactivating IT-currents. This switch induces TCD dynamics and entrains thalamo-cortical pathways, propagating low-frequency oscillations into the neocortex, and through the edge-effect, could serve as a possible mechanism underlying the diverse symptomatology found after ischaemic stroke.

Results

EEG results

Power spectra idiosyncratic to TCD display a slowing and occasional amplification of the dominant [alpha] -peak, increased power in the range of 1 to 30 Hz, and increased thalamo-cortical coherence [20,28]. Our stroke patients' EEG power spectra were consistent with this signature of TCD (Fig 1A). They displayed a significant increase in spectral power over the lower frequencies (Fig 1B) within the [delta]-band (p = .005) and [theta]

-band (p = .027) compared to healthy matched controls. Within the higher frequency ranges, we observed no significant increases in the [alpha] -band (p = .078), while [Beta]-band power was significantly attenuated (p < .001). This latter finding is opposite from the EEG characteristics generally reported for TCD, suggesting a pathology-specific modulation of TCD in the [Beta]-band. The dominant [alpha] -peaks also shifted significantly (p < .001) towards a lower frequency for patients (7.9±0.51 Hz, range 6.1-10.4 Hz) compared to controls (9.7±0.56 Hz, range 8.3-11.9 Hz), whereas these peaks did not significantly differ (p = .378) in power (Fig 1C).

Fig 1. EEG power spectra.

Wavelet convolution power spectra of the resting-state EEG data averaged over time and electrodes. (A) Mean spectra for patients (red) and controls (green), with the shaded area representing the 95%-confidence intervals. (B) Mean spectral energy (MSE) values (±s.e.m.) per frequency band, with significant differences in the [delta] (1-4 Hz), [theta]

(4-8 Hz), and [Beta] (12-30 Hz) band. (C) Locations of the individual dominant [alpha] -peaks, showing a significant slowing into a lower frequency for patients (7.9±0.3 Hz) as compared to controls (9.7±0.3 Hz). Statistically significant differences marked as **p < .01, *p < .05, Wilcoxon rank-sum test.

http://dx.doi.org/10.1371/journal.pcbi.1005048.g001

The analyses thus far spanned the average over all electrodes, but by arranging the spectral power differences over all frequencies against electrode positions, the topological distribution of differences between stroke patients and control subjects becomes evident (Fig 2A). This distribution in mean spectral energy (MSE) for the four different bands revealed localised increases in the [delta] and [theta]

-band, predominantly at the lateral ipsi-lesional electrodes (F8, C4, T4 and T6), and to a lesser extent contra-lesional (T3), together with the occipital electrodes (O1 and O2) (Fig 2B and 2C). The significant decrease of [Beta]-band power for patients was not localised, but was instead distributed across all scalp electrodes, while there were no differences between patients and controls in the [alpha] -band (Fig 2D and 2E). Individual patient comparisons with the control group average can be found in the SI (S1 Fig).

Fig 2. Topographical distribution of average differences between patients minus controls.

Patients suffered a lesion in either of the two hemispheres, and in order to prevent averaging out hemispheric differences, all non-midline electrodes were swapped laterally for patients with a left-hemispheric lesion. Thus, the right hemisphere in the diagram should be considered ipsi-lesional while the left hemisphere is contra-lesional. (A) Average spectral power differences between patients minus controls, organised per electrode and marked with an I for ipsi-lesional and C for contra-lesional, with clear boundaries separating the different bands around 8 and 12 Hz. (B-E) Topographical distribution of binned power differences for the four separate frequency bands: [delta], [theta]

, [alpha] and [Beta] respectively. Differences in the [delta] and [theta]

-band are qualitatively similar, where patients display an increase over the lateral ipsi-lesional electrodes (F6, C4, T4 and T6) together with the contra-lesional electrode T3 and occipital electrodes O1 and O2. There are no differences in the [alpha] -band, and the decrease of [Beta]-band activity is spread across all scalp electrodes.

http://dx.doi.org/10.1371/journal.pcbi.1005048.g002

Modelling thalamo-cortical dynamics post-stroke

In order to identify the mechanistic origins of stroke-induced TCD, we built a spiking model of the TCS, based on our previous work [29], and calibrated with the EEG data (see methods for further details). Our model captures the reciprocal connectivity between a cortical area (CRX), a modality-specific thalamic relay nucleus (SP), a multi-modal and non-specific thalamic nucleus (NSP), and the inhibitory thalamic reticular nucleus (TRN) (Fig 3). With this model we wanted to: (i) show that cortical lesions produce TCD-like dynamics post-ictal within the thalamus, as observed after peripheral deafferentation in modelling neurogenic pain [29]; (ii) explain the downward shift of dominant [alpha] -peaks in the cortical power spectra; and (iii) explain the stroke-specific decrease of spectral power within the [Beta]-band.

Fig 3. Schematic representation of the TCS model architecture.

The thalamic component of the model comprises specific (SP; corresponding to the ventral posterior complex of the thalamus in the somato-sensory domain), non-specific (NSP; comparable to the intralaminar nuclei) and the thalamic reticular nucleus (TRN). The cortical layer includes both excitatory (Exc.) pyramidal cells and inhibitory (Inh.) interneurons. Excitatory connections are marked with red lines terminating in arrowheads, while inhibitory connections are blue with circles indicating their synaptic terminals. The simulated stroke is implemented by removing the cortico-cortical connections from a specific cortical area, marked in black. All the model's parameter values can be found in Tables 1-3.

http://dx.doi.org/10.1371/journal.pcbi.1005048.g003

Table 1. Network connectivity values, based on [29].

http://dx.doi.org/10.1371/journal.pcbi.1005048.t001

Table 2. Thalamic neuron parameter values.

http://dx.doi.org/10.1371/journal.pcbi.1005048.t002

Table 3. Cortical neuron parameter values.

http://dx.doi.org/10.1371/journal.pcbi.1005048.t003

During a state of wakefulness, most thalamic neurons function under relatively depolarised conditions and, with inactive IT-currents, respond to excitatory input with sustained firing of unitary spikes [30]. These dynamics are captured in the model pre-lesion, with most neurons in a tonic firing mode and little to no coherence across the nuclei (Fig 4A-4D). We then approximated the conditions of stroke by lesioning thirty percent of the cortical population within the model, partially removing driving inputs to the thalamic nuclei. We observed that these conditions result in hyperpolarisation relative to baseline, where T-type Ca2+-channels are de-inactivated and neuronal discharges become LTS bursts rather than single spikes (Fig 4B-4D). As the recurrent cortical projections predominantly terminate onto the non-specific and reticular nuclei, bursting is initially confined to a subcircuit within the NSP and TRN. However, bidirectional divergent connectivity between these nuclei promotes the propagation of bursting activity to non-affected regions. Once a critical mass of neurons enter a bursting regime, dynamic recruitment of the TRN leads to hyperpolarised conditions for an increasing number of neurons in the NSP. This alternating interplay between the NSP and TRN synchronises neural activity across the different nuclei in the [theta]

-range (Fig 4A). The SP however is driven largely by peripheral sensory inputs, and a cortical lesion attenuates excitatory input for this nucleus to a much lesser extent. These neurons remain functional under relatively depolarised conditions in a tonic firing mode, and only a minor subset of neurons display bursting behaviour (Fig 4D).

Fig 4. Model results.

(A) Rasterplot of spikes for a single run of the model. From top to bottom: 100 neurons of the three thalamic nuclei TRN, NSP, SP plus the excitatory cortical population (CRX). A lesion was introduced after 10 seconds of real-time simulation time by deafferentating 30% of the cortical neurons, resulting in a loss of driving inputs for the thalamic neurons marked in grey. Arrows indicate entrained oscillations across the different nuclei. (B-D) Inter-spike intervals (ISI) for the NSP, SP and TRN, respectively. Most thalamic neurons pre-lesion are in a tonic-firing mode, with evenly distributed ISIs. After the lesion, the NSP and TRN switch to a bursting regime, with the majority of neurons having an ISI below 20 ms indicating intra-burst spikes, followed by longer inter-burst intervals. (E) Power spectra for the NSP (solid line), SP (dotted line) and TRN (dashed line) in the lower frequency range. All thalamic nuclei display a substantial increase of spectral power in the low-frequency range, peaking between 6 and 8 Hz, with the strongest oscillations present in the NSP and TRN. (F) Peak frequencies of the cortical power spectra in the range of 6-13 Hz, with an overall downward shift into a lower frequency post-lesion (9.2±0.38 Hz) compared to pre-lesion (9.9±0.44 Hz). Line colours indicate the highest corresponding frequency. (G) Average power spectra of the cortical component with and without the extrasynaptic increase of GABA post-lesion, with shaded areas representing the 95%-confidence intervals. The lesion affects [gamma]-band power in both models post-lesion, but with an extrasynaptic increase of GABA, [Beta]-band power is further suppressed from 23 Hz onwards, comparable to the decrease for patients in their EEG. (H) Peak frequencies for different lesion sizes, averaged over 10 simulations per step. As the lesion size increases, the dominant [alpha] -peaks slow down towards a lower frequency post-lesion. Statistically significant differences marked as **p < .01, Wilcoxon-rank sum test.

http://dx.doi.org/10.1371/journal.pcbi.1005048.g004

Synchronised bursting is characterised on a population level by a substantial increase in low-frequency oscillations post-lesion for all thalamic nuclei, peaking between 6 and 8 Hz, most prevalent in the NSP and TRN (Fig 4E). With the propagation of such oscillations through TC efferents, dominant [alpha] -peaks of the cortical power spectrum are driven from 9.9±0.44 Hz towards a lower frequency of 9.2±0.38 Hz post-lesion (p = .002) (Fig 4F), coherent with slowing of the [alpha] -peaks found in the EEG spectra of patients. Parametrisation of the lesion size parameter shows its linear relationship with this shift, where an increased number of lesioned neurons correspond with a greater slowing of dominant peaks (Fig 4H).

During an ischaemic episode, the brain's response is to increase extrasynaptic levels of GABA as a neuroprotective mechanism, inhibiting all glutamate mediated neuronal activity in order to suppress rising levels of excitotoxicity, and prevent additional tissue damage [31]. To capture this extrasynaptic increase of GABA, a constant inhibitory current was applied to all cortical neurons post-lesion following a slow onset. The cortical lesion in the model resulted in a decrease of [gamma]-band activity, but adding this inhibitory current further suppressed [gamma]-band power (Fig 4G). Moreover, high-[Beta]-band activity also decreased, creating divergence between the pre- and post-lesion spectra from 23 Hz and upwards (Fig 4G).

Discussion

We have addressed the question whether thalamo-cortical dysrhythmia develops after ischaemic stroke. We hypothesised that cortical lesions, for instance due to ischaemic stroke, lead to an attenuation of the neocortical drive onto the thalamus, changing thalamic cellular and circuit dynamics such that thalamo-cortical dysrhythmia (TCD) emerges. Our analysis of the EEG of stroke patients confirms our hypothesis and shows the characteristic TCD shifts in the power spectrum. The electrophysiological data display high similarity with the well-established EEG characteristics found in the wide range of disorders previously linked with TCD [23,25,28].

The most prominent overlap is the slowing of the dominant [alpha] -peak, together with an increase in spectral power over the lower frequencies, both in the [delta]-band and [theta]

-band. These results are consistent with changes in the EEG power spectra generally observed in stroke. Rapid appearance of high-amplitude oscillations in the [delta]-band is highly typical of ischaemic stroke, and is often preceded by both a slowing of the [alpha] -peaks and attenuation of [Beta]-band oscillations in cases of intermediate to severe ischemia (see [32] for a review). This [Beta]-band attenuation, also observed in our patients, is opposite from what is generally reported for TCD, and could be related to increased levels of extrasynaptic GABA during the acute phase of stroke [31,33]. Indeed, approximating this increase in our computational model not only further suppressed [gamma]-band power, but also decreased high-[Beta]-band power from 23 Hz upwards, comparable to what we observed in the stroke patients.

Our patients' EEG were acquired during the acute phase, but are consistent with Dubovik et al. [34], who previously reported a similar downward shift in [alpha] -peaks, and comparable modulations in all frequency bands in EEG recorded 3 months post-stroke. The authors emphasise changes in the spatial distribution of [alpha] -band power, their coherence, and relationship to the symptomatology, supporting the notion of diaschisis following stroke [4]. Here we determine more precisely the similarities between such stroke-induced diaschisis and thalamo-cortical dysrhythmia, identifying a possible substrate underlying the dynamic changes to the TCS after stroke.

We subsequently modelled the dynamics of the TCS by which cortical insults lead to the genesis of pathological thalamic low-frequency oscillations and TCD-like phenomena. Using a spiking model of the thalamo-cortical system, we managed to identify core candidate mechanisms responsible for the attenuation found in the power spectra. Indeed, with a lesion in the cortex, excitatory inputs to thalamic nuclei are deafferentated, resulting in disfacilitation and cell membrane hyperpolarisation that in turn leads to the de-inactivation of T-type Ca2+-channels and a switch to LTS bursting. While initially confined to a subcircuit, activity spreads through divergent connectivity in the intrathalamic pathways, enforcing synchronised oscillations on a population level in the high-[theta]

range. The influx of such slow oscillations through thalamic efferents into the cortex drives the dominant [alpha] -peaks into a lower frequency. We predict that this entrainment of thalamo-cortical circuits will lead to further effects such as those observed in TCD.

While the propagation of low-frequency thalamic oscillations into the neocortex is common in a state of drowsiness or sleep, they have pathological properties in an awake state through an ''edge-effect'' [20,23]. Oscillations influence local interneurons and create a ring of imbalance in cortico-cortical inhibition surrounding the target area, resulting in waves of hyperactivity in neighbouring regions. During the acute phase of stroke, this mechanism could be involved in positive post-stroke symptoms that have a relatively short onset, such as post-stroke pain or fatigue, when similar pathways involved in neurogenic pain are entrained [23,25]. In contrast, TCD could additionally be associated with rapid emerging negative symptoms. For example, given the important role of the TRN in visual information processing and attention, a strong increase of reticular inhibition within the visual thalamo-cortical circuit could disrupt top-down inhibitory networks between prefrontal areas and thalamus, resulting in negative symptoms such as hemispatial neglect [11,35]. On a network level, prolonged pathological hyperactivity could affect the dynamics of cortical circuits during the chronic phase. While less understood, plasticity could lead to a dynamic reorganisation of cortical networks, resulting in chronic symptoms with a much longer onset, like depression, that typically develop after several months [18,36].

The approach to TCD has thus far been to examine subcortical drives onto the thalamus, either sensory input coming from peripheral nerve fibres (e.g., with neurogenic pain), or inhibitory projections from the basal ganglia (e.g., with Parkinson's disease) in the pallido-thalamic tract [37,38]. Our results suggest that TCD can extend into a much broader range of disorders, where any form of cortical tissue damage could lead to the development of TCD-like dynamics within the TCS. Cortical trauma can disrupt cortical function over spatially distributed areas through diaschisis, altering the physiological state of distant regions, and modulate their activity [34]. This process not only affects cortico-cortical connections, but also has a direct impact on thalamic function through cortico-thalamic pathways. While this change initially affects thalamic function locally, once it reaches the dynamical attractor state of TCD, pathological oscillations will persistently propagate through circuits of the neocortex, affecting increasingly more distal areas. Consistent with this notion is a recent link between TCD and severe sepsis and septic shock [39]. Sepsis is a severe inflammatory response to an infection, leading to sepsis-associated encephalopathy due to brain-cell damage, impaired neurotransmission, and neurodegeneration [40-42]. As a result, patients often suffer from a rapid decline in cognitive function, including impaired memory, attention, and visual-spatial processing. It has been proposed that the neuropathological processes following sepsis additionally lead to the development of TCD, as the patients' EEG display the distinctive slowing of the dominant [alpha] -peak into high-[theta]

range [39].

TCD is generally identified using either non-invasive EEG or MEG, a method that has been validated with intracranial recordings in patients undergoing stereotactic surgery [16,43,44]. While our current study does not include direct thalamic recordings, EEG allowed us to identify TCD post-stroke, but the exact relationship between TCD as a mechanism and the genesis of non-specific symptoms remains unclear, and thus requires further investigation. One direction is connecting extensive clinical assessment of stroke patients, covering cognitive domains such as neglect and attention, chronic pain, and mood-related impairments, with TCD-specific characteristic EEG features. Using high-density recordings, different source localisation techniques could furthermore be used to localise specific features such as [alpha] -peak shifts and low-frequency power modulations in relation to specific symptoms. Once this link between TCD and stroke-associated symptoms becomes more evident, it would open new avenues of interventions proven effective to treat TCD. This includes stereotactic thalamotomy or tractotomy, together with deep brain stimulation, used to treat Parkinson's disease and neurogenic pain syndrome [26,43,44].

Conclusions

By combining the analysis of electrophysiological (scalp EEG) recordings from stroke patients with a computational model of thalamo-cortical circuitry, we showed that characteristic, post-stroke EEG abnormalities can be accounted for in terms of TCD. We also identified the substrate of diaschisis and thalamo-cortical dysrhythmia, where stroke not only affects cortical function locally, but furthermore perturbs distributed cortical and thalamo-cortical networks (and thus brain regions distant to the lesion), making it a prime target for the emerging field of network medicine [45,46]. By understanding the adverse impact of stroke on brain network function in terms of TCD, we advance an alternative perspective on the genesis of diaschisis and the non-specific symptomatology of stroke. Similar disruptions within the TCS are present in other thalamo-cortical disorders, suggesting a common underlying mechanism among different neuropathologies, which can produce different symptoms depending on the affected thalamo-cortical circuit(s).

Methods

Clinical data

Approval to carry out the study was obtained from the local University and Hospital Human Research Ethics Committees. Written informed consent from each patient or substitute decision-maker was obtained. Participants for this study included stroke patients (N = 21, 11 female; mean age 72 years, range 38-85), all suffering from acute middle cerebral artery (MCA) stroke and recruited from the Royal Brisbane and Women's Hospital in Brisbane, Australia. Stroke was assessed using acute computed tomography (CT) scans, followed by magnetic resonance imaging (MRI) in six cases. EEG data were acquired at the patient's bedside in the acute phase, approximately 69 hours (range 21-99) after stroke onset and used previously to optimise electrode placement in stroke prognostics [47]. In case patients woke up with stroke symptoms, time of stroke onset was defined as the midpoint between bedtime and time of waking up. Patient demographic information can be found in the SI (S1 Table). EEGs were recorded using a NicOne Brain Monitor (Natus Medical Inc.), recording at a sampling rate of 500 Hz and using 19 Ag/Ag-Cl electrodes (Nicolet; Natus Medical Inc.), placed according to the international 10-20 system. The reference data came from healthy age-matched controls (N = 17, 8 female; mean age 68 years, range 60-80) with no cognitive impairments or history of depression and/or anxiety, previously gathered for [48]. EEGs were recorded using an elasticised quick cap with 32 Ag/Ag-CL electrodes (Neuromedical supplies), sampling at 500 Hz and digitised by a Neuroscan Synamps amplifier.

Topology model

A schematic of the model's architecture is shown in Fig 3. Thalamo-cortical dynamics were simulated using the open-source IQR neural network simulator [49] over a total of 40 separate simulations, each with a random initialisation of thalamo-cortical and cortico-cortical connectivity. The temporal integration of each time step was 1 ms with a total recording time of 20 seconds real-time. The model's network builds on [29] with similar thalamic circuitry, but we extended the Poissonian cortical input to a layer that includes both excitatory and inhibitory neurons and displays cortical dynamics (see neuron model below).

The thalamic component consists of three classes of nuclei: (i) a modality-specific relay nucleus (SP), corresponding to the ventral posterior complex of the thalamus in the somato-sensory domain; (ii) a non-specific higher-order nucleus (NSP), comparable to the intralaminar layer; and (iii) the inhibitory reticular (TRN) nucleus [30,50]. The excitatory SP and NSP nuclei are distinguished in terms of their distinct connectivity patterns within the model. The SP is characterised by parallel one-to-one connections (e.g., a single source neuron is connected to a single target neuron) with the TRN and is driven primarily by peripheral sensory inputs [51]. Conversely, the NSP receives most of its input from cortical layer VI and has strong divergent one-to-many connections (e.g., a single source neuron is connected to all target neurons with a given probability) with the TRN and cortex [52,53]. These anatomical properties are captured by the connectivity values between a source and target neurons of the model (Table 1). The random peripheral input for the SP and NSP consisted of Poisson spike trains with a given spiking probability P per iteration. We tuned the values of PSP = 0.5 and PNSP = 0.35 such that the nuclei at resting-state fire within a physiological range, while maintaining dominance of peripheral input for the SP [15].

Neuron model

Characteristic membrane properties of thalamic neurons include polarisation-dependent inactivation of T-type Ca2+-channels, capable of producing low-threshold calcium spikes under hyperpolarised conditions. The de-inactivation of such calcium currents is modelled by adding a slow variable, h, to the classical conductance-based leaky integrate-and-fire dynamics, suggested by [54]. TC neurons were modelled as single-compartment cells, where changes in the membrane potential depend on an input current, Iin, calcium currents, IT, and a constant conductance leak current IL: (1)

The calcium currents depend on the inactivation level, h, that relaxes to zero at depolarised levels when V > Vh, or approaches unity under hyperpolarised conditions, given time constants [tau]h- and [tau]h+: (2) (3)

The input current, Iin, is dependent on all excitatory minus inhibitory synaptic currents, gE/I, calculated as the sum over the multiplication of the connectivity matrix, Wij, with a dichotomous spiking vector, sj, that indicates spiking neurons (eqs 4 and 5). The result is weighted with a gain, kE/I, and follows exponential decay dynamics with time constants [tau]E/I. All parameter values can be found in Table 2 and were chosen based on [29,54].

(4)(5)

The cortical population includes 800 excitatory and 200 inhibitory neurons and were modelled as quadratic integrate-and-fire cells [55], according to a set of differential equations that capture the membrane potential, v, and a state variable, u: (6) (7)

A cell is considered to spike when the membrane potential reaches a threshold of 30 mV, after which v [arrow left] c and u [arrow left] u + d capturing after-spike repolarisation. Parameters a, b, c and d were chosen to produce firing behaviour similar to regular spiking (RS) neurons (Table 3), in concordance with [55]. We explored the parameter space for the intrinsic input, I, to find the optimal combination of values (IE = 6.7 and II = 2.7) that reproduces power spectra with dominant [alpha] -peaks comparable with those found in the healthy EEG data.

After simulating the thalamocortical system in a healthy state, a structural lesion was approximated by deafferentating 30% of adjacent excitatory and inhibitory cortical neurons, disconnecting all connections to and from these neurons. The impact of different lesion sizes on network behaviour can be found in Fig 4H. To retain comparability with the human EEG data, a noise constraint was applied by omitting simulations (4 out of 40) with dominant [alpha] -peaks pre-lesion outside boundaries given by the minimum and maximum peak frequencies of the control group's EEG.

Data analysis

Signal processing and analyses were performed offline using MatLab (Mathworks, Natrick, Ma, USA) and the toolbox EEGLAB [56], together with custom in-house scripts. All data were bandpass filtered between .5 and 35 Hz (12dB/octave) and re-referenced to the common average of all electrodes. The first 45 epochs of 2048 ms were selected that included no clear artefacts or extreme values of ±75 MV. Only the following 17 electrodes common across both groups were included for further analyses: F3, F4, F7, F8, Fz, C3, C4, Cz, P3, P4, Pz, T3, T4, T5, T6, O1 and O2.

Wavelet convolution was used to decompose both the EEG and model time-series data into their underlying oscillatory components [57]. At first, a set of complex Morlet wavelets were computed, scaling logarithmically from 1 to 35 Hz for the EEG and from 1 to 80 Hz for the model, normalised by the wavelet's maximum value. With an equal amount of cycles, higher frequencies span a shorter time-window than lower frequencies, hence the amount of cycles used scaled in accordance with their frequency. Next, the fast Fourier transformation (FFT) of the signal was multiplied in complex space with the FFT of each wavelet, and taking the inverse FFT gives an analytical signal in the time-domain. Taking the square of the real part of this signal provides the amplitude at each time-frequency point. Input for the convolution were pre-processed time-series signals of the EEG, described in detail in [47], while a local field potential (LFP) was used for the model, calculated as the combined membrane potentials for all non-lesioned neurons per population. The result of the convolution was averaged over time and binned into four commonly used frequency bands: delta ([delta]; 1.0-4.0 Hz), theta ([theta]

; 4.2-7.9 Hz), alpha ([alpha] ; 8.3-11.9 Hz) and beta ([Beta]; 12.4-30.6 Hz), giving the mean spectral energy (MSE) per band.

Statistical analysis

All statistical testing was done using the non-parametric Wilcoxon rank-sum test, with significant differences marked as **p < .01, *p < .05. Values reported in text are always means ±s.e.m. and shaded areas in figures represent the 95%-confidence intervals.

Supporting Information

S1 Table. Clinical description of patients.

(DOCX)

S1 Fig. Topographical distribution of individual differences in spectral power.

A difference was computed between each patient individually minus the control group average in the combined frequency range of [delta] and [theta]

(1-8 Hz). Stroke patients show high heterogeneity in affected regions. There are either unilateral or bilateral increases (i.e., patient 2 versus 3) that can be localized or spread across the scalp (i.e., patient 7 versus 12), affecting different cortical areas such as frontal regions (patient 14), visual areas (patient 2), motor areas (patient 20) or parietal regions (patient 12). Despite these individual differences, we do observe global changes across all patients, such as a shift in the [alpha] -peaks, suggesting the development of TCD.

(EPS)

Acknowledgments

The authors would like to thank X. Arsiwalla for useful discussions on the results, M. Boffo and K. Nikolaou for their comments on earlier versions of the manuscript, and S. Wierenga for assisting with designing the figures.

Author Contributions

1. Conceived and designed the experiments: JBGvW PFMJV RZ.

2. Performed the experiments: JBGvW RZ.

3. Analyzed the data: JBGvW.

4. Wrote the paper: JBGvW PFMJV SF RZ.

5. Clinical data collection: SF. Conceived hypothesis: PFMJV.

References

Bullmore E, Sporns O. Complex brain networks: graph theoretical analysis of structural and functional systems. Nat Rev Neurosci. 2009;10: 186-198. doi: 10.1038/nrn2575. pmid:19190637

Sporns O. Network attributes for segregation and integration in the human brain. Curr Opin Neurobiol. 2013;23: 162-171. doi: 10.1016/j.conb.2012.11.015. pmid:23294553

Fornito A, Zalesky A, Breakspear M. The connectomics of brain disorders. Nat Rev Neurosci. 2015;16: 159-172. doi: 10.1038/nrn3901. pmid:25697159

Carrera E, Tononi G. Diaschisis: past, present, future. Brain. 2014;137: 2408-2422. doi: 10.1093/brain/awu101. pmid:24871646

Brouns R, De Deyn PP. The complexity of neurobiological processes in acute ischemic stroke. Clin Neurol Neurosurg. 2009;111: 483-495. doi: 10.1016/j.clineuro.2009.04.001. pmid:19446389

Dirnagl U, Iadecola C, Moskowitz MA. Pathobiology of ischaemic stroke: an integrated view. Trends Neurosci. 1999;22: 391-397. doi: 10.1016/S0166-2236(99)01401-0. pmid:10441299

Lo EH, Dalkara T, Moskowitz MA. Neurological diseases: Mechanisms, challenges and opportunities in stroke. Nat Rev Neurosci. 2003;4: 399-414. doi: 10.1038/nrn1106. pmid:12728267

Ferro JM, Martins IP, Caeiro L. Behavioral neurology of stroke. In: Textbook of stroke medicine. Cambridge University Press; 2012. pp. 178-193

Appelros P. Prevalence and predictors of pain and fatigue after stroke: a population-based study. Int J Rehabil Res. 2006;29: 329-333. doi: 10.1097/MRR.0b013e328010c7b8. pmid:17106351

O'Donnell MJ, Diener H-C, Sacco RL, Panju AA, Vinisko R, Yusuf S. Chronic Pain Syndromes After Ischemic Stroke: PRoFESS Trial. Stroke. 2013;44: 1238-1243. doi: 10.1161/STROKEAHA.111.671008. pmid:23559265

Ringman JM, Saver JL, Woolson RF, Clarke WR, Adams HP. Frequency, risk factors, anatomy, and course of unilateral neglect in an acute stroke cohort. Neurology. 2004;63: 468-74. doi: 10.1212/01.WNL.0000133011.10689.CE. pmid:15304577

Robinson RG. Neuropsychiatric consequences of stroke. Annu Rev Med. 1997;48: 217-29. doi: 10.1146/annurev.med.48.1.217. pmid:9046957

Jones EG. The thalamus. 2nd ed. Cambridge University Press; 2007

Sherman SM, Guillery RW. Exploring the thalamus and its role in cortical function. 2nd ed. Cambridge, MA: MIT Press; 2006

Guillery RW, Sherman SM. Thalamic Relay Functions and Their Role in Corticocortical Communication. Neuron. 2002;33: 163-175. doi: 10.1016/S0896-6273(01)00582-7. pmid:11804565

Sarnthein J, Jeanmonod D. High Thalamocortical Theta Coherence in Patients with Parkinson's Disease. J Neurosci. 2007;27: 124-131. doi: 10.1523/JNEUROSCI.2411-06.2007. pmid:17202479

Sarnthein J, Jeanmonod D. High thalamocortical theta coherence in patients with neurogenic pain. Neuroimage. 2008;39: 1910-1917. doi: 10.1016/j.neuroimage.2007.10.019. pmid:18060808

Leuchter AF, Hunter AM, Krantz DE, Cook I a. Rhythms and blues: modulation of oscillatory synchrony and the mechanism of action of antidepressant treatments. Ann N Y Acad Sci. 2015;1344: 78-91. doi: 10.1111/nyas.12742. pmid:25809789

De Ridder D, Vanneste S, Langguth B, Llinás R. Thalamocortical Dysrhythmia: A Theoretical Update in Tinnitus. Front Neurol. 2015;6: 124. doi: 10.3389/fneur.2015.00124. pmid:26106362

Llinás RR, Ribary U, Jeanmonod D, Kronberg E, Mitra PP. Thalamocortical dysrhythmia: A neurological and neuropsychiatric syndrome characterized by magnetoencephalography. Proc Natl Acad Sci. 1999;96: 15222-15227. doi: 10.1073/pnas.96.26.15222. pmid:10611366

Steriade M. Impact of network activities on neuronal properties in corticothalamic systems. J Neurophysiol. 2001;86: 1-39. pmid:11431485

Llinás RR. Bursting of Thalamic Neurons and States of Vigilance. J Neurophysiol. 2006;95: 3297-3308. doi: 10.1152/jn.00166.2006. pmid:16554502

Llinás R, Urbano FJ, Leznik E, Ramírez RR, van Marle HJF. Rhythmic and dysrhythmic thalamocortical dynamics: GABA systems and the edge effect. Trends Neurosci. 2005;28: 325-333. doi: 10.1016/j.tins.2005.04.006. pmid:15927689

Jeanmonod D, Magnin M, Morel A. Low-threshold calcium spike bursts in the human thalamus. Brain. 1996;119: 363-375. doi: 10.1093/brain/119.2.363. pmid:8800933

Sarnthein J. Increased EEG power and slowed dominant frequency in patients with neurogenic pain. Brain. 2006;129: 55-64. doi: 10.1093/brain/awh631. pmid:16183660

Lozano AM, Lipsman N. Probing and Regulating Dysfunctional Circuits Using Deep Brain Stimulation. Neuron. 2013;77: 406-424. doi: 10.1016/j.neuron.2013.01.020. pmid:23395370

Stern J, Jeanmonod D, Sarnthein J. Persistent EEG overactivation in the cortical pain matrix of neurogenic pain patients. Neuroimage. 2006;31: 721-731. doi: 10.1016/j.neuroimage.2005.12.042. pmid:16527493

Schulman JJ, Cancro R, Lowe S, Lu F, Walton KD, Llinás RR. Imaging of Thalamocortical Dysrhythmia in Neuropsychiatry. Front Hum Neurosci. 2011;5: 69. doi: 10.3389/fnhum.2011.00069. pmid:21863138

Henning Proske J, Jeanmonod D, Verschure PFMJ. A computational model of thalamocortical dysrhythmia. Eur J Neurosci. 2011;33: 1281-1290. doi: 10.1111/j.1460-9568.2010.07588.x. pmid:21323765

Jones EG. Thalamic circuitry and thalamocortical synchrony. Philos Trans R Soc B Biol Sci. 2002;357: 1659-1673. doi: 10.1098/rstb.2002.1168. pmid:12626002

Green AR, Hainsworth AH, Jackson DM. GABA potentiation: a logical pharmacological approach for the treatment of acute ischaemic stroke. Neuropharmacology. 2000;39: 1483-1494. doi: 10.1016/S0028-3908(99)00233-6. pmid:10854894

Finnigan S, van Putten MJAM. EEG in ischaemic stroke: Quantitative EEG can uniquely inform (sub-)acute prognoses and clinical management. Clin Neurophysiol. 2013;124: 10-19. doi: 10.1016/j.clinph.2012.07.003. pmid:22858178

Clarkson AN, Huang BS, MacIsaac SE, Mody I, Carmichael ST. Reducing excessive GABA-mediated tonic inhibition promotes functional recovery after stroke. Nature. 2010;468: 305-309. doi: 10.1038/nature09511. pmid:21048709

Dubovik S, Pignat J-M, Ptak R, Aboulafia T, Allet L, Gillabert N, et al. The behavioral significance of coherent resting-state oscillations after stroke. Neuroimage. 2012;61: 249-257. doi: 10.1016/j.neuroimage.2012.03.024. pmid:22440653

Wimmer RD, Schmitt LI, Davidson TJ, Nakajima M, Deisseroth K, Halassa MM. Thalamic control of sensory selection in divided attention. Nature. 2015;526: 705-709. doi: 10.1038/nature15398. pmid:26503050

Gaete JM, Bogousslavsky J. Post-stroke depression. Expert Rev Neurother. 2008;8: 75-92. doi: 10.1586/14737175.8.1.75. pmid:18088202

Sarnthein J. Thalamic theta field potentials and EEG: high thalamocortical coherence in patients with neurogenic pain, epilepsy and movement disorders. Thalamus Relat Syst. 2003;2: 231-238. doi: 10.1016/S1472-9288(03)00021-9

Magnin M, Morel A, Jeanmonod D. Single-unit analysis of the pallidum, thalamus and subthalamic nucleus in parkinsonian patients. Neuroscience. 2000;96: 549-564. doi: 10.1016/S0306-4522(99)00583-7. pmid:10717435

Götz T, Baumbach P, Huonker R, Kranczioch C, Witte OW, Debener S, et al. Slowed peak resting frequency and MEG overactivation in survivors of severe sepsis and septic shock. Clin Neurophysiol. 2016;127: 1247-1253. doi: 10.1016/j.clinph.2015.07.021. pmid:26391681

Davies NWS, Sharief MK, Howard RS. Infection-associated encephalopathies--their investigation, diagnosis, and treatment. J Neurol. 2006;253: 833-845. doi: 10.1007/s00415-006-0092-4. pmid:16715200

Widmann CN, Heneka MT. Long-term cerebral consequences of sepsis. Lancet Neurol. 2014;13: 630-636. doi: 10.1016/S1474-4422(14)70017-1. pmid:24849863

Golzari SE, Mahmoodpoor A. Sepsis-associated encephalopathy versus sepsis-induced encephalopathy. Lancet Neurol. 2014;13: 967-968. doi: 10.1016/S1474-4422(14)70205-4. pmid:25231515

Jeanmonod D, Magnin M, Morel A, Siegemund M. Surgical control of the human thalamocortical dysrhythmia: I. Central lateral thalamotomy in neurogenic pain. Thalamus Relat Syst. 2001;1: 71-79. doi: 10.1017/S1472928801000036

Magnin M, Jeanmonod D, Morel A, Siegemund M. Surgical control of the human thalamocortical dysrhythmia: II. Pallidothalamic tractotomy in Parkinson's disease. Thalamus Relat Syst. 2001;1: 81-89. doi: 10.1017/S1472928801000024

Stam CJ. Modern network science of neurological disorders. Nat Rev Neurosci. 2014;15: 683-695. doi: 10.1038/nrn3801. pmid:25186238

Cramer AOJ, Waldorp LJ, van der Maas HLJ, Borsboom D. Comorbidity: A network perspective. Behav Brain Sci. 2010;33: 137-150. doi: 10.1017/S0140525X09991567. pmid:20584369

Schleiger E, Sheikh N, Rowland T, Wong A, Read S, Finnigan S. Frontal EEG delta/alpha ratio and screening for post-stroke cognitive deficits: The power of four electrodes. Int J Psychophysiol. 2014;94: 19-24. doi: 10.1016/j.ijpsycho.2014.06.012. pmid:24971913

Finnigan S, Wong A, Read S. Defining abnormal slow EEG activity in acute ischaemic stroke: Delta/alpha ratio as an optimal QEEG index. Clin Neurophysiol. 2016;127: 1452-1459. doi: 10.1016/j.clinph.2015.07.014. pmid:26251106

Bernardet U, Verschure PFMJ. iqr: A Tool for the Construction of Multi-level Simulations of Brain and Behaviour. Neuroinformatics. 2010;8: 113-134. doi: 10.1007/s12021-010-9069-7. pmid:20502987

Guillery RW. Anatomical evidence concerning the role of the thalamus in corticocortical communication: a brief review. J Anat. 1995;187: 583-92. pmid:8586557

Darian-Smith C, Tan A, Edwards S. Comparing thalamocortical and corticothalamic microstructure and spatial reciprocity in the macaque ventral posterolateral nucleus (VPLc) and medial pulvinar. J Comp Neurol. 1999;410: 211-234. doi: 10.1002/(SICI)1096-9861(19990726)410:2<211::AID-CNE4>3.0.CO;2-X. pmid:10414528

Van der Werf YD, Witter MP, Groenewegen HJ. The intralaminar and midline nuclei of the thalamus. Anatomical and functional evidence for participation in processes of arousal and awareness. Brain Res Rev. 2002;39: 107-140. doi: 10.1016/S0165-0173(02)00181-9. pmid:12423763

Saalmann YB. Intralaminar and medial thalamic influence on cortical synchrony, information transmission and cognition. Front Syst Neurosci. 2014;8: 83. doi: 10.3389/fnsys.2014.00083. pmid:24847225

Smith GD, Cox CL, Sherman SM, Rinzel J. Fourier analysis of sinusoidally driven thalamocortical relay neurons and a minimal integrate-and-fire-or-burst model. J Neurophysiol. 2000;83: 588-610. pmid:10634897

Izhikevich EM. Simple model of spiking neurons. IEEE Trans Neural Networks. 2003;14: 1569-1572. doi: 10.1109/TNN.2003.820440. pmid:18244602

Delorme A, Makeig S. EEGLAB: an open source toolbox for analysis of single-trial EEG dynamics including independent component analysis. J Neurosci Methods. 2004;134: 9-21. doi: 10.1016/j.jneumeth.2003.10.009. pmid:15102499

Cohen MX. Analyzing neural time series data: theory and practice. Analyzing neural time series data: theory and practice. Cambridge, MA: MIT Press; 2014